2.1. Simulations of reionisation

We model transmission through the IGM using radiative transfer simulations of reionisation run with FlexRT, the radiative transfer (RT) code described in Cain & D’Aloisio (Reference Cain and D’Aloisio2024). Our simulation setup and approach is the same as that described in Cain et al. (Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024). – we refer the reader to Section 3 of that work for details, and summarise salient aspects here. In FlexRT, the redshift evolution of the ionising photon emissivity from all galaxies, $\dot{N}_{\gamma}$ , is free to be adjusted at all redshifts to obtain a desired reionisation history and/or calibrate the simulation to match one or more observables (see also e.g. Kulkarni et al. Reference Kulkarni, Keating, Haehnelt, Bosman, Puchwein, Chardin and Aubert2019; Asthana et al. Reference Asthana, Haehnelt, Kulkarni, Aubert, Bolton and Keating2024.) We divide $\dot{N}_{\gamma}$ between halos in the simulation by assuming that the ionising output of an individual halo is proportional to its UV luminosity, $\dot{N}_{\gamma} \propto L_\mathrm{UV}$ . We assign UV luminosities to halos by abundance-matching to the UV luminosity function measured by Adams et al. (Reference Adams2024). In this work, we calibrate $\dot{N}_{\gamma}$ at $z \lesssim 7$ to match the mean transmission of the Ly $\alpha$ forest at $z \leq 6$ measured by Bosman et al. (Reference Bosman2022) (see Section 3.2 of Cain et al. Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024).

We consider three reionisation histories, all of which complete reionisation at $z \approx$ 5–5.5, as required by our calibration to the Ly $\alpha$ forest. The models differ significantly in their early stages – that is, when reionisation starts and how quickly it progresses. We use the late start/late end and early start/late end models studied in Cain et al. Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024, alongside a third model which starts reionisation even earlier than the latter. Since all our models end reionisation late, for brevity we refer to these as the late start, early start, and very early start models, respectively.Footnote ^c In Figure 1, we show the reionisation history (left), the mean Ly $\alpha$ forest transmission, $F_{\text{Ly}\alpha}$ , at $4.8 \leq z \leq 6$ (middle), and the CMB electron scattering optical depth (right). We denote the redshift of JADES-GS-z13-1-LA in the left panel. The ionised fraction in the late start, early start, and very early start at $z = 13$ is $\lt$ 1%, $\approx$ 5%, and $\approx$ 25%. These scenarios differ considerably in their Ly $\alpha$ transmission properties around star-forming halos at $z \gtrsim 10$ , despite all being calibrated to match the same Ly $\alpha$ forest transmission properties at low redshift.

Figure 1.

Properties of the reionisation models used in this work. Left: the volume-weighted mean ionised fraction vs. redshift (reionisation history), with the redshift of JADES-GS-z13-1-LA indicated. The late start, early start, and very early start models have ionised fractions of $\lt 1\%$ , $\approx 5\%$ , and $\approx 25\%$ at $z = 13$ , respectively. Middle: mean transmission of the Ly $\alpha$ forest at $4.8 \lt z \lt 6$ compared to measurements from Becker & Bolton (2013) and Bosman et al. (2022). Right: CMB electron scattering optical depth. The late start model is slightly more than $1\sigma$ below the Planck measurement, while the early start case is within $1\sigma$ of the fiducial measurement and the re-analysis by de Belsunce et al. (2021). The very early start case is in more than $4\sigma$ tension with the fiducial Planck result.

2.2. Modelling Ly $\alpha$ transmission statistics

We follow the procedure described in Section 3.3 of Cain et al. (Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024) to model Ly $\alpha$ transmission through the IGM on the red side of line systemic. We run an Eulerian hydrodynamics simulation of the IGM using the RadHydro code of Trac & Pen (Reference Trac and Pen2004, Reference Trac and Pen2006) with the same initial conditions used in the dark matter N-body simulation used to generate the halos used in the RT simulations. This run has a box size of 200 $h^{-1}$ Mpc and $N = 2048^3$ uni-grid gas cells, for a spatial resolution of $\Delta x = 97$ $h^{-1}$ kpc. While this resolution is insufficient to capture the physics at play in setting the properties of the emerging Ly $\alpha$ line and absorption by the CGM of galaxies, it is enough to capture the effect of large-scale gravitational inflows around halos and local density fluctuations in the IGM around galaxies. The former play a particularly important role in setting Ly $\alpha$ transmission near the line centre (Park et al. Reference Park2021).

We model the Ly $\alpha$ absorption line profile using the analytic approximation given in Tepper-García (Reference Tepper-Garca2006). At $z = 13$ , we trace 50 randomly oriented sightlines around halos with $-19 \lt M_\mathrm{UV}\lt-18$ , for a total of $\sim$ 50 000 sightlines. We integrate a distance of 200 $h^{-1}$ Mpc away from the halos, sufficient to converge on the damping wing absorption on the red side of systemic in the neutral IGM. To avoid including absorption arising from within the un-resolved halos themselves, we set the start of each sightline 500 $h^{-1}$ kpc away from the halo centre. Because the spatial resolution of the simulation is too low to capture the integration over the central line profile in regions with high inflow velocities, we artificially boost the spatial resolution in the line integration by a factor of $4\times$ , and use a cloud-in-cell scheme to interpolate grid quantities at intermediate pointsFootnote ^d (as described in Cain et al. Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024; Gangolli et al. Reference Gangolli, D’Aloisio, Cain, Becker and Christenson2024). We find this procedure produces converged transmission spectra in nearly all cases.

In Figure 2, we show the median IGM transmission ( $T_\mathrm{IGM}$ ) as a function of wavelength ( $\lambda$ ) for each of our reionisation models, using the same line styles as Figure 1. The shaded regions indicate $1\sigma$ sightline-to-sightline scatter. The vertical dashed line denotes the systemic Ly $\alpha$ wavelength. We see a marked difference between models, especially close to systemic. The late start model has $T_{\mathrm{ IGM}} \lt 20\%$ within a few Å of systemic, and declines gradually to 0 at line centre. In the early start and very early start models, $T_\mathrm{IGM}$ jumps to 30% and 60% close to systemic, respectively, reflecting the presence of ionised bubbles around most halos in those models. We thus expect required intrinsic Ly $\alpha$ brightness to produce a given observed brightness to differ considerably between these scenarios, especially for lines emitted close to systemic.

Figure 2.

Median Ly $\alpha$ transmission as a function of wavelength in each of our reionisation models. The colours and line styles are the same as in Figure 1, and the shaded regions denote $1\sigma$ sightline-to-sightline scatter. We see a substantial difference in IGM damping wing transmission between reionisation models, especially close to systemic Ly $\alpha$ .

2.3. Spectrum model

The available NIRSpec/PRISM spectrum for JADES-GS-z13-1-LA shows only one clear emission line – Ly $\alpha$ – with no other lines conclusively detected. SED fitting of JADES-GS-z13-1-LA indicates a young, metal-poor stellar population; however, standard SED fitting codes cannot model the attenuated Ly $\alpha$ emission line and strong Ly $\alpha$ damping wing that we observe in this object’s spectrum, so we use the phenomenological model described in Witstok et al. (Reference Witstok2025).

The spectrum of JADES-GS-z13-1-LA also has an unusual UV continuum turnover around $\lambda=1\,335$ Å. We model this using a DLA, as discussed above. However, if there were a DLA in front of the Ly $\alpha$ emission source, we expect that the DLA would completely attenuate the Ly $\alpha$ emission line, so the Ly $\alpha$ emission source cannot be in the same place as the UV continuum source. Witstok et al. (Reference Witstok2025) suggests two morphological models of JADES-GS-z13-1-LA in which the continuum is attenuated by a DLA, but the emission line is not: Ly $\alpha$ photons escape either through ionisation cones or diffusion through an inhomogeneous ISM. It is also possible that the Ly $\alpha$ photons diffused directly through the DLA and were redshifted in the process, as described in Dijkstra et al. (Reference Dijkstra2014), but our calculations indicate that this would redshift the Ly $\alpha$ line significantly more than we observe.

Following Witstok et al. (Reference Witstok2025), we model the spectrum as a Gaussian intrinsic Ly $\alpha$ emission line plus a power law UV continuum,Footnote ^e with the continuum subject to DLA absorption and the whole spectrum subject to IGM absorption. The DLA transmission, $T_{\text{DLA}}(\lambda)$ , is modelled using the analytic fit prescribed in Tepper-García (Reference Tepper-Garca2006) with the HI column density $N_\textit{HI}$ as a free parameter. The IGM transmission, $T_{\mathrm{IGM}}(\lambda)$ , is extracted from the simulation as described in Section 2.2. We parameterise the continuum flux using two free parameters: a characteristic wavelength $\lambda^*$ and the UV slope $\beta_\mathrm{ UV}$ :

(1) \begin{equation}F_{\text{cont}}(\lambda) = \Big(\frac{\lambda}{\lambda^*}\Big)^{\beta_\mathrm{ UV}}\,10^{-21}\,\text{erg}\,\text{s}^{-1} \text{cm}^{-2}\,{\unicode{x00C5}}^{-1}.\end{equation}

We parameterise the intrinsic emission line using three free parameters: the velocity offset $\Delta{v}$ , the intrinsic equivalent width EW, and the velocity width $\sigma$ . We set the line amplitude so that the line has equivalent width EW. Using these variables, the line flux $F_{\text{Ly}\alpha}(\lambda)$ is given by a Gaussian profile centred at $\lambda_{\text{Ly}\alpha}+\Delta\lambda$ with standard deviation $\sigma_\lambda$ , where $\Delta \lambda$ , $\sigma_\lambda$ are the wavelength-space equivalents of $\Delta v$ , $\sigma$ , respectively. The observed spectrum is modelled by:

(2) \begin{equation}F_{\text{obs}}((z+1)\lambda) = \big[F_{\text{Ly$\alpha$}}(\lambda) + F_{\text{cont}}(\lambda)T_{\text{DLA}}(\lambda)\big]T_\mathrm{IGM}(\lambda).\end{equation}

where the redshift z is taken as a free parameter since we have no other emission lines with which to precisely constrain it. The observed spectrum is convolved with the JWST NIRSpec PRISM broadening kernel as reported in Marshall et al. (Reference Marshall2025) to obtain the final modelled observations. We use the dynamic nested sampler from the dynesty package to compute posterior distributions of parameters using the random walk methodFootnote ^f (Skilling Reference Skilling, Fischer, Preuss and Von Toussaint2004; Speagle Reference Speagle2020). We assume a standard Gaussian likelihood to handle measured uncertainties in the observed spectrum.

In Figure 3, we show the spectrum of JADES-GS-z13-1-LA (blue line), including $1\sigma$ uncertainties denoted by the shaded regions. The peak on the left is the Ly $\alpha$ emission line, and the recovery of the flux to the right is well-fit by a damping wing profile (see annotations). The black curve shows the maximum-likelihood fit assuming the early start reionisation model. The bottom panel shows the residuals of the fit. We see that, as in Witstok et al. Reference Witstok2025, Equations (1)–(2) capture well the essential features of the spectrum near Ly $\alpha$ . Note that the UV continuum parameters ( $\lambda^*$ , $\beta_\mathrm{UV}$ , $N_{\text{HI}}$ ) are constrained by spectral features significantly redward of Ly $\alpha$ systemic. The power-law continuum parameters $\lambda^*$ , ${\beta}_{UV}$ are constrained by flux at wavelengths with $\lambda\gtrsim 1\,335$ Å redward of the damping wing absorption, and the DLA column density is constrained by flux at wavelengths with $1\, 250\,$ Å $\lesssim \lambda \lesssim 1\, 335\,$ Å, where the DLA absorption dominates over that from the IGM. Therefore, they do not vary significantly when we change the IGM transmission curve.

Figure 3.

Top: An example fit of the spectrum model for the early start model, obtained from the median of the highest-likelihood fits from each sightline’s dynesty posterior distribution for the early start reionisation model. The blue histogram indicates the PRISM SED from Witstok et al. (Reference Witstok2025), to which we fit the model in our analysis. The shaded regions indicate $1\sigma$ uncertainties for each spectral bin. The black line indicates the model, and the red dashed line indicates the model with an additional CII^* emission line with $\text{EW}_{\text{CII*}}=15$ Å. Bottom: residuals from the fit. The shaded regions indicate the same $1\sigma$ uncertainties as above. The red dashed histogram indicates residuals when the potential CII^* line is included in the model. Note that the residuals are consistent with continuum fluctuations except around $\lambda=1\,335$ Å. Including the CII^* emission line smooths out the residuals and makes them consistent with continuum fluctuations.

One notable deviation between the spectrum of JADES-GS-z13-1-LA and our fiducial best-fit is the elevated flux at the edge of the damping wing feature, around $\lambda = 1\,335$ Å. Because its width is exactly what we would expect from a sharp emission feature of negligible spread broadened by the NIRSpec PRISM, we posit that some of the flux may be from a CII^* $\lambda1335$ emission line.Footnote ^g The red-dashed curve is a second fit to the spectrum that includes this line, and we see that it fits this spectral feature well. The inferred rest-frame EW for the CII^* line is 15 Å, with a velocity offset $\Delta v_{\text{CII*}}$ that varies somewhat with reionisation history.Footnote ^h We find that although including a CII^* $\lambda1335$ emission line produces somewhat higher likelihoods, it only negligibly influences the other parameters, so the CII^* line can be safely ignored when modelling the Ly $\alpha$ emission line and DLA damping wing feature. An important caveat is that that the PRISM spectrum of JADES-GS-z13-1-LA does not similarly suggest possible [CII] $\lambda2326$ or CIII] $\lambda1909$ emission lines, which we would expect to see if we observe a CII^* $\lambda1335$ emission line in an AGN spectrum (Moy Rocca-Volmerange Reference Moy and Rocca-Volmerange2002; Humphrey Binette Reference Humphrey and Binette2014). Further, we estimate the peak signal-to-noise of the inferred line (if it is there) to be close to unity, so we do not claim a statistically significant detection.

Table 1.

Priors and posteriors for all parameters from dynesty modelling of JADES-GS-z13-1-LA. We give results for all three of our reionisation histories. Each prior is a uniform prior over the given interval. Note that only $\Delta v$ and EW change significantly with reionisation history.

4.1. The observed distribution of LAEs

We use the Fiducial EW and $\Delta v$ distributions given in Equations (3)–(4) of Qin & Wyithe (Reference Qin and Wyithe2025) to model the PDFs of EW and $\Delta v$ for the $z = 13$ galaxy population, which we denote as $p_{\text{obs}}(\text{EW})$ and $p_{\text{obs}}(\Delta v)$ . These add redshift evolution to the $z\sim6$ distributions given in Mason et al. (Reference Mason, Treu, Dijkstra, Mesinger, Trenti, Pentericci, de Barros and Vanzella2018) upon which they are based. The original distribution is derived from the De Barros et al. (Reference De Barros2017) and Pentericci et al. (Reference Mason, Treu, Dijkstra, Mesinger, Trenti, Pentericci, de Barros and Vanzella2018) samples, and the redshift-evolved distribution includes samples from Witstok et al. (Reference Witstok2024) and Tang et al. (Reference Tang, Stark, Topping, Mason and Ellis2024).Footnote ^j For the purpose of modelling probabilities, we use $M_\mathrm{ UV}=-18.7$ (Witstok et al. Reference Witstok2025a) and a dark matter halo mass of $1.06\times 10^{10}$ M $_{\odot}$ , which we derive from our simulation’s mass-luminosity relations, as derived via abundance matching (see Section 2.1). The former determines the EW distribution and the latter determines the $\Delta v$ distribution (Qin & Wyithe Reference Qin and Wyithe2025).

Setting $p(\text{EW},\Delta v|\text{SFG})=p_{\text{obs}}(\text{EW})p_{\text{obs}}(\Delta v)$ assumes that the EW and $\Delta v$ are statistically independent. In reality, they are likely correlated (Erb et al. Reference Erb2014), albeit with significant intrinsic scatter. Part of this correlation is accounted for in the assumed $M_\mathrm{UV}$ dependence of the intrinsic PDF (see Mason et al. Reference Mason, Treu, Dijkstra, Mesinger, Trenti, Pentericci, de Barros and Vanzella2018a). Since our distributions are derived from galaxies in a small UV magnitude range, we do not expect $M_\mathrm{UV}$ dependence to introduce a significant correlation. However, some additional correlation may arise from Ly $\alpha$ RT effects within galaxies, which we are not able to account for here. As in Qin & Wyithe (Reference Qin and Wyithe2025), we use our PDF of observed LAE properties to derive a PDF of intrinsic LAE properties by accounting for IGM attenuation of the $z\geq6$ LAE sample (from Tang et al. Reference Tang, Stark, Topping, Mason and Ellis2024; Witstok et al. Reference Witstok2024; and references therein), from which our observational EW distribution is derived. To this end, we use the (properly normalised) distribution of intrinsic EW and $\Delta v$ for LAEs:

(5) \begin{equation} p(\text{EW},\Delta v|\text{SFG}) \propto p_{\text{obs}}(\text{EW}\cdot \mathcal{T}_6(\Delta v))p_{\text{obs}}(\Delta v).\end{equation}

Here we introduce the attenuation factor $\mathcal{T}_6(\Delta v)$ , which denotes the ratio between the observed EW and the intrinsic EW of a Ly $\alpha$ emission line with velocity offset $\Delta v$ emitted at $z\sim6$ . Multiplying the intrinsic EW of an emission line by $\mathcal{T}_6(\Delta v)$ gives the corresponding observed EW. We compute $\mathcal{T}_6(\Delta v)$ from our simulations, so it additionally depends on the reionisation history. As a necessary simplification, we compute the function $\mathcal{T}_6$ using the average transmission curve for each reionisation history, and use the same function $\mathcal{T}_6$ for all sightlines in each reionisation history. Most sightlines at $z\sim 6$ have an attenuation factor of order unity, so this is not significantly different from computing attenuation factors for each individual sightline. Note that Equation (5) breaks the assumption of independence between the EW and $\Delta v$ PDFs.

4.2. Modeling an AGN distribution

While Ly $\alpha$ EW has been measured for statistical samples of quasars at $z \sim 5-6$ (Bañados et al. Reference Asthana, Haehnelt, Kulkarni, Aubert, Bolton and Keating2016; Gloudemans et al. Reference Gloudemans2022), these objects are typically 5–10 magnitudes brighter than JADES-GS-z13-1-LA. As such, we cannot directly estimate the intrinsic Ly $\alpha$ emission properties for AGN of comparable brightness, as we do for SFGs. Therefore, to estimate $p(\text{EW},\Delta v|\text{AGN})$ , we assume that the shape of the EW distribution is the same for AGNs as for SFGs, but shifted higher by a factor proportional to the ratio of the mean ionising photon production rates of AGN and SFGs. This implicitly assumes that the gas in the ISM/CGM of galaxies being in photo-ionisation equilibrium, such that Ly $\alpha$ emission from recombinations is proportional to the local absorption of ionising photons. It also makes the assumption that the radiative transfer physics of an AGN host galaxy is identical to that of an SFG. A lack of data on AGN-powered LAEs forces us to use this simplified model, which is a simple way to implement the central assumption that AGNs produce more ionising photons than SFGs.

Using this assumption, we estimate an EW distribution for AGN-powered LAEs by horizontally rescaling the SFG EW distribution by a factor of the ratio

\begin{align*}r_{\xi_{\text{ion}}}=\xi_{\text{ion,AGN}}/\xi_{\text{ion,LAE}}=2.88, \end{align*}

where

\begin{align*}\xi_{\text{ion,AGN}}\approx10^{25.9}\,\text{Hz}\,\text{erg}^{-1}\end{align*}

is the average ionising production efficiency of AGNs (Matsuoka et al. Reference Matsuoka2018) and

\begin{align*}\xi_{\text{ion,LAE}}\approx10^{25.4}\,\text{Hz}\,\text{erg}^{-1}\end{align*}

is the average ionising production efficiency of LAEs (Simmonds et al. Reference Simmonds2023; see also Saxena et al. Reference Saxena2024) This amounts to setting

(6) \begin{equation} p(\text{EW},\Delta v|\text{AGN})\propto p(\text{EW}/r_{\xi_{\text{ion}}},\Delta v|\text{SFG}).\end{equation}

JADES-GS-z13-1-LA is inferred to have

\begin{align*}\xi_{\text{ion}}\approx10^{26.5}\,\text{Hz}\,\text{erg}^{-1},\end{align*}

significantly greater than the average for either AGNs or LAEs (Witstok et al. Reference Witstok2025a). However, because AGNs generally have greater $\xi_{\text{ion}}$ than ordinary LAEs, this high $\xi_{\text{ion}}$ favours the AGN hypothesis. This approach quantifies our earlier suggestion that the extreme inferred properties of JADES-GS-z13-1-LA may indicate AGN activity.

4.3 Results

As our prior probability, $P(\text{AGN})$ , we use the overall AGN fraction of galaxies in the JADES survey, $f_{\text{AGN}}=20 \pm 5\%$ , as reported in Scholtz et al. (Reference Scholtz2025). This is the probability that a randomly chosen galaxy hosts an AGN, so it is the probability that JADES-GS-z13-1-LA hosts an AGN before accounting for its emission properties. $P(\text{SFG})$ is simply given by $1-P(\text{AGN})$ .

The dependence of the probability on IGM neutral fraction is shown in Figure 5. The large black points show our fiducial values $P(\text{AGN}|D)$ , and the magenta line denotes the boundary between a preference for the AGN and SFG scenarios. To show the importance of our assumptions in estimating these probabilities, we also plot results with varying values of $r_{\xi_{\text{ion}}}$ and $f_{\text{AGN}}$ in Figure 5. A lower value of $r_{\xi_{\text{ion}}}=2.00$ , which is suggested by some results which find a stronger upward redshift evolution of $\xi_{\text{ion}}$ for SFGs (Papovich et al. Reference Papovich2025), decreases the estimated $P(\text{AGN}|D)$ by more than 10%. However, even with lower values of $r_{\xi_{\text{ion}}}$ and $f_{\text{AGN}}$ , our main conclusions hold. In particular, the late start model somewhat favours the AGN scenario even with the most conservative parameters (the blue downward triangles), and the very early start model strongly favours the SFG scenario even with the most optimistic parameters (the red upward triangles). The early start model has the most uncertainty, probably because (i) it is at a transitional point of reionisation where sightline to sightline (and halo to halo) variation is particularly high, and (ii) the probability is further away from the minimum (prior)Footnote ^k and maximum (100%) values that it can attain, so the distribution of probabilities is wider.

Figure 5.

Values of $P(\text{AGN}|D)$ , including uncertainties, for varying values of $r_{\xi_{\text{ion}}}$ and $f_{\text{AGN}}$ . Our fiducial value, for $r_{\xi_{\text{ion}}}=2.88$ and $f_{\text{AGN}}=0.20$ , is shown in black. We compare this to probabilities obtained with a lower $r_{\xi_{\text{ion}}}$ of $2.00$ to show how $r_{\xi_{\text{ion}}}$ affects the estimate. We also compare this to probabilities obtained with priors at the upper and lower 1 $\sigma$ limits from Scholtz et al. (Reference Scholtz2025). These priors, along with our fiducial prior, are shown on the plot as dotted grey lines.

Owing to the dearth of available data on high-redshift AGN-driven LAE populations and the general difficulty of detecting Type 2 AGNs, our result can only be considered a rough estimate. Calhau et al. (Reference Calhau2020) notes that AGNs are systematically under-detected because obscured AGNs are often difficult to detect in the radio and X-ray bands. Therefore, since our calculation is based on empirically derived AGN fractions, it may be an underestimate.

If JADES-GS-z13-1-LA is not an AGN, our results (at face value) favour an ionised fraction in excess of a couple percent at $z = 13$ . Recent observations tend to support a scenario in which the bulk of reionisation started relatively late at $z \lt 10$ (Cain et al. Reference Cain, Lopez, D’Aloisio, Munoz, Jansen, Windhorst and Gangolli2024), disfavouring our very early start scenario. However, reionisation could still have a long, extended tail well above $z = 10$ (Park et al. Reference Park, Shapiro, Komatsu, Iliev, Ahn and Mellema2013), or that the early stages of reionisation may have been non-monotonic, with a ‘bump’ of ionisation of high redshifts (Cen Reference Cen2003; Ahn & Shapiro Reference Ahn and Shapiro2021). These types of scenarios often involve some early contribution to reionisation, at the 5–10% level, from Pop. III stars (Tan Reference Tan2025), and/or contributions to reionisation from mini-halos (Norman et al. Reference Norman, Chen, Wise and Xu2018). At face value, such early reionisation scenarios are generally disfavoured by constraints on the electron scattering optical depth $\tau_\textit{CMB}$ from the low- $\ell$ polarisation of the CMB measured by Planck (Wu et al. Reference Wu, McQuinn, Eisenstein and Iršic2021b). However, it has been recently noted that the latest DESI-DR2 data (DESI Collaboration Reference Collaboration2025) may prefer a higher value of $\tau_\textit{CMB}$ in combination with the CMBFootnote ^l (Sailer et al. Reference Sailer, Farren, Ferraro and White2025; Jhaveri, Karwal, & Hu Reference Jhaveri, Karwal and Hu2025) (although see Cain et al. Reference Cain, Engelen, Croker, Kramer, D’Aloisio and Lopez2025) for difficulties with this scenario). Forthcoming CMB experiments such as LiteBIRD (Matsumura et al. Reference Matsumura2014) will help resolve this question.

4.4. Corroborating the AGN scenario

A few additional lines of evidence support the possibility that JADES-GS-z13-1-LA hosts an AGN. On their own, they are tentative evidence at best, but together with the relatively high estimated probability of the AGN scenario in both of our more plausible reionisation histories, they collectively paint a compelling picture of this possibility.

As mentioned in Section 2.3 and discussed in Witstok et al. (Reference Witstok2025), the UV turnover in the JADES-GS-z13-1-LA spectrum is best explained by a DLA with HI column density $\sim10^{23}\,\text{cm}^{-2}$ . This is consistent with the column density of the torus that surrounds an AGN (Peca et al. Reference Peca2021). This would account for both the obscuration of the AGN and the damping wing feature. AGN obscured by dense HI in their torus often display significant variability in HI column density on short timescales (Liu et al. Reference Liu2017; Cox et al. Reference Cox, Torres-Albà, Marchesi, Gianolli, Zhao, Ajello, Pal and Silver2025). If follow-up observations of JADES-GS-z13-1-LA detected clear evidence of changes in $N_\textit{HI}$ , this would be compelling evidence for the torus obscuration scenario.

In Section 2.3, we also report a tentative, potential detection of a CII^* doublet emission line at $\lambda=1\,334,1\,335$ Å (Figure 3). This potential line was not mentioned in Witstok et al. (Reference Witstok2025), perhaps because it coincides with the peak of the DLA feature. This emission line is often detected in quasar spectra (Vanden Berk et al. Reference Vanden Berk2001; Maiolino et al. Reference Maiolino2024). On the other hand, this line is rarely detected in SFG spectra, and when it is detected in these spectra it is usually an absorption feature (Berg et al. Reference Berg2022). If higher-resolution followup spectroscopy confirms the presence of a strong CII^* emission line, especially one without a resonant absorption feature, it would be compelling evidence that this galaxy hosts an AGN.

For JADES-GS-z13-1-LA, Witstok et al. (Reference Witstok2025) reports a Ly $\alpha$ luminosity of $\sim$ $10^{43.4}\,\text{erg}\,\text{s}^{-1}$ , and we infer that the Ly $\alpha$ luminosity is $\gtrsim$ $10^{43}\,\text{erg}\,\text{s}^{-1}$ . At this luminosity, a relatively high fraction of LAEs host AGNs (Konno et al. Reference Konno, Ouchi, Nakajima, Duval, Kusakabe, Ono and Shimasaku2016; Calhau et al. Reference Calhau2020), although this fraction has not been determined beyond $z\sim6$ . Additionally, Baek & Ferrara (Reference Konno, Ouchi, Nakajima, Duval, Kusakabe, Ono and Shimasaku2013) find that, in dark matter halos of masses near $\sim$ $10^{10}$ M $_{\odot}$ , SFGs cannot produce the Ly $\alpha$ emission we observe for this galaxy, while Compton-thin AGNs (that is, AGNs with obscuring column densities of $\sim$ $10^{22-24}\,\text{cm}^{-2}$ ) can. This statement is related to our previous argument about emission features, but qualitatively corroborates it with independent results that link extreme Ly $\alpha$ emission and AGN activity.