Towards Self-consistent Modeling of Quasi-Stars and their NLTE Spectra: Implications for Little Red Dots

abstract · abstract

Quasi-stars, massive envelopes surrounding rapidly accreting black hole (BH) seeds, are promising candidates for explaining the rapid formation of supermassive black holes in the early Universe and may account for some of the JWST-discovered Little Red Dots (LRDs). We present a self-consistent framework combining stellar evolution calculations with non-local thermodynamic equilibrium (NLTE) spectral synthesis to test this connection. Our models drop the polytropic assumptions of previous work: the envelope is described entirely by mixing-length theory, with saturated convection and radiative energy setting the luminosity at an inner boundary a few hundred Schwarzschild radii from the BH. This removes the BH-mass limits of earlier models, and we follow the evolution up to BH-to-total mass ratios approaching unity. We then use TLUSTY to compute the first converged NLTE atmospheric models for the surface conditions and primordial hydrogen-helium composition expected in quasi-stars. Our models produce cool photospheres, prominent Balmer breaks, and Balmer and Paschen absorption features. NLTE effects substantially modify line strengths, increasing Balmer absorption depths by $\sim 35-55\%$ relative to LTE models, while having a minor effect on the continuum. Comparison with observed LRD spectral energy distributions shows broad agreement in continuum shape, although our models do not reproduce their commonly observed strong emission lines. Under the quasi-star interpretation, comparisons with individual LRDs suggest characteristic quasi-star masses of approximately $10^5-10^7 M_\odot$. These results establish a physically motivated connection between quasi-star evolution and observable spectra, providing a framework for testing quasi-stars as progenitors of LRDs while identifying limitations of current atmospheric models.

S5: V Summary / block 1 · html

The incorporation of NLTE effects into our quasi-star spectra strengthen Balmer line absorption by a factor on the order of ∼ 35 − 55 % \sim 35-55\% .

S4.SS1: IV.1 Quasi-star Evolutionary Model / block 1 · html

Figure 1 shows the evolution of various physical properties for a 10 5 ​ M ⊙ 10^{5}M_{\odot} quasi-star, comparing the results of our new model with the polytrope-based ones from [ 18 ] and the purely analytic model from CB24. The surface properties of our current model largely agree with those of [ 18 ] , although deviations begin to manifest when the latter nears its limiting black hole mass. Interestingly, despite our new model making no polytropic assumptions in its derivation, its mid-envelope values of p gas / p rad p_{\rm gas}/p_{\rm rad} and γ ad \gamma_{\rm ad} are closer to the predictions of CB24’s polytropic model than [ 18 ] ; namely, γ ad \gamma_{\rm ad} tends toward 4/3 over time in our new model, whereas the old model asymptotically diverged from it.

S4.SS1: IV.1 Quasi-star Evolutionary Model / block 1 · html

Figure 3 : Several structural profiles for our fiducial quasi-star. We show the results at M BH = 0.5 ​ M ⋆ M_{\rm BH}=0.5M_{\star} , as well as the comparison between our new and old models at M BH ≈ 0.17 ​ M ⋆ M_{\rm BH}\approx 0.17M_{\star} .

S4.SS2: IV.2 Quasi-star Spectra / block 1 · html

α MLT \alpha_{\rm MLT} is set in TLUSTY via the HMIX0 parameter. For the sake of rigor, we aimed to match it with the value used for the corresponding MESA model. However, we found that models with higher values of HMIX0 were generally easier to converge than those with lower values. As a result, a few of our TLUSTY models do not match the mixing length parameters of their MESA models (see Table 1 ). Fortunately, the value of HMIX0 has little influence on our spectra, in contrast to its significance for the MESA model. This is because convection becomes inefficient near the surface of a quasi-star, especially in the spectrum-forming layers; indeed, what matters for the spectrum is the ( T eff , log ⁡ g ) (T_{\rm eff},\log g) pair, and this is provided by MESA models that have the intended α MLT \alpha_{\rm MLT} .

S4.SS3: IV.3 Comparison with LRDs / block 1 · html

Figure 11 : Comparison of several LRDs with the midlife spectrum of a “varying” 10 6 ​ M ⊙ 10^{6}M_{\odot} model with α MLT = 1 \alpha_{\rm MLT}=1 , both with and without attenuation ( A V = 2 A_{V}=2 ). Flux is normalized at 0.52 ​ μ ​ m 0.52\mu\mathrm{m} . Strong attenuation generally results in a poorer match between our quasi-star spectrum continua and the observations from LRDs.

S4.SS3: IV.3 Comparison with LRDs / block 1 · html

Figure 12 : Comparison of LRDs with midlife quasi-star models, without flux normalization. Results are shown for both α MLT = 1 \alpha_{\rm MLT}=1 and 2, as well as with “varying” (solid) and “fixed” (dashed) opacity treatments.

S5: V Summary / block 1 · html

We have presented significant developments in the evolutionary modeling of quasi-stars with MESA , the first non-LTE atmospheric models based on numerical simulations of quasi-stars from first principles, and the first successfully converged non-LTE TLUSTY models at the conditions and compositions expected of quasi-stars. We have provided modeled spectra for a multitude of quasi-star evolutionary constructions, including different compositions, mixing length parameters, opacity treatments, and a model with the surface density inversion numerically suppressed.

S5: V Summary / block 1 · html

Our main results can be summarized as follows:

S5: V Summary / block 1 · html

The upper limit on the BH mass for older models (CB24, Hassan et al. 18 ) appears to be caused by the enforcement of saturated convection out to high radii. Our current model, which utilizes saturated convection for the luminosity near the BH but models the rest of the quasi-star with MESA ’s mixing length theory, can describe the full lifetime of the quasi-star up to M BH / M ⋆ ∼ 1 M_{\rm BH}/M_{\star}\sim 1 .

S4: IV Conclusions & Discussion / block 1 · html

We have presented new MESA quasi-star evolutionary models generated with the MESA-QUEST framework and compared their synthetic spectra to JWST observations of LRDs. Our results demonstrate that late-stage quasi-stars, with black holes comprising ≳ \gtrsim 10% of their total mass and photospheric temperatures near 6000 K, naturally produce the defining continuum features of LRDs: large Balmer breaks, flat or red near-infrared slopes, and weak ultraviolet emission. These properties emerge self-consistently from the quasi-star’s radiation-supported envelope, and do not require an old stellar population or extreme dust reddening.

S2.SS1: II.1 Quasi-star Modeling / block 1 · html

Previous works modeling the accretion and resulting luminosity report only a factor of < 2 <2 deviation from the object’s Eddington limit throughout the quasi-star’s lifetime [ 6 , 13 ] . The dimensionless scaling factor α \alpha serves as a proxy for a range of physical effects that enhance or suppress accretion such as rotation, magnetic fields, photon-trapping, etc., which we plan to implement and explore more thoroughly in future work. Here we adopt a simple scaling factor α \alpha and use α = 1 \alpha=1 as a first approximation.

S4: 4 Conclusions / block 1 · html

As a pathway to understand the physical origin of the rest-optical emission from LRDs, we have calculated stellar evolution models of rapidly accreting stars. We start from low mass-proto-stars that experience a net mass gain and follow their evolution until they become SMS and reach the GRI, which leads to the formation of so-called QS ( Begelman et al., 2006 ; Begelman et al., 2008 , where a stellar envelope is supported by a central BH, see) . We considered maximum mass gain rates of 0.001, 0.1 and 1 M ⊙ {\rm M}_{\odot} /yr, and metallicities from Z = 0 Z=0 to 0.01. We show that our scenario works independently of metallicity.

S4: 4 Conclusions / block 1 · html

We are aware that the employed mass gain rates are in agreement with the mass-inflow rates onto simulated haloes from the cosmological simulations. For example, Cenci and Habouzit (2025) reports inflow rates range from ∼ 10 − 3 \sim 10^{-3} to ≳ 1 \gtrsim 1 M ⊙ {\rm M}_{\odot} /yr at z ≃ 4 z\simeq 4 , which increase with higher redshifts. The mass growth threshold can also be reached in dense and compact forming massive star clusters with stellar surface densities similar to those estimated for LRDs (typically > 10 5 >10^{5} M ⊙ {\rm M}_{\odot} pc ) − 2 {}^{-2}) , where both gas-accretion and collisions can concur to the formation of SMS, with no limitation in terms of metallicity, contrary to the DCBH scenarios ( Gieles et al., 2018 ; Lahén et al., 2025 ; Baggen et al., 2026 ; Ramírez-Galeano et al., 2025 , and references therein) .

S3.SS3: III.3 Driving Mechanism / block 1 · html

These results indicate that quasi-stars are driven primarily by the κ \kappa -mechanism in the He II ionization zone. This occurs despite their extremely cool effective temperatures and metal-poor composition, which preclude driving by the deep iron opacity bump ( T ∼ 2 × 10 5 T\sim 2\times 10^{5} K) found in massive stars. Consequently, their pulsational physics aligns more closely with supermassive analogs of classical Cepheids rather than purely H-driven Mira variables [ 17 , 37 ] .

S4: IV Discussion / block 1 · html

The recent discovery of the lensed Little Red Dot RX1 [ 38 ] provides the first compelling observational counterpart to a pulsating quasi-star phase. Our linear stability analysis confirms that such envelopes are generically unstable to radial pulsations, yielding fundamental modes with periods of ∼ 30 \sim 30 – 180 180 years and first-overtone modes with periods of ∼ 10 \sim 10 – 30 30 years. The latter is remarkably consistent with the ∼ 30 \sim 30 -year variability cycle (in the co-moving frame) inferred from the ∼ 130 \sim 130 -year baseline provided by the lensed images. Crucially, Zhang et al. [38] report that RX1 traces a counter-clockwise hysteresis loop in the luminosity–temperature plane, a behavior that is difficult to reconcile with stochastic accretion-driven variability, but is a well-known hallmark of thermodynamic work cycles associated with stellar pulsations.

S5: V Discussion / block 1 · html

The picture of LRDs as quasi-stars carries with it an evolutionary sequence that includes stages prior to and following the phenomena we have focused on in this paper. We have illustrated this sequence schematically in Fig. 1 . The SMS stage should last only ∼ 2 − 5 % \sim 2-5\% of a late-stage quasi-star lifetime and should thus be rarest on the sky. They should also be ∼ 10 − 100 \sim 10-100 times less luminous on average, i.e. L ∼ 10 42 − 10 43 L\sim 10^{42}-10^{43} erg s -1 , with color temperatures ranging from ∼ 5000 \sim 5000 K at lower masses to ∼ 10 4 \sim 10^{4} K as masses approach 10 5 ​ M ⊙ 10^{5}M_{\odot} [ 17 ] . Early-stage quasi-stars should be more common than SMSs ( ∼ 10 % \sim 10\% compared to late-stage quasi-stars) and — at constant envelope mass — should evolve from hotter (several × 10 4 \times 10^{4} K) to cooler color temperatures as the BH grows inside [ 11 ] . However, this evolution could be obscured by the diversity of quasi-star masses and mass ratios. Efforts should be made to find the immediate precursors to LRDs. Little is known about the fate of a quasi-star once the BH swallows most or all of the envelope, e.g., whether the residual gas disperses or is able to settle into a steady state with a much lower envelope mass if matter continues to be supplied from the environment at a sufficient rate [ 32 , 12 , 21 , e.g., resembling one of the stand-alone models as in] . Alternatively, if stripped of their envelopes these immediate post-quasi-stars could be the first proper AGNs in the universe — examples of which may have already been discovered [ 9 , 34 , e.g.] .

S2: II Late-stage quasi-stars / block 1 · html

implying Δ r / R ∗ ≈ 0.2 A 23 − 2 / 15 β − 1 14 / 75 M 6 17 / 150 \Delta r/R_{*}\approx 0.2A_{23}^{-2/15}\beta_{-1}^{14/75}M_{6}^{17/150} . A numerical solution assuming fiducial parameters M 6 = β − 1 = 1 M_{6}=\beta_{-1}=1 and calculating n e / n n_{e}/n self-consistently gives nearly identical results, which is not surprising since the ionization fraction is 0.96.

S2.SS3: II.3 Envelope binding energy and collapse energetics / block 1 · html

As this is a key intermediary result of our manuscript, with important implications on the plausibility that collapsed SMDSs lead to late stage quasi-star like configurations, we validate it analytically in Appendix A , where we derive useful closed form approximations for the binding energy of a radiation dominated shell using a modified form of the virial theorem for a truncated envelope.

S3.SS2: III.2 Envelope inflation energetics / block 1 · html

Thus the desired inflation requires an energy input comparable to, but slightly below, the original envelope binding energy. Motivated by this result, in what follows we will assume Δ ​ E infl ∼ | E bind | \Delta E_{\rm infl}\sim|E_{\rm bind}| .

S5: 5 Summary / block 1 · html

In the quasi-star model presented in this paper, the rest-UV radiation observed in LRDs is mostly due to the contribution of a low-mass ( M ∗ ∼ 10 8 ​ M ⊙ M_{\ast}\sim 10^{8}\ M_{\odot} ) star-forming galaxy, in the assumption of a SSP with stellar age lower than 1 Gyr and low metallicity ( Z = 0.1 ​ Z ⊙ Z=0.1\ Z_{\odot} ). Different hypotheses would clearly bring to different estimates of the stellar mass. More in detail, the assumption of a lower metallicity ( Z = 0.05 ​ Z ⊙ Z=0.05\ Z_{\odot} ) would produce an increase of nearly 0.5 dex in the estimated stellar mass.

S5: 5 Summary / block 1 · html

Within these assumptions on the mass ratio between the SMBH and the quasi-star, and on the properties of the host galaxy, we obtain that the SMBH at the centre of the quasi-stars represent between 1% and 0.1% of the total mass of their host galaxy. This value is lower than those obtained through the line broadening and by dynamical considerations, but is still in slight tension with the local relation by Reines and Volonteri (2015) .

S2.SS1: II.1 Simulation and Initial Conditions / block 1 · html

At z = 1100 z=1100 , when baryons and photons start to decouple, we include an additional uniform baryon matter field with the same resolution as PDM in the initial condition set up, resulting in a total number of N = 2 × 256 3 N=2\times 256^{3} particles with a mass of ∼ 100 ​ ( 18 ) ​ M ⊙ \sim 100(18)~M_{\odot} for PDM (gas) particles. With this setting, we both approximate the initial structure of PDM around the PBH, and let the gas collapse onto this PBH-seeded halo. In running these simulations, the softening length of PDM and gas is set to ϵ PDM = ϵ gas ∼ 0.01 ​ L / N PDM 1 / 3 ≃ 0.01 ​ h − 1 ​ kpc \epsilon_{\rm PDM}=\epsilon_{\rm gas}\sim 0.01\,L/N_{\rm PDM}^{1/3}\simeq 0.01~h^{-1}\rm kpc , and the initial chemical abundances are assigned to the values predicted for the intergalactic medium (IGM) at z = 1100 z=1100 , as summarized in [ 36 ] .

S2.SS1: II.1 Simulation and Initial Conditions / block 1 · html

In addition, we also consider the relative streaming motion between gas and PDM particles by assigning a universal velocity offset. In our earlier work [ 135 ] , we found that streaming can enhance, rather than delay, galaxy formation—contrary to the conclusion generally reached within the standard Λ \Lambda CDM framework [ 106 , 105 , e.g.,] . This enhancement arises from a displacement in the center of mass (away from the PBH) caused by the relative motion between baryons and dark matter, leading to star formation in the dense gas structures that form in the wake of the PBH. These wake-driven overdensities promote gravitational collapse and thus play a catalytic role in early star formation. To systematically explore the impact of streaming velocity on black hole formation, we implement a velocity offset parameterized as a multiple of the root-mean-square streaming velocity σ b ​ χ = 30 ​ km \sigma_{\rm b\chi}=30~\rm km\, s -1 . Specifically, we consider four cases with velocity offsets of 0.4, 0.8, 1.2, and 1.6 σ b ​ χ \sigma_{\rm b\chi} , corresponding to simulation runs labeled as PBH_LW_wstr_fd005 , PBH_LW_str_fd005 , PBH_LW_mstr_fd005 and PBH_LW_sstr_fd005 , respectively.

S4: IV Lyman-Werner Radiation from the Companions / block 1 · html

The above estimates for the LW radiation field do not account for self-shielding effects [ 134 ] due to a column of gas between the blue and red components. We test whether this affects our conclusions by estimating the timescale for molecular hydrogen dissociation in this column. For a red component exposed to an intrinsic LW radiation field J 21 , LW J_{\rm 21,LW} produced by its companion(s), the corresponding LW photon production rate of the companions can be written as N ˙ LW = 3.5 × 10 48 ​ ( 4 ​ π ) ​ ( d eff / kpc ) 2 ​ J 21 , LW \dot{N}_{\rm LW}{=}3.5\times 10^{48}(4\pi)(d_{\rm eff}/\rm kpc)^{2}J_{\rm 21,LW} photons s -1 . Assuming the LW photons dissociate H 2 \rm H_{2} with 10% efficiency, we can estimate the propagation of the H 2 dissociation front through a sphere with radius d eff d_{\rm eff} , following equation (2) of Sullivan et al. [117] , as

S6: VI Conclusions / block 1 · html

In this work, we have identified candidate UV-bright companions to a substantial fraction of LRDs (43% increasing to 72% when including tentative companions), assumed a physical connection between the components, and used their luminosities and projected separations to infer the local Lyman-Werner radiation fields experienced by the compact red sources. The inferred radiation fields are sufficiently intense to plausibly suppress molecular cooling and enable rapid collapse to very high densities, in line with established theoretical models. While this agreement does not constitute a proof of the scenario, it demonstrates that the picture is internally consistent and observationally viable.