The Physical Nature of the Off-centered Extended Emission Associated with the Little Red Dots

We follow the procedures in Chen et al. (2025) for PSF model construction. For each target, field stars observed with the same telescope on-sky position angle are selected by identifying sources with point-like morphology using SExtractor (Bertin & Arnouts, 1996) and comparing their SEDs with stellar spectral templates. Empirical PSF models are then generated using the pixel-based algorithm PSFEx (Bertin, 2013). The signal-to-noise ratio of the field-star images varies with wavelength, as the LW filters generally have lower signal-to-noise ratio compared to the SW filters because of the decline in stellar flux toward the redder bands. The signal-to-noise ratio of the medium-band filters is also lower that that of the broad-band filters owing to their narrower bandpass. In order to achieve PSFs with more uniform quality across different bands, we construct hybrid PSF models that combine the central core of an empirical PSF with the outer wings of a theoretical PSF generated using WebbPSF (Perrin et al., 2015). For the two reddest medium bands (F460M and F480M), even the inner core of their PSFs proved impossible to constrain, and we therefore excluded them from further consideration.

We use GalfitS (R. Li & L. C. Ho, in preparation), a parametric image modeling code that simultaneously fits multi-band images with different resolutions and pixel scales. For each source, the fit is conducted on a $1\arcsec\times 1\arcsec$ image cutout centered on the point source. The input error map contains three major components: the background, readout noise provided by the weight maps in UNCOVER DR3, and Poisson noise estimated from the mosaics. The final error map is generated adopting the procedures described in the DAWN JWST Archive²²2https://dawn-cph.github.io/dja/index.html. To remove nearby contaminants for MSA13821, we follow Chen et al. (2025) to create a $4\arcsec\times 4\arcsec$ cutout image in the F277W band that serves as the detection image, perform source detection using the Photutils.detect_source package with nstd = 0.7, npix = 15, and contrast = 0.01, and then measure and remove the local sky background with all sources masked. For MSA10686 and MSA41225, extended wings of foreground galaxies are visible within their $1\arcsec\times 1\arcsec$ cutout images. As these potentially can affect the measurement for the extended off-centered component of the LRD, especially in the LW filters, we do not mask but simultaneously fit and subtract them from the cutout images in all bands. Specifically, we model the light profiles of the nearby galaxies with a Sérsic (1968) function and the LRDs are described as a point source. Then, the nearby galaxies, LRDs, and a uniform background model are fit together in cutout images large enough to include the centers of all nearby galaxies, which in our case is $2\arcsec\times 2\arcsec$ for MSA10686 and $4\arcsec\times 4\arcsec$ for MSA41225. During the fitting, $R_{e}$, $n$, and the positions of the Sérsic and point source profiles are fixed across different bands, allowing only the brightness of the light profiles and background to vary. The extended wings from nearby galaxies as well as the local background in the final $1\arcsec\times 1\arcsec$ cutout images are subsequently subtracted using the best-fit models. Note that in the above fitting, we do not include an additional component to describe the off-centered blobs, because the morphology and brightness of these low-surface brightness components hardly can be constrained in the relatively large image cutouts. Since our goal in this initial step is to remove contamination from the extended wings of much brighter nearby galaxies, the effect of these much fainter blobs should be minimal.

For the final fitting to characterize the off-centered blobs, we adopt the cutout images with background subtracted and nearby sources masked. The cutout size is $1\arcsec\times 1\arcsec$ for all targets except for the much brighter MSA45924, for which we use a $2\arcsec\times 2\arcsec$ cutout to include the prominent PSF diffraction spikes. We model the light profile of the LRD with a point source and adopt the Sérsic profile to represent the off-centered blobs. Initial values of structural parameters are set according to the values provided in Chen et al. (2025); the fit assumes that the parameters are constant in different bands. All the fits use the Adam Optimizer as the optimization algorithm to locate the position of maximum likelihood in parameter space. We run this process 20,000 times. We then perform nested sampling around the position of maximum likelihood to infer parameter uncertainties from the posterior distribution. Figure 1 illustrates the fit for MSA10686. Results for the other three sources are shown in Appendix A.

To minimize potential bias in the SED modeling of the off-centered blobs, we adopt stricter constraints compare to Chen et al. (2025) in determining their flux upper limits. We first run a set of additional fits using only a point source to model the LRD itself. We then evaluate the goodness-of-fit in each band using the Bayesian information criterion,

For the bands that satisfy the aforementioned BIC selection criterion, we further compare the measured flux of the off-centered blobs with the detection limit of the same band. The best-fit fluxes are only taken as solid detection if they surpass the detection limit by more than the uncertainty of the best-fit flux. We report the final photometric measurements in Tables 1 and 2.

We first estimate photometric redshifts using EAzY (Brammer et al., 2008), aiming to determine whether the off-centered blobs are physically associated with the LRD point sources, or just galaxies at different redshift projected close to the line-of-sight. We use the standard set of 12 templates in EAzY generated from the stellar population synthesis code of Conroy et al. (2010). To account for the possible high-redshift nature of these blobs, we also include templates from Larson et al. (2023) that are optimized for typical galaxies in the range $4\lesssim z\lesssim 7$. Upper limits during the fits are treated as data points with both their values and errors equal to 0.5 times the actual upper limits. This configuration allows the model values in these bands to vary between zero and the upper limit with the $1\,\sigma$ confidence, with preference given to model fluxes closer to one half of the detection upper limits. The consequences of this treatment of the upper limits will be explored further in Section 4.2. No magnitude priors are adopted. We found that all four sources exhibit bimodal redshift probability distribution functions, with one peak located near the spectroscopic redshift of the LRDs and the other at a lower redshift.

To quantify the probability that the off-centered blobs lie at redshifts discordant with those of the LRDs, we follow the procedures in Bloom et al. (2002), Berger (2010), and Wu et al. (2022), assuming a uniform spatial distribution of galaxies such that the overlapping probability can be described by a Poisson distribution,

Our photometric redshift analysis suggests that, apart from MSA13821, the physical association of the off-centered blobs and their nearby LRD is ambiguous: there exists two possible solutions, one similar to the spectroscopic redshift of the LRD point source and another at lower redshift. To evaluate further the robustness of these two solutions and to explore the physical properties of these potential companion galaxies, we model the observed multi-band photometry with the SED fitting module in GalfitS. Specifically, model SEDs are calculated with the theoretical stellar templates from the Binary Population and Spectral Synthesis (BPASS V2.3; Byrne et al. 2022) code assuming a stellar initial mass function from Kroupa (2001). We adopt a continuous star formation history, which includes a stellar population undergoing continuous star formation since 200 Myr after the Big Bang, as well as a single stellar population formed in a burst that lasted for 100 Myr. In practice, the age of the single stellar population $t$ and the fraction of the total stellar mass formed during the burst ($f\_burst$) are treated as free parameters. Nebular emission from young stars is computed using the spectral synthesis code CLOUDY (Ferland et al., 1998; Chatzikos et al., 2023), assigning the same gas-phase metallicity as the stellar metallicity. The combined templates are then corrected for dust extinction using the Calzetti et al. (2000) extinction curve and absorption from the intergalactic medium using the formalism of Meiksin (2006) to produce the final galaxy SED model. Parameter ranges during the fitting are set as follows: $f\_burst$ between $10^{-2}$ to 1, metallicity $Z$ between 0.05 to 1 $Z_{\odot}$, $A_{V}$ from 0 to 5 mag, and stellar mass $M_{*}=10^{6.5}-10^{12}\,M_{\odot}$. We perform two fits, one at the spectroscopic redshift of the LRD point source and another at the best-fit photometric redshift. In each run, the age of the stellar population varies between $t=0$ and the cosmic age minus 200 Myr at the corresponding redshift.

Figure 3 presents the SED fits for both the high- and low-redshift solutions using stellar models. We slightly modify the manner in which we calculate the total $\chi^{2}$ of the fitting results. For the bands with solid detections, we sum up the residual value of $({\rm data}-{\rm model})^{2}/{\rm error}^{2}$ as in the conventional method, while for the bands with non-detections we only sum the residual values when the model fluxes exceed the upper limits. For MSA10686 and MSA13821, the low- and high-redshift solutions differ very little in terms of the total $\chi^{2}$ ($\Delta\chi^{2}<2$; $\Delta p\propto e^{\frac{\Delta\chi^{2}}{2}}<2.72$). By contrast, the low-redshift solution for MSA41225 has $\chi^{2}$ lower by 18.94, which translates to a factor of $\sim 1.3\times 10^{4}$ in terms of probability, while for MSA45924 the high-redshift solution, with $\Delta\chi^{2}=170.56$, gives a much better description to the data. Therefore, for all sources except MSA41225, considering the low overlapping probability (typically $\sim 10^{-6}-10^{-8}$), our results favor the high-redshift solutions consistent with the spectroscopic redshift of the LRD point sources over the low-redshift solutions given by EAzY, firmly suggesting that the off-axis blobs are physically associated with the LRDs. In terms of MSA41225, we are less certain about its redshift due to the relatively large probability deficit of the high-z solution implied by the $\Delta\chi^{2}$. Therefore, in this and the following sections, we proceed under the assumption that the blobs are at same redshift as their LRD point sources, and restrict our analysis to MSA10686, MSA13821, and MSA45924, whose stellar models at the LRD redshift are equally good or show significant improvement compare to their low-redshift photometric redshift solutions. The best-fit parameters assuming the LRD redshift are reported in Table 3.

The [O III] $\lambda\lambda 4959,\,5007$ and ${\rm H}\beta$ complex falls within the wavelength coverage of the medium-band filter F300M for MSA10686 and F360M for MSA13821. Both bands are detected unambiguously (Figure 3). The best-fit stellar SED model attributes the localized flux excess mostly to the [O III] doublet, with an [O III] $\lambda 5007$ equivalent width of 138 Å for MSA10686 and 256 Å for MSA13821. However, the SED models for these two targets underpredict the observed F300M flux by $22\%$ and the F360M flux by $40\%$. Assuming that all of the flux deficit originates from [O III] and that the two lines have the same width and an intensity ratio of $\lambda 5007/\lambda 4959=2.98$ (Storey & Zeippen, 2020), the true equivalent width of [O III] $\lambda 5007$ should be 276 Å for MSA10686 and 738 Å for MSA13821, akin to so-called extreme emission-line galaxies (EELGs; e.g., van der Wel et al. 2011). Figure 4 compares the stellar mass, SFR, and effective radius of the off-centered blobs, within the context of the stellar model, with a sample of high-redshift EELGs at $4<z<9$ selected using empirical color-color criteria derived from templates based on local metal-poor starburst galaxies (Llerena et al., 2024). High-redshift EELGs lie systematically above the star-forming main sequence at $z=4.5$ (Koprowski et al., 2024), and they generally follow the stellar mass-size relation at $z=5$ (Morishita et al., 2024). The off-centered blobs, if they are indeed galaxies, have properties that are consistent with those of high-redshift EELGs.

While the stellar model can generally reproduce the observed SED, it underpredicts the fluxes from the [O III]$+\rm H\beta$ complex, indicating a potentially different, nonstellar source of ionization. We explore the possibility that the off-centered emission originates from gas photoionized by the radiation from the nearby LRD point source. We calculate a grid of spectral templates using CLOUDY, adopting an open geometry with the inner and stop radius set to $r_{0,{\rm c}}-r_{{\rm e,c}}$ and $r_{0,{\rm c}}+r_{{\rm e,c}}$, respectively, with $r_{{\rm e,c}}$ the effective radius of the off-centered emission and $r_{0,{\rm c}}$ the distance from the LRD point source. The output spectrum is especially sensitive to the shape of the ionizing spectrum in the extreme-UV, which contributes the majority of the ionizing photon budget. We construct the input ionizing spectrum by interpolating the observed multi-band photometry of the LRD point source. As we cannot directly measure the continuum shorter than the Lyman limit due to intergalactic medium absorption, we assume that it follows the same power law redder than the Lyman limit³³3We have verified that under this approximation and case B recombination the number of Lyman continuum photons suffices to produce the observed ${\rm H}\alpha$ luminosities reported in Greene et al. (2024).. We calculate the spectral slopes by fitting a power-law function to all the detected broad-band photometry of the LRD point sources whose pivot wavelength is shorter than 3500 Å, as the spectral turnover for LRDs generally occurs near the Balmer limit. Only broad-band photometry is used to avoid possible contamination of emission lines to the medium bands, but we exclude F070W and F090W because they are affected by intergalactic medium absorption. Expressing the UV power law continuum as $f_{\nu}\propto\,{\nu}^{\alpha}$, $\alpha=-0.77$, $-0.91$, $-0.41$, and $-2.12$ for MSA10686, MSA13821, MSA41225, and MSA45924, respectively. The fraction of the ionizing radiation that enters the blobs is calculated using a covering factor estimated from the solid angle of a sphere centered on the LRD whose radius and position equal the effective radius and the center of the off-centered blob. We compute spectra with different combinations of the density and metallicity, with ${\rm log}\,(n/{\rm cm}^{-3})=1-4$ in steps of 1 dex and $Z=0.005$, 0.02, 0.2, 0.4, and $1\,Z_{\odot}$ assuming solar elemental abundances. We then correct the output spectra for intergalactic medium absorption following the formalism (Meiksin, 2006) used for our stellar model, before interpolating them as a two-variable function of $\log n$ and $Z$ using the Python package Scipy.interpolate.RegularGridInterpolator to obtain the final nebular spectral models. Least-squares fitting is performed to search for the nebular parameters that produce a spectrum bearing the closest resemblance to the observed SEDs. The fitting includes a free parameter $f$, which varies between 1 and 10, to scale the total flux of the model spectrum, to account for potential effects from non-spherical geometry as well as inhomogeneity in density and metallicity distribution, which are not considered in our CLOUDY calculations.

The fitting results for the nebular model are shown in Figure 5 and reported in Table 3. Overall, the nebular model describes well the SEDs of MSA10686 and MSA13821, especially in the band with a flux upturn due to strong line emission from [O III] and ${\rm H}\beta$. For MSA13821, the nebular model achieves a better fit ($\Delta\chi^{2}=9.42$) compared to the stellar model. The result from the nebular model of MSA10686 has similar goodness-of-fit as the stellar model, with slightly higher $\chi^{2}$ ($\Delta\chi^{2}=2.16$). The best-fit models for MSA10686 and MSA13821 suggest the off-center, extended blob arises from low-density ($n\approx 10\,{\rm cm}^{-3}$), metal-poor ($Z\approx 0.05-0.1\,Z_{\odot}$) gas. Note that while the best-fit density of MSA10686 reaches the lower limit of our parameter grid, we are unable to explore results of lower density, as it causes catastrophic failure in CLOUDY due to large negative optical depth in line transition.

By contrast, the nebular model failed to fit the SED of MSA45924. Most of the mismatch between data and model occurs in the bluest end of the SED sampled by F070W and F090W, where the turnover in the nebular spectrum due to two-photon decay underpredicts the observed flux.

Off-centered, extended emission is associated with a significant fraction of LRDs. In their systematic study of 99 photometrically selected LRDs at $z\approx 4-8$, Rinaldi et al. (2025) find that $\sim 30\%$ show a variety of complex, asymmetric UV morphologies beyond the dominant central point source. This is consistent with the results of Chen et al. (2025), based on much meager statistics. What is the physical nature of this asymmetric, resolved component? Assuming that the off-centered emission near LRDs is of stellar origin, we show that their inferred SFRs, stellar masses, and effective radii qualitatively resemble those of high-redshift EELGs. Under this interpretation, the off-centered blobs may be regarded as companion galaxies, which are presumably interacting or merging with a large fraction of the LRDs. While LRDs lack a galaxy component (Chen et al., 2025) massive enough to be commensurate with the local scaling relations between black hole and galaxy mass (Kormendy & Ho, 2013; Greene et al., 2020), gas tidally stripped from these star-forming galaxies could be the fueling channel that maintains these black holes in an active state. Alternatively, if the off-centered blobs comprise stars of the LRD host galaxies themselves, then their global asymmetry can arise from dynamical processes that can dislodge the black hole to kpc-scale distances, such as gravitational recoil (Bekenstein, 1973; Campanelli et al., 2007) or gravitational slingshot effects (Hoffman & Loeb, 2007), indicating that the central black holes in the LRDs are produced by merging from multiple lower mass black holes. The stellar scenario offers the best interpretation for the off-centered emission in MSA45924.

Extended structures observed in broad-band images, especially near luminous AGNs like the LRDs, are not necessarily of stellar origin. We explore an alternative possibility, in which the off-centered blobs constitute nebular gas directly ionized by the radiation field emanating from the adjacent LRD. Our photoionization calculations show that this scenario is viable, for two of the four sources studied here. Using the observed SED of the LRD as the ionizing radiation source and the relative displacement between the off-centered component and the point source, the nebular model can very well explain the observed flux and SED shape of MSA10686 and MSA13821, even better than the stellar model in the case for MSA13821. The nebular model has notable success in accounting for the high equivalent width of the [O III] line detected in the medium-band filter F300M for MSA10686 and F360M for MSA13821. This is reminiscent of the situation in high-redshift galaxies whose hard ionizing radiation field produces sufficient nebular emission to affect the broad-band flux (e.g., Cameron et al., 2024; D’Eugenio et al., 2025; Kataria et al., 2025; Solimano et al., 2025). For these two targets, the best-fit nebular model gives a density of $\sim 10\,{\rm cm}^{-3}$, with radius-averaged electron temperature of $\sim 10^{4}\,{\rm K}$. Since the best-fit density already reaches the lowest value of our photoionization calculations, the real density could be lower. The low density and low surface brightness of the nebular emission qualitatively resembles the warm ionized medium ubiquitously detected in local spiral galaxies (e.g., Wang et al. 1998). Interestingly, although under the nebular scenario the extended emission does not require necessarily the presence of any galaxy, the best-fit metallicity for MSA10686 and MSA13821 are both $Z\approx 0.05-0.1\,Z_{\odot}$, much higher than that of the intergalactic medium at cosmic noon. According to Simcoe et al. (2004) and McQuinn (2016), only $\lesssim 10\%$ of sight lines with O IV or C IV absorption at $z\approx 2.5$ have metallicities higher than $0.01\,Z_{\odot}$. At redshifts closer to the sources in our sample, the 95% upper limit of the intergalactic medium [Mg/H] ratio, defined as the logarithmic number ratio of magnesium atoms to hydrogen atoms relative to the solar ratio, is $-3.73$ at $z=6.47$ (Tie et al., 2024), suggesting $Z<10^{-3}Z_{\odot}$. The nebular gas could have been chemically enriched by stars in a preexisting, albeit undermassive (Chen et al., 2025), host galaxy. or by star formation connected with the accretion disk (e.g., Artymowics et al., 1993; Zhang et al., 2025).

Future JWST observations with spatially resolved spectra, such as those that can be obtained with the integral-field units of NIRSpec or via slitless spectroscopy with NIRCam or NIRISS, can provide more definitive diagnostics of the physical nature of the extended, off-centered emission associated with LRDs. For example, our nebular model for MSA1068 predicts prominent [S II] $\lambda\lambda 6716,\,6731$ emission, similar to the situation seen in the warm ionized medium in the local Universe that shares similar density and temperature as our best-fit model spectrum. Detection of [S II] and measurement of the [S II]/${\rm H}\alpha$ ratio for the extended emission would provide smoking gun support for their nebular nature.

We estimate the total mass of ionized gas in both models. Under the stellar scenario, in which the continuum is assumed to be produced by stars, photometry in F300M and F360M provides constraints on the total flux of the [O III]$+{\rm H}\beta$ complex in MSA10686 and MSA13821, respectively. We convert this value to ${\rm H}\beta$ luminosity using the relation between [O III]/${\rm H}\beta$ ratio and specific star formation rate derived in Dickey et al. (2016). The gas mass is then calculated using the dereddened ${\rm H}\beta$ luminosity, assuming $n_{e}=100\,{\rm cm}^{-3}$, $T=10^{4}\,{\rm K}$, and an effective recombination coefficient of $\alpha_{{\rm eff},{\rm H}\beta}=3.03\times 10^{-14}\,{\rm cm}^{-3}$ (Draine, 2011),

If the off-centered blobs near MSA10686 and MSA13821 are galaxies, their estimated ionized gas mass is ${\rm log}\,(M_{\rm gas}/M_{\odot})=6.96\pm 0.05$ and $6.93\pm 0.12$, respectively⁴⁴4We bootstrap all parameter uncertainties through each step of the calculation and take the final standard deviation as the error.. Under the nebular scenario, the broad-band emission is attributed to the recombination continuum. In this case, we estimate the total ionized gas mass using the best-fit electron density and total emission volume of our CLOUDY model, multiplied by the best-fit scaling factor $f$. For MSA10686 and MSA13821, their best-fit nebular templates correspond to ${\rm log}\,(M_{\rm gas}/M_{\odot})=8.71\pm 0.28$ and $7.90\pm 0.54$, higher than the prediction from the stellar model by more than an order of magnitude.

If the off-centered gas clouds on kpc scales are infalling within the halo of the LRDs, their gas mass provides a rough estimate of the LRD halo mass accretion rate. The infall time of the gas clouds onto the LRD halo can be estimated as

These estimates depends on our highly incomplete knowledge of the halo masses of LRDs. Clustering analysis of 27 low-luminosity broad-line AGNs discovered by JWST at $5<z<6$, among them LRDs, constrain their typical halo mass to $M_{h}\approx 3\times 10^{11}\,M_{\odot}$ (Arita et al., 2025). On the other hand, Jiang et al. (2025) recently proposed a model for the formation of seed black holes in LRDs via the gravothermal collapse of the cores of self-interacting dark matter (SIDM) halos. Within this framework, core collapse occurs most efficiently within a relatively narrow range of halos with masses $M_{h}\approx 10^{8}-10^{10}\,M_{\odot}$, significantly lower than that suggested by the clustering analysis of Arita et al. (2025). In Figure 6, we present the value of $t_{\rm dyn}$ and $\dot{M}_{\rm acc}$ as functions of $M_{\rm vir}$ for halos with $c=3$, which is typical for the conventional structure formation framework assuming cold dark matter (CDM), and halos with $c=15$ as proposed by Jiang et al. (2025) in the SIDM scenario.

For typical CDM halos with $c=3$ at $z=6$, if the $3\times 10^{11}\,M_{\odot}$ halo mass estimated by Arita et al. (2025) can be applied to typical LRDs, the infall time of a $10^{8}\,M_{\odot}$ gas cloud starting at 1 kpc scale is $t_{\rm dyn}\approx 11\,{\rm Myr}$, and the accretion rate $\dot{M}_{\rm acc}\approx 9\,M_{\odot}\,{\rm yr}^{-1}$. Alternatively, adopting the lower halo mass ($10^{8}-10^{10}\,M_{\odot}$) predicted for the SIDM halos with $c=15$, the infall time increases to $t_{\rm dyn}\approx 12-64\,{\rm Myr}$, and the accretion rate drops to $\dot{M}_{\rm acc}\approx 2-8\,M_{\odot}\,{\rm yr}^{-1}$. For an LRD at $z=6$, if the halo accretion rate is sustained to $z=4$, the accreted gas would form a galaxy with a total stellar mass $M_{*}\approx 10^{9}\,M_{\odot}$ assuming a star formation efficiency of 0.1. The stellar mass would be even higher if stars form more efficiently at high redshift (e.g., Inayoshi et al. 2022; Ceverino et al. 2024). The large stellar mass assembled would gradually render the LRD spatially resolved, which explains the rapid decline in the number density of compact sources with the characteristic “V-shape” SEDs at $z<4$ reported by both observations (Kocevski et al., 2025; Ma et al., 2025) and theoretical predictions (Inayoshi, 2025).