Recent James Webb Space Telescope observations hint at unexpectedly intense cosmic star formation in the early Universe, often attributed to enhanced star formation efficiencies (SFEs). Here, we analyse the SFE in THESAN-ZOOM, a novel zoom-in radiation-hydrodynamic simulation campaign of high-redshift (z greater than or similar to 3) galaxies employing a state-of-the-art galaxy formation model resolving the multiphase interstellar medium (ISM). The halo-scale SFE (is an element of(halo)*) the fraction of baryons accreted by a halo that are converted to stars follows a double power-law dependence on halo mass, with a mild redshift evolution above M-halo greater than or similar to 10(95) M-circle star. The power-law slope transitions from similar to 2/3 to similar to 1/3 as halo mass increases, which hints at a transition from energy-driven to momentum-driven outflow. is an element of(halo)* is a factor of 2-3 larger than commonly assumed in empirical galaxy formation models at M-halo less than or similar to 10(95) M-circle star. On galactic (pkpc) scales, the Kennicutt-Schmidt relation of neutral gas is universal in THESAN-ZOOM, following Sigma(SFR) alpha Sigma(2)(gas), indicative of a turbulent energy balance in the ISM maintained by stellar feedback. The rise of is an element of(halo)* with halo mass can be traced primarily to increasing gas surface densities in massive galaxies. These results are is an element of(halo)* primarily increasing gas galaxies. robust against variations in numerical resolution and star formation and feedback models, depending mainly on the total feedback momentum budget. Although the increase in is an element of(halo)* with redshift is modest, it is sufficient to explain the large observed number density of UV-bright galaxies at z >= 12. However, reproducing the brightest sources at M-UV less than or similar to-21 may require extrapolating the SFE beyond the halo mass range covered by THESAN-ZOOM.

Line intensity mapping (LIM) has garnered attention as a powerful cosmological probe, with current-generation instruments such as SPHEREx (Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Explorer) capable of mapping the evolution of large-scale structure during the epoch of reionization (EoR). Lyman-alpha (Ly alpha) emission in the EoR is strongly shaped by resonant absorption from neutral hydrogen in the diffuse intergalactic medium (IGM), which transforms galactic sources into a low surface-brightness background. In this work, we leverage the state-of-the-art THESAN cosmological simulations to produce high-resolution theoretical predictions for future Ly alpha LIM studies, constructing continuous light-cones for line-of-sight cosmological integrations. We assess the contributions of recombination, collisional excitation, and unresolved HII regions to the total Ly alpha spectral intensity. In addition, we explore the IGM in absorption at different redshifts using damping wing analysis. We produce channel maps exploring spatial fluctuations across redshift bands probe-able by LIM instruments. We find that the slope of the absorption-included Ly alpha fluctuation power spectrum at smaller scales (k greater than or similar to 2 & times; 10(-2) arcsec(-1)) steepens toward lower redshift, and that our emission-only Ly alpha power spectrum lies above the SPHEREx sensitivity, whereas the absorption-included signal is similar to 7-8 orders of magnitude lower-providing a conservative lower limit on inhomogeneity signatures and highlighting the importance of including resonant scattering in the future. We also find that including outflows in a simple toy model boosts power by four orders of magnitude for a channel map spanning z is an element of [5.5, 6.5]. We identify limitations in our analysis and propose next steps, including incorporating the effects of resonant Ly alpha scattering and line interlopers, as well as larger simulation volumes.

Understanding the sources and evolution of cosmic reionization remains a central challenge in astrophysics, with the escape of ionizing Lyman-continuum (LyC) photons from early galaxies representing a major uncertainty. In this work, we use more than 35,000 galaxy realisations from the THESAN-ZOOM cosmological radiation-hydrodynamic simulations to identify indirect diagnostics of the LyC photon escape fraction (f_esc) and the LyC photon escape rate (Ṅ_ion,esc) across the redshift range z=3-16. We train random forest regression models using these diagnostics to predict both quantities. We present four models: two trained with the full set of simulation-derived indicators to predict f_esc and Ṅ_ion,esc, and two restricted to observables accessible to JWST photometric surveys. We find the 10-to-100Myr star-formation rate ratio (SFR_10 / SFR_100) and the gas-to-stellar mass ratio (M_gas / M_*) to be the strongest diagnostics of f_esc, suggesting a strong relationship between ionizing photon escape and gas clearing through bursty star formation. In contrast, rest-frame UV (1500 Å) absolute magnitude (M_UV) dominates Ṅ_ion,esc prediction. Motivated by the strong predictive power of M_UV, we combine observed UV luminosity functions with derived Ṅ_ion,esc - M_UV relations to construct histories of reionization. These are consistent with observational constraints, avoiding the recently reported crisis in the ionizing photon budget. Our analysis suggests that the bulk of reionization occurred rapidly after z ≈ 8, driven by UV-bright galaxies, with the M_UV < -17 populations providing the dominant contribution.

Metal absorbers represent a powerful probe of galaxy feedback and reionization, as highlighted by both observational and theoretical results showing an increased abundance of low-ionized metal species at higher redshifts. The origin of such absorbers is currently largely unknown because of the low number of galaxy counterparts detected, suggesting that they might be surrounded by low-mass faint sources below the current detection threshold. We use the THESAN-ZOOM radiation hydrodynamic simulations to investigate the connection between properties of neutral oxygen (O-1) absorbers and galaxies across redshift z approximate to 5-8 . We find that the virialized gas in haloes becomes progressively ionized with cosmic time, leading to a decrease of ti 0.2 in the covering fraction of neutral oxygen, while the total oxygen covering fraction remains constant. Comparing the O-1 line density obtained from our covering fractions with the trend suggested by blind quasar observations, we determine that the observable absorbers (N-OI greater than or similar to 10(13) cm(-2)) are not confined to haloes: atz >= 5 the majority (greater than or similar to 60 per cent ) arise beyond Rvir, consistent with recent James Webb Space Telescope (JWST) results. Close to O-1 absorbers, low-mass galaxies (M-star <= 10(8) M-circle dot) are more commonly found, while, depending on the simulated environment, we do not exclude the possibility of nearby more massive star-forming sources (>= 5 M-circle dot yr(-1)) similar to those suggested by the latest ALMA observations. These results establish O-1 absorbers as sensitive tracers of the evolving ionisation structure around faint galaxies to be probed by forthcoming deep spectroscopic surveys.

It has been suggested that "Little Red Dots" (LRDs) might be accreting black holes enshrouded by dense gas in a nearly closed geometry, which completely covers the central black hole, leading to an atmosphere-like structure known as the "black-hole star" (BH^⋆). We test this scenario by analysing new JWST spectroscopy in the far ultraviolet (FUV, rest-frame) of the prototypical LRD Abell2744-QSO1, at z=7.04. We found the presence of broad Lyα emission with an FWHM of ∼ 1000 km/s, and detections of OI, CIV, and/or FeII emission lines. The NIRCam imaging and NIRSpec slit images indicate that the low-velocity component (v≲ 200 km/s) of Lyα is likely spatially extended, but the high-velocity component (v≳ 200 km/s) of Lyα remains unresolved. Based on the multi-component kinematics and flux of Lyα relative to Balmer lines, we conclude that the observed line profile is unlikely to be broadened by subsequent resonant scattering through the interstellar medium. This suggests that the high-velocity component of Lyα originates in the broad-line region, although resonant scattering in the dense gas likely makes Lyα broader than Hα as observed. The nebular features of this LRD indicate that there is at least one relatively optically thin direction where Lyα can escape from the broad-line region (BLR). We also found indications that photons from the BLR are powering fluorescence of FeII and OI on a larger physical scale. The FUV features thus challenge the fully-covered geometry interpretation and suggest that there are "holes" in the BH^⋆, or the absorbing medium is simply clumpy.

Isolated star-forming galaxies require inflows of fresh gas from the surrounding medium to sustain episodes of star formation over time. However, there are very few direct detections of accretion onto external galaxies. Studies in absorption can only observe along limited sightlines, while those in emission can have difficulty distinguishing inflowing gas in the foreground of the galactic disk from similarly Doppler-shifted outflowing gas in the background. We explore the possibility of using the Balmer decrement (H alpha/H beta) in low-inclination systems as a diagnostic for disentangling the flow geometry in disk-like galaxies. We leverage mock spatial-spectral observations of an isolated Milky Way-mass galaxy simulated using the radiation-hydrodynamics code AREPO-RT and post-processed with the Monte Carlo radiative transfer code COLT. We find that gas components located in front of the disk exhibit systematically lower Balmer decrements than gas embedded in or behind the disk, with a mean front-back offset of Delta(H alpha/H beta) approximate to -0.14. The ability to differentiate between the disk and far-side components is limited by the extremely clumpy, multiphase dust distribution along the line of sight introducing substantial scatter. Overall, the results provide a useful observational diagnostic of inflow and outflow in dusty face-on galaxies.

Monte Carlo radiative transfer (MCRT) is widely used to model Lyman-alpha (Lya) resonant-line transport, but convergence is difficult to assess in optically thick media where photons undergo many scatterings before escape. This is especially important for internal quantities such as radiative acceleration and the force multiplier, which depend on momentum deposition throughout the gas rather than only on emergent spectra. We study the convergence of Lya MCRT momentum-transfer estimators in static spherical clouds. We first establish diffusion-limit benchmarks for radial acceleration profiles and integrated force multipliers, then develop a moment-based framework for diagnosing convergence from the photon-packet contribution distribution. This framework separates three distinct questions: whether the estimator converges to the correct mean, how large its finite-sampling uncertainty is, and whether the estimated uncertainty is itself stable. We apply this hierarchy to the direct event-based scattering estimator, a gradient-of-energy-density estimator, and a divergence-of-radiation-pressure estimator. Zeroth-order convergence is assessed with profile comparisons, integrated force-multiplier bias, and finite-group relative error. First-order convergence is quantified with fractional error, the photon number required to reach a target precision, and the corresponding runtime requirement. Second-order convergence is tested with the coefficient of variation of variance, which measures the reliability of the variance estimate used in the first-order diagnostics. Core-skipping prescriptions, source geometry, estimator construction, and spatial resolution enter this hierarchy in different ways. Our results provide a practical convergence framework for internal Lya MCRT force calculations and show why statistical precision, computational cost, and physical accuracy must be evaluated separately.

Ultra-faint dwarf galaxies (UFDs) are among the smallest and oldest galaxies in the Universe and are widely regarded as relics of cosmic reionization. To investigate how reionization quenches star formation and shapes the diversity of UFDs, we present a suite of eight cosmological zoom-in simulations of isolated UFDs with present-day halo masses of ∼10^9 M_⊙. The simulations are performed with the radiation-magnetohydrodynamic galaxy formation framework Realistic ISM modeling in Galaxy Evolution and Lifecycles (RIGEL), coupled to realistic large-scale radiation fields extracted from the THESAN reionization simulation. Despite residing in similar z=0 halos, the simulated galaxies span nearly two orders of magnitude in stellar mass and broadly reproduce the observed luminosities, sizes, metallicities, and stellar kinematics of Local Group UFDs. We find that reionization quenches star formation through a two-stage process. The arrival of the ionization front rapidly photoionizes the diffuse circumgalactic and intergalactic gas, suppressing further gas accretion onto the galaxy. Star formation nevertheless continues for several hundred Myr using the surviving self-shielded gas reservoir and ceases only after this gas is consumed or dispersed. Within 500 Myr after reionization, less than 40

The central engines of Little Red Dots (LRDs) may be “black hole stars" (BH*s), early stages of black hole growth characterized by dense gas envelopes. So far, the most direct evidence for BH*s comes from a handful of sources where the host galaxy is completely outshone as suggested by their remarkably steep Balmer breaks. Here we present a novel scheme to disentangle BH*s from their host galaxies assuming that the [OIII]5008Å line arises exclusively from the host. Using a sample of 98 LRDs (z 2-9) with high quality NIRSpec/PRISM spectra, we demonstrate that the host-subtracted median stack displays a Balmer break >2× stronger than massive quiescent galaxies, with the rest-optical continuum resembling a blackbody-like SED (T_eff 4050 K, log(L_bol) 43.9 erg s^-1, R_eff 1300 au). We measure a steep Balmer decrement (Hα/Hβ>10) and numerous density-sensitive features (e.g., FeII, HeI, OI). These are hallmark signatures of dense gas envelopes, providing population-level evidence that BH*s indeed power LRDs. In the median LRD, BH*s account for ∼20% of the UV emission, ∼50% at the Balmer break, and ∼90% at wavelengths longer than Hα with the remainder arising from the host. BH*s preferentially reside in low-mass galaxies (M_⋆ 10^8 M_⊙) undergoing recent starbursts, as evidenced by extreme emission line EWs (e.g., [OIII]5008Å 1100Å, CIII] 12Å), thereby favoring BH* origins linked to star-formation. We show V-shaped LRD selections are biased to high BH*/host fractions (≳60% at 5500Å) – less dominant BH*s may be powering JWST's blue broad-line AGN. We find BH*s are so commonplace and transient (duty cycle ∼1%, lifetime ∼10 Myrs) that every massive black hole may have once shone as a BH*.

Observations of the redshifted 21-cm line during the Epoch of Reionization will open a new window to probe the intergalactic medium during the formation of the first stars, galaxies, and black holes. A particularly promising route to an initial detection is to cross-correlate tomographic 21-cm maps with spectroscopically confirmed Lyman-α emitters (LAEs). High-redshift LAEs preferentially reside in ionized bubbles that are strongly anticorrelated with the surrounding neutral regions traced by 21-cm observations. In this work, we study the prospect of detecting such a cross-correlation signal by stacking 21-cm image cubes around LAEs using a current-generation 21-cm instrument – the Hydrogen Epoch of Reionization Array (HERA). Our forecast adopts a realistic mapping pipeline to generate foreground-free 21-cm image cubes. The statistical properties of these images, arising from the complex instrumental response, are carefully accounted for. We further introduce a physically motivated signal template calibrated on the THESAN radiation-hydrodynamic simulations, which connects the cross-correlation amplitude to the global neutral fraction. Our results show that a sample of 50 spectroscopically confirmed LAEs is sufficient to begin constraining the reionization history. These results represent an important preparatory step toward joint analyses of 21-cm experiments with upcoming wide-area, high-redshift galaxy surveys from Euclid and the Nancy Grace Roman Space Telescope.

We use the THESAN radiation-hydrodynamics simulations to investigate how Lyman-$α$ emitters (LAEs) trace ionized bubble sizes during the Epoch of Reionization. We generate realistic LAE catalogs by combining accurate intrinsic Ly$α$ production and intergalactic transmission with an empirical model for dust absorption and gas outflows. By calibrating to observationally-constrained Ly$α$ luminosity functions, we reproduce the rapid decline in Ly$α$ visibility toward higher redshifts while revealing mild tensions in LAE fractions near the end of reionization. Before the midpoint of reionization, galaxies within larger line-of-sight bubbles ($\gtrsim 10$ cMpc) have higher observed Ly$α$ luminosity and equivalent width (EW), demonstrating that the evolving LAE fraction provides a practical statistical tracer for bubble size. These correlations weaken as percolation progresses and the IGM becomes increasingly ionized. In LAE selected samples with $L_{\text{Ly}α} > 10^{41.5}\ \text{erg s}^{-1}$, Ly$α$ properties correlate with bubble size more strongly than UV magnitude, especially at $z \gtrsim 7$. This simulation-based framework maps LAE selections to bubble-size statistics, clarifies biases in more idealized models, and will supply public catalogs to interpret current and forthcoming JWST and narrow-band LAE surveys in terms of the evolving topology of reionization.

Several Lyman-alpha (Lyα) emitters deep into the reionization era exhibit double-peaked Lyα emission profiles, raising the question of how the intergalactic medium can transmit photons blueward of the Lyα resonance at such high redshifts. To investigate this, we compute Lyα transmission along sightlines originating from galaxies in the Cosmic Dawn III simulation and identify cases that closely reproduce the observed double-peaked emission. In these cases, the sightlines intersect highly underdense voids located a few comoving megaparsecs from the source galaxy. These voids allow photons emitted blueward of Lyα to redshift through resonance without scattering while traversing them. The low opacity arises because the neutral hydrogen density scales with the square of the underlying gas density under ionization equilibrium, making sufficiently underdense regions with ≲30 % of cosmic mean density highly transmissive. Such voids naturally occur in the fluctuating cosmic density field, even in the vicinity of galaxies, and can also be associated with transmissive spikes in the Lyα forest. We find that the global probability of observing double-peaked emission is ∼0.5 % during reionization at an 80% global ionization fraction, while no cases are found at 60% ionization. We also find that this probability depends sensitively on the local ionizing background intensity, increasing by a factor of ∼10^3 for a tenfold increase in intensity. Recent accretion episodes of supermassive black holes can further increase this probability. These results suggest that double-peaked Lyα emission in high-z galaxies can serve as a sensitive probe of the ionizing background during the late stages of cosmic reionization.

The distribution of gas-phase metals within galaxies encodes the impact of stellar feedback on galactic evolution. At high-redshift, when galaxies are rapidly assembling, feedback-driven outflows and turbulence can strongly reshape radial metallicity gradients. In this work, we use the FIRE-2, SPICE, Thesan and Thesan Zoom cosmological simulations -- spanning a range of stellar feedback from bursty (time-variable) to smooth (steady) -- to investigate how these feedback modes shape gas-phase metallicity gradients at $3 10^{9}~{\rm M_\odot}$. These results demonstrate that bursty stellar feedback provides sufficient turbulence to prevent strong negative gradients from forming, while smooth stellar feedback does not generically allow for efficient radial redistribution of metals thereby keeping gradients steep. Finally, we compare with recent observations, finding that the majority -- but, notably, not all -- of the observed gradients may favor a bursty stellar feedback scenario. In all, these results highlight the utility of high-resolution observations of gas-phase metallicity at high-redshift as a key discriminator of these qualitatively different feedback types.

James Webb Space Telescope (JWST) has revealed the apparent evolution of the black hole (BH)-stellar mass (M-BH-M-*) relation in the early Universe, while remaining consistent with the BH-dynamical mass (M-BH-M-dyn) relation. We predict BH masses for z>3 galaxies in the high-resolution thesan-zoom simulations by assuming that the M-BH-M-dyn relation is fundamental. Even without live BH modelling, our approach reproduces the JWST-observed M-BH distribution, including overmassive BHs relative to the local M-BH-M-* relation. We find that M-BH-M-* declines with M-*, evolving from similar to 0.1 at M-*=10(6) M-circle dot to similar to 0.01 at M-* = 10(10.5) M-circle dot. This trend reflects the dark matter (f(DM)) and gas fractions (f(gas)), which decrease with M-* but show little redshift evolution down to z = 3, resulting in small M-*/M-dyn ratios and thus overmassive BHs in low-mass galaxies. We use prospector-derived stellar masses and star formation rates to infer f(gas) across 48 022 galaxies in the JWST Advanced Deep Extragalactic Survey at 3 < z < 9, finding excellent agreement with our simulation. Our results demonstrate that overmassive BHs would naturally result from a fundamental M-BH-M-dyn relation and be typical of the gas-rich, dark matter-dominated nature of low-mass, high-redshift galaxies. Such overmassive BHs may strongly influence early galaxy formation, and we caution that our approach does not include the self-consistent BH-galaxy co-evolution required for a complete understanding.

Patchy reionization couples the ionized-bubble morphology to the underlying density field, making the CMB Thomson optical depth sensitive to both the global ionization history and anisotropic fluctuations on the sky. Using the large-volume radiation-hydrodynamical Lumina simulation, we compute τ_ CMB in two ways: (i) from global volume- and mass-weighted ionization histories, and (ii) from explicit line-of-sight integrations through on-the-fly light cones. We find that the sightline-averaged optical depth in the light cone, ⟨ τ_ LOS⟩ = 0.0550, exceeds the value inferred from a global volume-weighted history, τ_ CMB,V = 0.0515, by ≈ 7%. This enhancement is largely captured by the global mass-weighted prediction, τ_ CMB,m = 0.0544, indicating that precision comparisons to CMB optical-depth constraints should use mass-weighted electron fractions or explicit light-cone integration rather than volume-weighted ionized fractions alone. The excess optical depth accumulates primarily near z_ LOS = 8.0^+1.9_-1.3, where the combination of high physical density and strong ionization-field patchiness is greatest. The resulting τ_ LOS field is non-Gaussian and exhibits ≳ 5% sightline-to-sightline scatter, with fluctuations tracing rare early-ionized overdensities and large-scale structure. Coarse-graining experiments show that smoothing the ionization field on ≳ 3 cMpc scales suppresses the density-ionization correlation and biases τ_ CMB low relative to the resolved calculation. Finally, angular power spectra and real-space correlation functions decomposed into HII, HeII, and HeIII auto- and cross-contributions reveal scale-dependent departures from simple hydrogen-helium co-tracing and evolving characteristic scales with redshift.

The dynamical role of Lyman-α (Lyα) radiation pressure feedback has been debated for nearly a century, with recent analytical and 1D numerical studies highlighting its potential dominance over other stellar feedback processes at Cosmic Dawn. Despite this, no multi-dimensional Lyα radiation hydrodynamics (RHD) simulations have been performed to date. In this paper, we present the first 2D Lyα RHD simulations using Lydion, an RHD code with a novel M1 moment method for Lyα transfer, and self-consistent dust dynamics. Lydion yields a ∼𝒪(100) × speed-up compared to Monte Carlo radiative transfer in simple benchmarks, making 2D Lyα RHD feasible. We perform simulations of star clusters and isolated stars embedded in dense, metal-poor (Z/Z_⊙≤ 0.01) clouds, and find that Lyα feedback dramatically boosts outflows and dominates over feedback from direct and infrared radiation pressure. Lyα leakage through lower-column density channels, Doppler shifts, and Lyα photon destruction, while important, cannot prevent the build-up of strong Lyα radiation pressure in H II regions, leading to radiative forces ∼ (2 - 16) × L_ bol/c, and Lyα force multipliers M_ F∼ 10-60. Lyα feedback may not preclude efficient star formation, but raises the threshold gas surface density for this to occur. We conclude that nearly all galaxy and star formation simulations are currently missing the strongest source of radiation pressure feedback in dense and metal-poor environments.

Characterizing the mechanisms and galaxy properties conducive to the escape of ionizing (LyC) emission is necessary to accurately model the Epoch of Reionization and identify the sources that powered it. Using Hubble Space Telescope data, the Ly α and Continuum Origins Survey (LaCOS) is the first program to obtain uniform, multiwavelength subkiloparsec imaging for a large sample (42) of galaxies observed in LyC and enable statistically robust studies between LyC and resolved galaxy properties. Here, we characterize the morphology and galaxy merger properties of LaCOS galaxies and investigate their connection with the escape fraction of LyC emission f esc LyC . We find strong anticorrelations between f esc LyC and size ( r 20 , r 50 , and r 80 ) measured in filters containing emission from star-forming regions, and with the asymmetry and clumpiness in F150LP, a filter tracing UV continuum and Ly α . We find that ≥48% of LaCOS galaxies, and ≥41% of LaCOS LyC-emitters are visually classified as galaxy mergers. Galaxies robustly identified as mergers in LaCOS are at advanced stages of interaction, close to coalescence. The f esc LyC properties of robust mergers and low-probability mergers cannot be differentiated statistically, and we only find significant difference between the two populations in terms of their sizes and LyC luminosity: robust mergers having larger values. We conclude that (i) f esc LyC tends to be larger in galaxies with a small number of compact, centrally located, UV-emitting star-forming regions, (ii) mergers at advanced stages of interaction represent a sizable fraction of LyC-emitting samples at z ∼ 0.3, and (iii) mergers can facilitate the escape of LyC photons from galaxies.

Emission lines from ionised gas are key diagnostics of star formation, metallicity, and ionisation conditions in galaxies. Interpreting spatially resolved observations from integral-field surveys (e.g. MaNGA, MUSE, JWST/NIRSpec) and comparing them with hydrodynamical simulations requires 3D photoionisation models that handle realistic geometries, dust attenuation, and synthetic instrument output. We present a new photoionisation module for the Monte Carlo radiative transfer code SKIRT that predicts emission-line luminosities of ionised gas in 3D, combining pre-computed Cloudy tables for gas temperature and opacity with a direct calculation of ion fractions and line emissivities. The local ionising radiation field (1-6 Ryd) is characterised by log U and four spectral-shape ratios; Cloudy tables map these to temperature and opacity, converging through SKIRT's existing iteration cycle. An inline solver then determines ion fractions from the converged field and temperature and evaluates line emissivities. We validate against Cloudy on 60 spherical shell models and against COLT on a Milky Way-analogue galaxy. On the 1D grid, hydrogen recombination lines agree with Cloudy to within a few per cent (Halpha median ratio 0.97) and the forbidden lines to within 5

ABSTRACT The timing of cosmic reionization across Local Group (LG) analogues provides insights into their early histories and surrounding large-scale structure. Using the radiation-hydrodynamic simulation thesan-1 and its dark-matter-only counterpart thesan-dark-1, we track the reionization histories of all haloes, including 224 LG analogues within the proximity of any of the 20 Virgo-like clusters with halo masses above $10^{14}\, \text{M}_{\odot }$ at $z = 0$ and their environments. The statistically controlled samples quantify how the reionization redshift ($z_\text{reion}$) correlates with halo mass, local overdensity, and present-day pair properties. Even at fixed mass, haloes in denser regions ionize earlier, and increasing the overdensity smoothing scale systematically suppresses small-scale structure, including local variations and environmental gradients in $z_\text{reion}$. Virgo-like clusters accelerate reionization in their surroundings out to $\sim$5–10 cMpc, beyond which local overdensity again becomes the dominant factor. Within LG pairs, reionization timing offsets reach up to $\sim$150 Myr and correlate with present-day halo separation, reflecting sensitivity to large-scale structure rather than mass ratio in driving asynchronous reionization. The results support an extreme inside-out picture where clustered sources rapidly ionize their immediate neighbourhoods, while lower density regions self-ionize later and voids wait for external homogenization. These links between environment and reionization timing explain the influence of protoclusters and help interpret fossil records in LG dwarfs around the Milky Way. For Milky Way analogues, we find a reionization redshift as early (late) as $z_{\rm reion} = 12.7^{+2.0}_{-1.7}$ ($8.88^{+0.66}_{-0.70}$) when considered on 125 ckpc (500 ckpc) scales, with LG analogues following an inside-out reionization picture.

Understanding how galaxies and active galactic nuclei (AGN) jointly drive the reionization of the intergalactic medium (IGM) across cosmic time remains a major challenge in cosmology. We present Lumina, a large-volume radiation-hydrodynamic simulation that self-consistently follows the coupled evolution of the intergalactic medium, galaxies, and AGN through HI, HeI, and HeII reionization down to redshift z=3. Lumina evolves a cosmological volume of comoving side length L_box=500 cMpc with 2× 6000^3 resolution elements, corresponding to baryonic and dark-matter mass resolutions of 3.6× 10^6 M_⊙ and 1.9× 10^7 M_⊙, respectively. The simulation uses the moving-mesh code AREPO, combining the IllustrisTNG galaxy-formation model with a GPU-accelerated M1 radiation-transport solver in six frequency bins. The initial conditions employ separate transfer functions for baryons and dark matter and include their relative streaming velocity. Lumina predicts a late, predominantly stellar-driven hydrogen reionization, with the median sub-volume fully ionized by z≈ 5.2 and residual neutral HI patches persisting until z≈ 4.75. HeII reionization is driven self-consistently by AGN and is nearly complete by z=3. The simulation yields a Thomson-scattering optical depth in excellent agreement with Planck, an IGM thermal history and photoionization background broadly consistent with observational constraints, and a clear late-time thermal boost associated with HeII reionization. Its galaxy population remains consistent with the original IllustrisTNG project, while the larger volume improves statistics for rare objects, large-scale environments, and cosmic variance, enabling forward modelling of observables linking HI and HeII topologies to the evolving galaxy and AGN populations.