The James Webb Space Telescope JWST observations have revolutionized extragalactic research particularly with the discovery of little red dots LRD which we propose are dustreddened broadline active galactic nuclei AGNs. Their unique vshape spectral feature observed through JWSTNIRCam challenges us to discern the relative contributions of the galaxy and AGN. We study a spectral energy distribution SED model for LRDs from restframe UV to infrared bands. We hypothesize that the incident radiation from an AGN characterized by a typical SED is embedded in an extended dusty medium with an extinction law similar to those seen in dense regions such as Orion Nebula or certain AGN environments. The UVoptical spectrum is described by dustattenuated AGN emission featuring a red optical continuum at lambdagt4000 A and a flat UV spectral shape established through a gray extinction curve at lambdalt3000 A due to the absence of smallsize grains. There is no need for additional stellar emission or AGN scattered light. In the infrared the SED is shaped by an extended dust and gas distribution gammalt1 rhopropto rgamma with a characteristic gas density of simeq 10103rm cm3 which allows relatively cool dust temperatures to dominate the radiation thereby shifting the energy peak from near to midinfrared bands. This model unlike the typical AGN hot torus models can produce an infrared SED flattening that is consistent with LRD observations through JWST MIRI. Such a density structure can arise from the coexistence of inflows and outflows during the early assembly of galactic nuclei. This might be the reason why LRDs emerge preferentially in the highredshift universe younger than one billion years.
['Astrophysics - Astrophysics of Galaxies']
Little Red Dots Rapidly Growing Black Holes Reddened by Extended Dusty Flows
{'Little Red Dots: Rapidly Growing Black Holes Reddened by Extended Dusty Flows': '<!-- image -->  \n1 Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China  \n2 Department of Astronomy, School of Physics, Peking University, Beijing 100871, China  \n- 3 Global Center for Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan\n- 4 Department of Physics, School of Advanced Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan', 'ABSTRACT': 'The James Webb Space Telescope (JWST) observations have revolutionized extragalactic research, particularly with the discovery of little red dots (LRD), which we propose are dust-reddened broad-line active galactic nuclei (AGNs). Their unique v-shape spectral feature observed through JWST/NIRCam challenges us to discern the relative contributions of the galaxy and AGN. We study a spectral energy distribution (SED) model for LRDs from rest-frame UV to infrared bands. We hypothesize that the incident radiation from an AGN, characterized by a typical SED, is embedded in an extended dusty medium with an extinction law similar to those seen in dense regions such as Orion Nebula or certain AGN environments. The UV-optical spectrum is described by dust-attenuated AGN emission, featuring a red optical continuum at λ > 4000 ˚ A and a flat UV spectral shape established through a gray extinction curve at λ < 3000 ˚ A, due to the absence of small-size grains. There is no need for additional stellar emission or AGN scattered light. In the infrared, the SED is shaped by an extended dust and gas distribution ( γ < 1; ρ ∝ r -γ ) with a characteristic gas density of ≃ 10 -10 3 cm -3 , which allows relatively cool dust temperatures to dominate the radiation, thereby shifting the energy peak from near- to mid-infrared bands. This model, unlike the typical AGN hot torus models, can produce an infrared SED flattening that is consistent with LRD observations through JWST MIRI. Such a density structure can arise from the coexistence of inflows and outflows during the early assembly of galactic nuclei. This might be the reason why LRDs emerge preferentially in the high-redshift universe younger than one billion years.  \nKeywords: Galaxy formation (595); High-redshift galaxies (734); Quasars (1319); Supermassive black holes (1663)', '1. INTRODUCTION': "The James Webb Space Telescope (JWST) has dramatically advanced the study of high-redshift ( z > 4) active galactic nuclei (AGNs) by unveiling low-luminosity AGN populations more representative than luminous quasars (e.g., Onoue et al. 2023). Among these discoveries, very red and compact objects, so-called little red dots (LRDs), have attracted considerable attention (e.g., Barro et al. 2023; Kocevski et al. 2023; Labbe et al. 2023; Matthee et al. 2023). Spectroscopic observations on photometrically identified LRDs have  \ninayoshi@pku.edu.cn  \nrevealed that more than 70% exhibit broad hydrogen Balmer lines (e.g., Greene et al. 2024), indicating the presence of massive black holes (BHs). The abundance of these AGNs in cosmic volume (Kokorev et al. 2023; Kocevski et al. 2024; Akins et al. 2024) is significantly higher than what is expected from ground-based quasar surveys (e.g., Matsuoka et al. 2018, 2023; Niida et al. 2020; He et al. 2023), and provides insights on the cosmological evolution of quasar luminosity functions (Li et al. 2023), the impact on cosmic reionization (e.g., Dayal et al. 2024; Madau et al. 2024), and the implications to rapid spins of these high-redshift BHs through the measurement of the radiation-to-mass conversion efficiency (Inayoshi & Ichikawa 2024).  \nDespite their significance, the properties of these newly identified LRDs remain puzzling, particularly concerning the origin of their characteristic v-shaped spectral energy distribution (SED) (e.g., Barro et al. 2023; Kocevski et al. 2023; Labbe et al. 2023). These SEDs feature a red continuum in the rest-frame optical spectrum at ≳ 4000 ˚ A and a blue excess in the restframe UV range. Observational studies using SED template fitting suggested that the red optical part could either be due to dust-reddened AGNs or starburst galaxies, while the blue UV component might represent light from unobscured galaxies or scattered light from embedded AGNs.  \nFurthermore, rest-frame near- to far-infrared observations offer additional insights. Data from the JADES, COSMOS-Web and PRIMER-COSMOS surveys show that most LRDs are not detected in the JWST/MIRI F1800W bands, and even for the detected sources, the MIRI SEDs tend to flatten from the rest-frame optical to the near-infrared (NIR) at 1-3 µ m (e.g., Williams et al. 2024; Akins et al. 2024; see also Wang et al. 2024). Those observations appear to contradict expectations that the SED should be very red throughout the optical to NIR spectrum due to re-emission from hot dust components heated by the AGNs. Consequently, these results have led to arguments disfavoring the scenario accounting for dust-obscured AGNs. While a certain fraction of LRDs show extremely red SEDs from NIR to mid-infrared (MIR) (Lyu et al. 2024; P'erez-Gonz'alez et al. 2024), the reasons for such diversity at λ ≳ 1 µ m are not well understood. Alternatively, if considering dust-obscured galaxies instead of AGNs, cold dust associated with starbursts would dominate the far-infrared (FIR) light. However, current ALMA observations provide no strong evidence to support this scenario (Labbe et al. 2023; Williams et al. 2024; Akins et al. 2024; Casey et al. 2024).  \nIt is worth noting that LRDs are extremely compact (effective radius of R eff < 100 pc at z ∼ 5; Labbe et al. 2023; Furtak et al. 2023) and almost point-like objects even with supreme spatial resolution by JWST. A small fraction of red sources with a v-shaped SED show extended morphology in images of shorter-wavelength filters, which may indicate the presence of an underlying host galaxy (Killi et al. 2023) or diffuse ionized gas. However, most LRDs do not have detectable extended components, allowing us to place upper limits on the mass and size of the galaxies. Comparing the BH mass estimated from broad-line emission spectroscopy to the upper limits on galaxy mass, the BH-to-galaxy mass ratio in LRDs (e.g., Pacucci et al. 2023; Kokorev et al. 2023; Greene et al. 2024) is substantially higher  \nthan the empirical value observed in the nearby universe (e.g., Kormendy & Ho 2013). This implies the existence of an evolutionary stage where growing BHs occupy ∼ a few -10% of the total BH+galaxy system (for theoretical results see e.g., Inayoshi et al. 2022; Hu et al. 2022; Scoggins et al. 2023).  \nIn this paper, we study a SED model for LRDs from rest-frame UV to infrared bands. We hypothesize that the incident radiation from an AGN, characterized by a typical SED, is embedded in an extended dusty medium with an extinction law similar to those seen in dense regions such as Orion Nebula or certain AGN environments. The UV-optical spectrum is described by dustattenuated AGN emission, composed of a red optical continuum and a UV excess established through gray extinction. There is no need for additional stellar emission or AGN scattered light. In the infrared, the SED is shaped by an extended dust and gas distribution, which allows relatively cool dust temperatures to dominate the radiation, thereby shifting the energy peak from near- to mid-infrared bands. This model can produce an infrared SED flattening that is consistent with LRD observations through JWST MIRI. We finally provide discussion and open questions rising from our conclusions regarding (i) the temporal emergence of LRDs at z > 4, (ii) the universality of v-shaped SEDs, and (iii) the cosmological evolution of early BH populations.", "2. THE 'V'-SHAPE SED IN REST-FRAME UV/OPTICAL": "It is evident that dust extinction curves in AGN significantly deviate from those observed in the Milky Way (e.g., Maiolino et al. 2001b,a, 2004; Gaskell et al. 2004; Li 2008; Xie et al. 2017), showing a characteristic flat wavelength dependence at λ < 4000 ˚ A. This feature is believed to stem from a deficiency of small-size grains within the AGN environment, which can be attributed to dust destruction mechanisms such as thermal sublimation and sputtering (e.g., Laor & Draine 1993; Tazaki & Ichikawa 2020), and grain charging (e.g., Draine & Salpeter 1979). Recently, Tazaki et al. (2020) have applied these dust destruction processes specifically at the polar regions of 1 -10 pc directly irradiated by the AGN (e.g., Honig & Kishimoto 2017), and show the critical role of Coulomb explosions in flattening the extinction curve.  \nFigure 1 presents seven different dust extinction laws, normalized at the visual extinction level at 5500 ˚ A ( A λ /A V ): the Small Magellanic Cloud (SMC, Gordon et al. 2003), starburst galaxies (Calzetti et al. 2000), the interstellar dust in Milky Way with R V = 3 . 1 and 5.0 (Cardelli et al. 1989), high-redshift galaxies  \nFigure 2. The flux densities of dust-reddened AGNs normalized at 4000 ˚ A (solid curves), where the incident quasar SED is adopted from a SED template model with a removal of host galaxy and dust emission contributions (dashed curves, Temple et al. 2021). The six different extinction laws shown in Figure 1 are applied by setting a visual extinction of A V = 3 mag. The Markov's extinction curve is not used because of the lack of full wavelength coverage. The SED shaped with the Orion Nebula extinction law reproduces the characteristic v-shape SED of LRDs. The averaged photometric data of LRDs at z ∼ 5 . 5 and 7 . 5 are overlaid for comparison (Barro et al. 2023).  \n<!-- image -->  \nFigure 1. Extinction curves as a function of wavelength normalized by the visual extinction at 5500 ˚ A: the Small Magellanic Cloud (SMC, Gordon et al. 2003), starburst galaxies (Calzetti et al. 2000), the interstellar dust in Milky Way with R V = 3 . 1 and 5.0 (Cardelli et al. 1989), high-redshift galaxies at 6 < z < 12 (Markov et al. 2024), the Orion Nebula (Baldwin et al. 1991), and composite AGN spectra (Gaskell et al. 2004). Extinction curves in the absence of small-size dust grains, as observed in the Orion Nebula, are significantly flattened at λ < 2500 ˚ A, unlike the other models commonly used in galaxy SED fitting methodology (e.g. Calzetii's law in Labbe et al. 2023).  \n<!-- image -->  \nat 6 < z < 12 (Markov et al. 2024), the Orion Nebula (Baldwin et al. 1991), and composite AGN spectra (Gaskell et al. 2004). The SMC and Calzetti's laws are commonly used in galaxy SED fitting methodology (e.g., Carnall et al. 2018; Boquien et al. 2019). These models suggest that the extinction level nearly monotonically increases to shorter wavelengths. The extinction curve for the Milky Way interstellar dust also shows a similar behavior with a bump at λ ∼ 2175 ˚ A, which is considered to originate from small graphite grains and/or polycyclic aromatic hydrocarbon nanoparticles (Draine &Lee 1984). In contrast, the extinction law measured in the Orion Nebula, one of the brightest nebulae powered by the central massive star forming region, exhibits a curve that is significantly flattened in the far-UV wavelength range and shows the absence of the 2175 ˚ A bump. The relatively grey extinction curve can be explained by the deficit of small-size dust grains. Importantly, similar extinction laws are also observed in those of AGNs (Gaskell et al. 2004; see also Maiolino et al. 2001a, Xie & Ho 2022), which might preferentially remove small grains with sizes of a ≲ 0 . 1 µ m. Additionally, JWST NIRSpec observations indicate a flattening trend of dust extinction in star-forming galaxies at 6 < z < 11, potentially due to the formation of large-size dust grains within energetic supernova ejecta (Markov et al. 2024).  \nIn Figure 2, we show the attenuated flux using the five extinction laws. We here adopt the composite SED of lowz quasars from Temple et al. (2021) as the incident flux. This model provides empirical SED templates for unobscured quasars with E ( B -V ) < 0 . 3 by using a total of ∼ 1 . 5 × 10 4 SDSS quasars with redshifts 0 < z < 5. The SED covers from the rest-frame UV to NIR up to 3 µ m. In these templates, the contribution from starlight from the host galaxy is quantified using large quasar samples and subtracted to isolate the pure quasar SED. For this work, we generate the SED template at z ≃ 5, which excludes both the host galaxy and dust contribution. In contrast, the widely-used quasar SED template provided by Vanden Berk et al. (2001) is known to retain significant galaxy contributions, leading to a substantially redder spectrum at λ ≳ 5000 ˚ A compared to those in Temple et al. (2021). Although this discrepancy between the two SED models does not impact our analysis of the rest-frame UV-to-optical characteristics of LRDs, it introduces a serious inconsistency with rest-frame NIR observations due to inadequate subtraction of host galaxy contributions (e.g., MIRI photometry; P'erez-Gonz'alez et al. 2024). Therefore, for the purpose of accurately modeling LRD SEDs in the absence  \nof host galaxy effects, we intentionally use the model from Temple et al. (2021) as the incident flux.  \nFor all the cases, we assume a visual extinction A V = 3 mag and normalize the attenuated fluxes at λ = 4000 ˚ A. As reference, we overlay photometric SED data obtained from stacking LRD samples at z ∼ 5 . 5 (orange) and z = 7 . 5 (blue) (Barro et al. 2023). The attenuated flux density reproduces a red continuum component at λ ≳ 3000 -4000 ˚ A with a slope of β opt ≳ 0, where the slope is defined as β ≡ dln f λ / dln λ . However, the SED shape significantly differs in the rest-UV part, depending on the assumed extinction curves. Models with extinction laws observed in the Orion Nebula and AGN composite spectra flatten the SED with a UV slope of β UV ∼ -2. In contrast, other models result in strongly attenuated flux that fails to maintain a UV excess. In the two models, the incident UV flux is reduced by two orders of magnitude, i.e., 1% of the incident flux.  \nThe application of the extinction law derived from AGN composite spectra (Gaskell et al. 2004) was tested on a LRD (ID: J0647-1045) identified at z = 4 . 5321 by JWST NIRSpec (Killi et al. 2023). For this particular object, the inferred extinction curve matches that proposed by Gaskell et al. (2004). Our findings in this work further demonstrate that such gray extinction curves can explain the characteristic v-shaped SEDs of LRDs, as evidenced by accumulated large samples, without considering multiple components such as host galaxies or scattered AGN light.  \nPrevious observational studies have proposed that the observed excess in rest-frame UV ( β UV ∼ -2) may be attributed to the scattering of unobscured quasar light (Kocevski et al. 2023; Labbe et al. 2023; Greene et al. 2024). According to this scenario, the fraction of the scattered flux relative to the primary component heavily depends on both the covering factor of the scattering medium and the viewing angle, estimated to be about 1-3% to account for the observed UV characteristics of LRDs. However, this scenario does not account for the scattering fraction uniformly observed across most LRDs. Relying on such a coincidental process to explain the universal characteristics of an abundant AGN population is challenging 1 . In contrast, our model, which applies the dust opacity law observed in  \nFigure 3. The observed SEDs of LRDs across various redshifts of 4 ≤ z ≤ 7, where the intrinsic AGN luminosity is set to L bol = 10 46 erg s -1 . We adopt the composite SED template from Temple et al. (2021), and assume the Orion-like extinction law with A V = 3 mag. The 5 σ pointsource imaging depths of JWST survey programs are overlaid: COSMOS-Web (gray, Casey et al. 2023), CEERS (purple, Finkelstein et al. 2023), and JADES-Medium/Deep (blue solid and dashed, Eisenstein et al. 2023).  \n<!-- image -->  \nnearby or low-redshift AGNs, offers a more robust scenario for the LRD characteristics (see also discussion in Section 4.3), and furthermore bring us new insights on the AGN structure within parsec scales.  \nFigure 3 presents the SEDs of our LRD model at various redshifts from z = 4 to z = 7. For these models, we assume a bolometric luminosity of L bol = 10 46 erg s -1 as the incident radiation from the AGN, and calculate the monochromatic luminosity at rest-frame 1450 ˚ A using a bolometric correction factor of f bol , UV = 4 . 4 (Richards et al. 2006). Each case employs the OrionNebula attenuation law with a visual extinction A V = 3 mag, as shown by the red curve in Figure 2. For comparison, we overlay the 5 σ point-source imaging depths for each NIRCam broadband filter used in several JWST survey programs. The depths achieved by the CEERS and JADES surveys are sufficient to detect the characteristic SED shape of LRDs from the restframe UV to optical wavelengths. This is consistent with the observational fact that LRDs identified in these deep surveys exhibit bolometric luminosities similar to our models, L bol ∼ 10 45 -46 erg s -1 (e.g., Kokorev et al. 2024). However, the COSMOS-Web survey, which covers wider areas ( ∼ 0 . 5 deg 2 , approximately ∼ 20 times larger than other surveys with ∼ 100 arcmin 2 ) with relatively shallower depths (as shown in Figure 3), requires higher bolometric luminosities for the UV component to be detectable and to capture the characteristic SED  \nshape. Our SED models satisfy the selection criteria for LRDs; β UV < -0 . 37 and β opt > 0 (Kocevski et al. 2024; see also Greene et al. 2024).", '3. REPROCESSED INFRARED RADIATION FROM OBSCURING DUST': "Next, we analyze the SED of LRDs in the IR band, focusing on the re-emission from dust grains that obscure the AGN light from UV to optical wavelengths. Since UV photons emitted from the AGN heat the obscuring material, the dust temperature rises to a sublimation point of T sub ≃ 1500 K, and decreases with distance from the AGN. In scenarios where dust and gas are centrally concentrated, as in dusty tori (e.g., Ramos Almeida & Ricci 2017), the hot dust just behind the sublimation radius significantly contributes to the dust-grain mass budget and thus preferentially re-emits thermal heat into the NIR bands, λ ≃ 2 . 9 µ m ( T dust / 10 3 K) -1 . As a result, the NIR flux density can be ≳ 2 dex higher than the attenuated AGN emission at the optical band, resulting in an extremely red continuum from the optical to infrared regimes, as observed in dust-obscured AGN at lower redshifts (e.g., Alonso-Herrero et al. 2006; Dey et al. 2008; Ramos Almeida et al. 2011; Ichikawa et al. 2015, 2017; Toba et al. 2017; Lyu & Rieke 2022). However, a good fraction of high-redshift LRDs observed with MIRI broadband filters (e.g., F770W and F1800W) do not show such red color but rather exhibit a relatively flat SED from the rest-frame optical to the NIR wavelengths (Williams et al. 2024; Wang et al. 2024). This trend is also seen in stacking images of multiple LRDs (Akins et al. 2024). Nevertheless, there remain LRDs that show extremely red SEDs from NIR to MIR (Lyu et al. 2024; P'erez-Gonz'alez et al. 2024).  \nThe shape of the IR SED, reprocessed by dust grains heated by an AGN, depends on the radial distribution of the dust (e.g., Pier & Krolik 1992; Nenkova et al. 2008; Nikutta et al. 2021a). When the obscuring material does not form a concentrated density structure or is dispersed by bipolar outflows, the peak energy shifts to the MIR since the bulk of the re-emission comes from relatively cooler dust located farther from the dust sublimation radius (e.g., Barvainis 1987). NIR interferometry observations of nearby AGNs in recent decades (e.g., Kishimoto et al. 2011a,b) have found that the NIR emitting region appears to have a radius nearly twice as large as the expected sublimation radius, as traced by the NIR dust reverberation mapping observations (Suganuma et al. 2006; Koshida et al. 2014). This suggests that the dust density profile has a slightly lower level of concentration, producing a larger emissivity size than the sublimation  \nradius (see also Nikutta et al. 2021b). Additionally, MIR interferometry observations resolving the central ∼ pc scale also showed that the significant fraction of MIR dust emission originates from the polar region, suggesting the presence of bipolar dusty outflows (e.g., Jaffe et al. 2004; Honig et al. 2012, 2013; Honig 2019; Tristram et al. 2014). According to radiation transfer calculations that explain the IR SED of AGNs, including dust emission towards the polar region, the gas and dust density profiles with a low level of concentration ( ρ ∝ r -γ ; -0 . 5 < γ < 0 . 5) can reproduce the observed MIR bump along with a relatively flat NIR SED at λ ≃ 3 -5 µ m (Honig & Kishimoto 2017).  \nMotivated by these AGN IR observations and models, we focus not only on a compact dust structure such as a dust torus, but rather a more extended dust structure described by a radial density profile of ρ = ρ 0 ( r/r in ) -γ , where ρ 0 is the gas density at the sublimation radius r in , and the power-law index ranges 0 < γ < 2. With the underlying density distribution, the temperature of dust grains, T dust ( r ), can be determined by balancing dust heating from irradiation with cooling through thermal emission as  \n∫ L ν e -τ ν ( r ) 4 πr 2 κ abs ν dν = 4 π ∫ κ abs ν B ν [ T dust ( r )] dν, (1)  \nwhere L ν is the luminosity density of the central AGN at a frequency of ν , κ abs ν is the absorption opacity of dust, B ν ( T dust ) is the Planck function with a dust temperature T dust ( r ) at a distance of r from the center, and τ ν is the optical depth defined by  \nτ ν ( r ) = ∫ r r in ρκ ext ν dr ' , (2)  \nwhere κ ext ν is the sum of the absorption and scattering opacity.  \nWith a dust temperature profile, the luminosity density of dust re-emission is calculated by  \nL ν, dust = 4 π Ω ∫ r out r in κ abs ν B ν [ T dust ( r )] ρr 2 dr, (3)  \nwhere 0 < Ω ≤ 4 π is the solid angle covered by dust and r out is determined by a constraint on the column density so that the visual extinction reaches A V ≃ 3 mag to reproduce the characteristic red continuum spectra of LRDs in the rest-frame optical bands (see Section 2). In our calculation, we assume a full-sky covering case (i.e., C ≡ Ω / 4 π = 1), while this value will be determined when the relative cosmic abundance between LRDs and unobscured AGNs and the dependence of C on the AGN properties such as the Eddington ratio are better understood for LRDs (e.g., Ichikawa et al. 2019; Ricci et al.  \nFigure 4 presents multi-wavelength SEDs of a LRD from the rest-frame UV to far infrared with various combinations of ( n 0 , γ ). This wavelength coverage corresponds to the observed wavelengths frame from JWST NIRCam/MIRI to ALMA bands for sources at z = 6. For comparison, we overlay the stacking photometric data of observed LRDs reported by the JADES survey (thick curve, P'erez-Gonz'alez et al. 2024) and the COSMOS-Web survey (circles, Akins et al. 2024). Furthermore, the upper bounds of the mid- to far-IR flux densities due to non-detection in SCUBA and ALMA bands are shown (triangles, Akins et al. 2024; see also Labbe et al. 2023).  \n<!-- image -->  \nFigure 4. The flux densities of LRDs normalized at 4000 ˚ A for different values of n 0 and γ , where the incident SED is the same as in Figure 2. The stacked photometric data of LRDs are taken from the JADES survey (thick curve, P'erez-Gonz'alez et al. 2024) and the COSMOS-Web survey (circles Akins et al. 2024), and the upper bounds of the mid- to far-IR flux densities due to non-detection in SCUBA and ALMA bands are taken from Akins et al. (2024). Left panel: dependence of the SED on the density n 0 . The highest density case of n 0 = 10 3 cm -3 results in overly bright rest-frame NIR, violating the MIRI F770W constraints. Lower density cases reproduce SEDs consistent with the observed data points, while the lowest density of n 0 = 1 cm -3 yields FIR SED peak inconsistent with the upper bounds from the ALMA observations. Right panel: Variation of SEDs with different density slopes γ . As the density is less concentrated ( γ < 1), the energy peak shifts from the NIR to the MIR regime.  \n<!-- image -->  \n2017, 2023). Note that the equation (3) can be applied in the optically thin limit of dust re-emission in the dusty surrounding medium. This approximation is valid for a wide range of parameter combinations of the density value ( n 0 ≡ ρ 0 /m p ) and slope ( γ ) of the distribution we explore in our SED modeling for LRDs (see Section 4.1).  \nIn addition to a diluted density distribution, we consider the feedback effect on the dusty medium surrounding the central AGN (e.g., Chan & Krolik 2016; Namekata & Umemura 2016; Kudoh et al. 2023; Soliman & Hopkins 2023). Radiation pressure expels dust grains between the sublimation radius and the dust photosphere to UV radiation into a thin shell that accumulates at r = r ph , the radius of the photosphere where the optical depth measured from the center reaches unity. Within the thin layer at r = r ph , the optical depth transitions from 0 to 1, and the dust temperature varies depending on the detailed density distribution. For simplicity, we assume a constant temperature within this layer, corresponding to the inner edge of the thin shell. This temperature can be determined using Equation (1) with r = r ph and τ ν ( r ) = 0, setting an upper bound on the temperature and emission from the thin shell. In our model, feedback effects are incorporated by substituting r in with r ph in Equation (3). While detailed radiation hydrodynamic simulations can pinpoint the exact location of the photosphere, we evaluate r ph using the original unperturbed density profile, because the optical  \ndepth value within a radius of r remains unchanged if mass is conserved. As a result of the feedback, hot dust just behind the sublimation radius r in is pushed to larger radii and becomes cooler, shifting the NIR emission to MIR. Despite these change within r ph , the temperature profile and emission of the region outside r ph remain unaffected.  \nIn the left panel of Figure 4, we show the dependence of the SED shape on the gas density (1 ≤ n 0 / cm -3 ≤ 10 3 ) at the sublimation radius, while the slope of its radial profile is fixed to γ = 0. When the density is as high as n 0 = 10 3 cm -3 , the outer edge of the obscuring material is located at r out ≳ 30 r in , where the dust temperature ranges from T dust ≃ 1500 K to ∼ 200 K. Therefore, the reprocessed IR spectrum begins to rise  \nFigure 5. Opt-NIR slope β opt / NIR as the function of density power-law index γ and gas density n 0 . High density cases lead to bright rest-frame NIR emission and thus a steeper slope. The solid curve shows the boundary of β opt / NIR = -1, below which the NIR SED in our model is consistent wit the requirement from JWST/MIRI observations for LRDs. The dashed curves present constant 30 ≤ T out / K ≤ 120 and the temperature of Cosmic Microwave Background (CMB) at z = 6. Below T out ≃ 30 K, the SED becomes inconsistent with the non-detection of LRDs with ALMA observations. The allowed parameter space is approximated as 10 2 . 9 γ +0 . 83 ≲ n 0 / cm -3 ≲ 10 3 .  \n<!-- image -->  \nfrom λ ≃ 1 µ m and exhibits a near- to mid-infrared bump that covers λ ≃ 2 -15 µ m. In this case, the bright NIR bump makes the SED shape substantially redder from λ ≃ 6000 ˚ A to ∼ 2 µ m, consistently observed in AGNs with hot dust tori (Polletta et al. 2007). However, this extremely red SED color does not match the JWST MIRI photometric data for LRDs, which suggest a modest red color or even a flat spectral shape. We note that such a contradiction with the MIRI data appears as long as n 0 ≳ 10 3 cm -3 .  \nOn the other hand, as the density n 0 decreases and the outer radius expands to reach a sufficient attenuation level, the peak energy of the IR SED decreases because relatively cooler dust grains dominate both the mass and IR emissivity. The energy transfer from NIR to MIR significantly suppresses the NIR bump and eases the mismatch between our SED model and the JWST MIRI photometric data. For the cases with n 0 = 10 and 100 cm -3 , the SED shapes from the optical band to λ = 2 -3 µ m agree to the observed spectra of LRDs (or their stacked photometry based on the F444W and F770W filters for z ≃ 6 LRDs), and also the peak energy is shifted to the FIR band at λ ∼ 20 -50 µ m. As the density decreases to n 0 = 1 cm -3 , the IR SED peaks at  \nλ ∼ 100 µ m and exceeds the upper bounds obtained by ALMA observations.  \nWe also observe that our SED model with γ = 0 and n 0 < 10 3 cm -3 aligns with the upper bounds of the F1800W photometry. It is important to note that the F1800W upper limit is derived from stacking images covering less than 20% of the LRD samples observed with the F770W MIRI filter (Akins et al. 2024). Further investigations into the rest-frame MIR band, such as through F1800W, are crucial for enhancing our understanding of the nature of LRDs (see the PRIMER survey; Dunlop et al. 2021).  \nThe right panel of Figure 4 presents cases for a centrally concentrated density distribution with γ = 2 and γ = 1, alongside the case of γ = 0 and n 0 = 100 cm -3 (see the left panel). At the sublimation radius, the densities are set to n 0 = 4 × 10 4 cm -3 and 4 × 10 3 cm -3 for γ = 2 and 1, respectively, ensuring sufficient extinction ( A V ≃ 3 mag). When the density is centrally concentrated with γ = 1 -2, a significant fraction of the obscuring material is efficiently heated by the AGN, and thus the hot dust just behind the sublimation radius is responsible for the IR emissivity. As a result, the spectrum appears redder in the NIR wavelengths, similar to the high-density case shown in the left panel.", '4.1. Physical properties of gas obscuring LRDs': "In this section, we explore the parameter regions ( n 0 , γ ) that satisfy the conditions required for capturing the characteristic SED shape of LRDs. We consider the following three key conditions.  \n- 1. Dust attenuation: To reproduce the rest-frame UV-to-optical SED of LRDs, the dust attenuation is set to A V ≃ 3 mag. This level of attenuation corresponds to a hydrogen column density of N H ≃ 7 . 5 × 10 22 cm -2 ( Z/ 0 . 1 Z ⊙ ) -1 (Baldwin et al. 1991).\n- 2. Optical-NIR SED slope: The spectral slope between the rest-frame optical and NIR regions needs to be consistent with the JWST NIRCam/MIRI photometry, showing that the F444W -F770W colors for z ≃ 5 -6 LRDs are ≃ 0 . 5 mag in stacking images. To quantify the optical-to-NIR slope, we analyze flux densities at λ = 0 . 65 and 1 . 2 µ mfrom our SED model. An crucial consideration is the influence of H α line emission, which can cause an excess in the NIRCam filter-convolved flux density at λ ≃ 0 . 65 µ m. For LRD samples identified by slitless spectroscopy (Matthee et al. 2023), we adopt an average rest-  \nalent width (EW) of EW 0 , H α = 500 ˚ A. The prominent line emission enhances the filterconvolved flux densities within a filter with a bandwidth of ∆ λ filter by a factor of ℓ ≡ [EW 0 , H α (1 + z ) + ∆ λ filter ] / ∆ λ filter . Assuming that the H α line emission falls within the F444W, this excess factor of ℓ ≃ 1 . 35 applies for a z = 6, corresponding to an additional 0 . 33 mag. Therefore, we define the optical-to-NIR slope as  \nβ opt / NIR = log 10 [ f ν 2 / ( ℓf ν 1 )] log 10 ( λ 2 /λ 1 ) -2 , (4)  \nwhere λ 1 = 0 . 65 µ m, λ 2 = 1 . 2 µ m to measure the continuum flux f ν , and ν 1(2) = c/λ 1(2) . Here we impose β opt / NIR < -1 . 0 for this study.  \n- 3. ALMA non-detection: Given the non-detection of LRDs in ALMA observations, we impose a lower limit on the dust energy that re-emits the AGN light. This requires that the minimum dust temperature, which is analytically calculated in our SED model, exceeds a temperature threshold of T min = 30 K.  \nIn the analysis described in Section 3, we find that the SED peaks in the NIR regime and fails to meet the slope criterion of β opt / NIR < -1 when a centrally concentrated density profile with γ > 1 is assumed (see the right panel of Figure 4). We thus focus on the parameter region with γ ≤ 1 in the following discussion.  \nFigure 5 summarizes the characteristic properties of the SED across various combinations of n 0 and γ . The color contours indicate the optical-to-NIR slope as defined by Equation (4). Employing the criterion of β opt / NIR = -1, the parameter space is divided into two regions. In the high density region ( n 0 ≳ 10 3 cm -3 ), the NIR SED is influenced by re-emission from hot dust grains heated by the AGN, resulting in a steeper opticalto-NIR slope that is inconsistent with the observational requirement. We also present curves showing a constant dust temperature T out at the outer edge of the obscuring material at r = r out . In cases with lower density, the outer radius further extends to achieve A V ≃ 3 mag and T out also decreases. By imposing T out ≥ T min = 30 K, we ensure that the FIR part of the SED aligns with non-detection of LRDs with ALMA.  \nIn summary, the region between the two curves of β opt / NIR = -1 and T out = 30 K defines the allowed parameter space of ( n 0 , γ ). Here, the density range is approximated as  \n10 2 . 9 γ +0 . 83 ≲ n 0 / (cm -3 ) ≲ 10 3 , (5)  \nwhich is valid γ < 0 . 8. We also note that the density values are the lower side of the typical electron den-  \nity of the low redshift narrow line region (NLR) with n e ∼ 10 -10 6 cm -3 (e.g., Netzer 1990; Nesvadba et al. 2006; Liu et al. 2013; Harrison et al. 2014). Recently, Integral Field Unit (IFU) observations that spatially resolves NLRs of low-redshift AGNs show that the electron density at the inner ∼ 100 pc scales is found to be n e ∼ 10 2 -10 3 cm -3 , and decreases toward larger radii (Kakkad et al. 2018). For high-redshift star forming galaxies at z ∼ 6 observed by JWST/NIRSpec, the electron density in the nebular region is estimated as n e ≳ 10 2 cm -3 (e.g. Isobe et al. 2023; Abdurro'uf et al. 2024).  \nFinally, Figure 6 shows the compiled SED models that agree to the MIRI (F770W and F1800W) and ALMA constraints. We consider a LRD at z = 6 with a bolometric luminosity of L bol = 5 × 10 45 erg s -1 . The color code presents variations of gas density n 0 , while maintaining a constant radial distribution index γ = 0. For comparison, we also show the SED templates of AGN hot dust tori (Torus model), Apr 220, and Mrk 231 provided by (light blue curves Polletta et al. 2007), for which the same bolometric luminosity is assumed. The SED shapes of these compared templates violate observational constraints either at the NIR or FIR regimes.  \nWith the median value of the density in the log space ( n 0 ∼ 100 cm -3 ), the dust mass enclosed within r ≤ r out is estimated as M dust ≃ 1 . 2 (9 . 6) × 10 4 M ⊙ for γ = 0 and 0 . 2, respectively, assuming a dust-togas mass ratio of 0 . 001( Z/ 0 . 1 Z ⊙ ). These results agree well to the mean value of the dust mass ⟨ M dust ⟩ ≃ 1 . 6 +4 . 8 -0 . 9 × 10 4 M ⊙ estimated using NIRCam images of LRDs (Casey et al. 2024). Adopting the dust-to-stellar mass ratio of f ⋆, dust ∼ 10 -4 -10 -3 in a galaxy younger than one billion years (e.g., Valiante et al. 2009), the expected mass of stars that form these dust grains would be the order of ∼ 10 7 -9 M ⊙ . This estimate can be used in comparison to other independent stellar-mass measurements, e.g., by image decompositions and SED fitting.", '4.2. Host galaxy contributions?': 'Previous observations of LRDs have provided limited information about their host galaxies. While some LRDs exhibit extended morphology in short-wavelength images, achieving a robust confirmation of host galaxy features remains challenging. Detailed measurements of JWST/PRISM spectra have successfully identified a continuum break feature at λ ≃ 4000 ˚ A in two LRDs: MSAID13123 at z = 7 . 04 in Greene et al. (2024) and RUBIES-BLAGN-1 at z = 3 . 1 in Wang et al. (2024). These findings suggest the presence of an underlying host galaxy. However, if stellar light accounted for the  \nFigure 6. Observed SEDs of LRDs in our model assuming an AGN bolometric luminosity of L bol = 5 × 10 45 erg s -1 at z = 6. The colored curves presents variations in SEDs with gas density n 0 , ranging from 150 cm -3 (dark color) to 10 cm -3 (light color), while maintaining a constant radial distribution index γ = 0. The stacked photometric data for LRDs and upper bounds talen from Akins et al. (2024) are overlaid. For comparison, we also show the SED templates of AGN hot dust tori (Torus model), Apr 220, and Mrk 231 provided by (light blue curves, Polletta et al. 2007), for which the same bolometric luminosity is assumed. At the bottom of the plot, we show the wavelength coverage of related telescopes.  \n<!-- image -->  \ncontinuum flux across the entire optical band, the required stellar mass would approach M ⋆ ≳ 10 11 M ⊙ , comparable to that of the Milky Way. When combining such a high mass and their high abundance in the cosmic volume, this scenario implies an extraordinary stellar mass density that contradicts the standard cosmological framework (see also Inayoshi & Ichikawa 2024; Akins et al. 2024). As a possible solution, Wang et al. (2024) suggest that the stellar contribution on the SED is limited to the short-wavelength range at λ ≲ 4000 ˚ A and the AGN is responsible for the red optical continuum at the longer wavelengths. Under this assumption, the SED fitting indicates a minimal stellar mass of M ⋆ ∼ 10 9 M ⊙ , resolving the tension with the cosmological prediction.  \nHere, we emphasize that our approach does not rule out the existence of host galaxies in LRDs, but demonstrates that the SED model focusing only on AGNs successfully captures the characteristic spectral properties of LRDs: (1) the v-shaped SED in rest-frame UV to optical bands and (2) modest NIR emission. In our  \nAGN-only spectral model, as a natural consequence, the red optical continuum is attributed to the attenuated AGN light. This interpretation is indeed consistent with observations that the observed continuum and H α line fluxes for LRDs follow the empirical correlation observed in nearby AGNs (Greene & Ho 2005). In other words, if stellar light dominates the observed continuum flux of LRDs, the EW of broad H α emission line would be unusually high, compared to normal AGNs. This consistency is crucial, as it strengthens the interpretation of BH mass measurements using the single-epoch method.', '4.3. Why do LRDs emerge at high redshifts of z > 4 with a universal SED shape?': 'There is a critical question: Why do LRDs emerge at z > 4 with a uniform extinction level of A V = 3 mag? There is no clear answer to this question at this time, but some speculations are given below.  \nThe key assumption in our model is the relaxation of density concentration around the AGN. This idea implies that LRDs form in the early stages of galaxy and  \nAGN formation, characterized by dynamically evolving systems (e.g., clumpy gaseous structures, irregularly shaped and misaligned galactic disks, significant perturbations from mergers) before a well-organized dust torus structure is established (e.g., Ceverino et al. 2010, 2015; Hopkins et al. 2024). In the nascent stages of galaxy formation, the unsettled nature of the gas distribution permits a variety of attenuation levels and column densities. At higher redshifts, if the column density is typically higher and spans a broader range, numerous objects will exhibit A V ≫ 1. However, in the context of our proposed model, only those with moderate extinction levels with A V ∼ 1 -3 would reveal the characteristic SED shape, as seen in observed LRD samples (e.g., Barro et al. 2023; Matthee et al. 2023; Greene et al. 2024; Kocevski et al. 2024). With A V ≳ 5, which corresponds to A 3000 = 1 . 6 A V ≳ 8 at λ = 3000 ˚ A with the Orion-Nebula extinction curve, it is challenging to observe the rest-frame UV part of LRDs even with the JWST sensitivity. In this case, these objects appear to have only a red continuum without an UV excess, which are categorized as hot-dust-obscured AGNs observed at the lower-redshift universe (Alonso-Herrero et al. 2006; Ramos Almeida et al. 2011; Ichikawa et al. 2014; Hickox & Alexander 2018; Lyu & Rieke 2022). Our scenario seems more relevant in earlier epochs before a robust dust torus structure forms in its galactic nucleus. When a centrally concentrated dust torus is established, as seen in nearby AGNs, there is a dichotomy between type 1 and type 2 AGN populations due to the geometric arrangement (e.g., Ramos Almeida & Ricci 2017; Hickox & Alexander 2018) rather than a modest variation due to A V .  \nAnother critical element in our model is the flattening of dust extinction curves at shorter wavelengths, λ ≲ 3000 ˚ A. This characteristic may help explain why dust-reddened objects with v-shaped SEDs predominantly appear at high redshifts. In the early stages of protogalaxies, particularly those younger than one billion years ( z ≳ 6 -7), dust production is primarily driven by massive stars that expel heavy elements into the interstellar medium (ISM) through core-collapse supernovae (Todini & Ferrara 2001; Valiante et al. 2009). Due to effective dust destruction mechanisms that target smaller grains, the dust-size distribution tends to favor larger particles (Nozawa et al. 2006; Asano et al. 2013). These larger grains, once released into the ISM, undergo shattering processes and break into smaller fragments in a typical timescale of ≳ 0 . 5 -1 Gyr. As a result, the extinction curve at shorter wavelengths begins to evolve towards those observed in the ISM of the Milky Way. The temporal evolution in dust characteristics would be  \ncritical for shaping the unique spectral features of LRDs only at high redshifts.', '4.4. Implications to the cosmological BH evolution': 'In this paper, we propose that LRDs originate from rapidly growing BH embedded within dusty flows with an extended density distribution. In this framework, we reach two intriguing conclusions regarding the cosmological BH evolution.  \nFirst, the luminosity of LRDs is primarily powered by AGNs that host accreting BHs. These LRDs are more abundant than X-ray selected AGNs, suggesting an increase in the cosmic growth rate of BHs beyond z ≃ 6 (Kokorev et al. 2024; Kocevski et al. 2024; Akins et al. 2024). The BH accretion rate density derived from the bolometric luminosity functions of LRDs is significantly higher than that estimated from other AGN surveys. To reconcile this high value with the BH mass density at z ∼ 5, a radiative efficiency of ≳ 20% - twice the canonical 10% value - is required with a 3 σ confidence level (Inayoshi & Ichikawa 2024). This high radiative efficiency implies that the BHs are rapidly spinning in the LRD stage through violent mass accretion during the early galaxy assembly.  \nSecond, when the extinction level is accurately determined from the reddened optical continuum shape, the AGN continuum and broad-line luminosity can be more precisely estimated. These luminosity values are crucial for measuring BH mass using the single-epoch method. For instance, Matthee et al. (2023) calculated BH masses for LRDs using the H α luminosity as an indicator of AGN continuum luminosity, but without accounting for dust reddening correction. According to the empirical relationship in the single-epoch method, we derive that the Eddington ratio follows λ Edd ∝ L 0 . 394 H α (see Greene & Ho 2005; Matthee et al. 2023). Therefore, by applying a dust correction, the Eddington ratio increases by a factor of L ∼ 10 0 . 14 A V ∼ 2 . 6 for A V = 3. With this correction, the Eddington ratios of BHs reported by Matthee et al. (2023) approach unity, suggesting that a significant fraction of LRDs may be super-Eddington accretors. This result highlights the importance of accurate dust correction for high-redshift AGNs, particularly LRDs, to understand the true nature of their accretion dynamics (Lupi et al. 2024; see also Du et al. 2014).', '5. CONCLUSIONS': 'JWST observations have uncovered the presence of very compact and red objects at the high-redshift universe, referred to as little red dots (LRDs). Spectroscopic observations have confirmed the broad components of emission lines and the presence of AGNs  \nwithin LRDs. Despite the breakthrough discovery, the characteristic v-shape spectral feature observed through JWST/NIRCam challenges us to understand the origin of these LRDs and the contribution from AGNs powered by accreting BHs.  \nIn this work, we propose a SED model for LRDs spanning from rest-frame UV to infrared bands. We hypothesize that the incident radiation from an AGN, characterized by a typical SED, is embedded in an extended dusty medium, which has an extinction law similar to those seen in dense regions such as Orion Nebula or certain AGN environments. The UV-optical spectrum is described by dust-attenuated AGN emission with A V ≃ 3 mag, featuring a red optical continuum at λ > 4000 ˚ A and a flat UV spectral shape established through a gray extinction curve at λ < 3000 ˚ A, due to the absence of small-size grains. In the infrared, the SED is shaped by an extended dust and gas distribution ( γ < 1; ρ ∝ r -γ ) with a characteristic gas density of ≃ 10 -10 3 cm -3 at the dust sublimation radius, which allows relatively cool dust temperatures to dominate the radiation, thereby shifting the energy peak from near- to mid-infrared bands. This model, unlike the typical AGN hot torus models, can produce a infrared SED flattening  \nthat is consistent with LRD observations through JWST MIRI. In this scenario, there is no need for additional stellar emission or AGN scattered light.', 'ACKNOWLEDGMENTS': 'We greatly thank Changhao Chen, Seiji Fujimoto, Jinyi Shangguan, and Mingyang Zhuang for constructive discussions. We acknowledge support from the National Natural Science Foundation of China (12073003, 12003003, 11721303, 11991052, 11950410493), and the China Manned Space Project (CMS-CSST-2021-A04 and CMS-CSST-2021-A06). This work is also supported by Japan Society for the Promotion of Science (JSPS) KAKENHI (20H01939; K. Ichikawa). LCH was supported by the National Science Foundation of China (11991052, 12233001), the National Key R&D Program of China (2022YFF0503401), and the China Manned Space Project (CMS-CSST-2021-A04, CMSCSST-2021-A06). This work made use of Astropy: 2 a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration et al. 2022). This work use the dust extinction data from version 23.01 of CLOUDY , last described by Gunasekera et al. (2023).', 'REFERENCES': "Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ,  \n345, 245, doi: 10.1086/167900 Carnall, A. C., McLure, R. J., Dunlop, J. S., & Dav'e, R. 2018, MNRAS, 480, 4379, doi: 10.1093/mnras/sty2169 Casey, C. M., Akins, H. B., Kokorev, V., et al. 2024, arXiv e-prints, arXiv:2407.05094, doi: 10.48550/arXiv.2407.05094 Casey, C. M., Kartaltepe, J. S., Drakos, N. E., et al. 2023, ApJ, 954, 31, doi: 10.3847/1538-4357/acc2bc Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151, doi: 10.1111/j.1365-2966.2010.16433.x Ceverino, D., Dekel, A., Tweed, D., & Primack, J. 2015, MNRAS, 447, 3291, doi: 10.1093/mnras/stu2694 Chan, C.-H., & Krolik, J. H. 2016, ApJ, 825, 67, doi: 10.3847/0004-637X/825/1/67 Dayal, P., Volonteri, M., Greene, J. E., et al. 2024, arXiv e-prints, arXiv:2401.11242, doi: 10.48550/arXiv.2401.11242 Dey, A., Soifer, B. T., Desai, V., et al. 2008, ApJ, 677, 943, doi: 10.1086/529516 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89, doi: 10.1086/162480  \nAbdurro'uf, Larson, R. L., Coe, D., et al. 2024, arXiv e-prints, arXiv:2404.16201,  \ndoi: 10.48550/arXiv.2404.16201  \nAkins, H. B., Casey, C. M., Lambrides, E., et al. 2024, arXiv e-prints, arXiv:2406.10341,  \ndoi: 10.48550/arXiv.2406.10341  \nAlonso-Herrero, A., P'erez-Gonz'alez, P. G., Alexander,  \nD. M., et al. 2006, ApJ, 640, 167, doi: 10.1086/499800 Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, Earth, Planets and Space, 65, 213, doi: 10.5047/eps.2012.04.014  \nAstropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74 Baldwin, J. A., Ferland, G. J., Martin, P. G., et al. 1991, ApJ, 374, 580, doi: 10.1086/170146  \nBarro, G., Perez-Gonzalez, P. G., Kocevski, D. D., et al. 2023, arXiv e-prints, arXiv:2305.14418,  \ndoi: 10.48550/arXiv.2305.14418  \nBarvainis, R. 1987, ApJ, 320, 537, doi: 10.1086/165571  \nBoquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103, doi: 10.1051/0004-6361/201834156  \nCalzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692  \nXie, Y., Li, A., & Hao, L. 2017, ApJS, 228, 6, doi: 10.3847/1538-4365/228/1/6"}