On the sizes of ionized bubbles around the highest redshift galaxies. Spectral shapes of the Lyman-alpha emission from galaxies

We develop a new method to determine the distance between a high-redshift galaxy and a foreground screen of atomic hydrogen. In a partially neutral universe, and assuming spherical symmetry, this equates to the radius of a ionized 'bubble' (R_B) surrounding the galaxy. The method requires an observed Lya equivalent width, its velocity offset from systemic, and an input Lya profile for which we adopt scaled versions of the profiles observed in low-z galaxies. We demonstrate the technique in a sample of 23 galaxies at z>6, including eight at z = 7.2 - 10.6 recently observed with JWST. Our model estimates the emergent Lya properties, and the foreground distance to the absorbing IGM. We find that galaxies at z>7.5 occupy smaller bubbles (~0.5 - 1 pMpc) than those at lower-z. With a relationship that is secure at 99% confidence, we empirically demonstrate the growth of ionized regions during the reionization epoch for the first time. We independently estimate the upper limit on the Str\"omgren radii (R_S), and derive the escape fraction of ionizing photons budget necessary for reionization.


INTRODUCTION
Observations of Ly emission from galaxies have long been known to fulfill a key role in charting the history of cosmic reionizaton (see Dĳkstra 2014, for a review).Because nebular Ly emission can be redshifted from the systemic velocity (typically by a few hundred km s −1 ) it can be detected through the damping wing of the Gunn-Peterson trough (Miralda-Escudé & Rees 1998), albeit with possibly-significant attenuation.This makes the population Ly-emitting galaxies a very powerful tool to study the epoch of reionization (EoR), where an abundant Ly-emitting galaxy population produces strong line emission just redward of systemic velocity.
Purely photometric Ly measurements have been employed for EoR studies for almost two decades.The most common approaches have been to study flux deficit of Ly compared Corresponding author: Matthew J. Hayes matthew@astro.su.se with expectations.The first attempts performed differential comparison of the Ly luminosity function (LF) across redshifts (e.g.Malhotra & Rhoads 2004;Kashikawa et al. 2006), although see Dĳkstra et al. (2007) for further considerations.Developments of the technique study Ly in comparison to the UV continuum flux, either as the the 'volumetric escape fraction' (Hayes et al. 2011;Wold et al. 2017) or an evolution of the equivalent width distribution (Stark et al. 2011;Schenker et al. 2014;Cassata et al. 2015;Arrabal Haro et al. 2018).These methods, however, all rely upon ensemble of galaxies to derive one quantity, which is typically the average neutral fraction ( xHI ) at a given redshift (e.g.Mason et al. 2018).One cannot trivially derive higher order estimates of the reionization process, such as spatial variations, bubble size distribution, etc.For this kinematic/spectroscopic data are needed.
More nuanced estimates can be attained if the intrinsic Ly properties of a galaxy are known: specifically comparisons of the observed Ly EW and velocity profiles with their intrinsic values would lead directly to a measure of the distance between a galaxy and the foreground screen of absorbing H i. Little progress has been made because we need to know both the EW and velocity offset of emergent Ly, which depend on stellar conditions and radiative transfer effects in the interstellar and circumgalactic media (Verhamme et al. 2006;Dĳkstra et al. 2006;Laursen et al. 2009).Moreover, without systemic redshifts ( sys ) we cannot even begin to estimate the velocity offset with respect to that of the IGM (simply distance in an expanding universe).Mason & Gronke (2020) provide the singular exception to this, attempting to circumvent the problem using a double-peaked Ly-emitter and deriving a bubble radius of ≈ 0.7 pMpc in the  = 6.6 galaxy COLA1 (Matthee et al. 2018).
However this field is rapidly changing, and over recent years systemic redshifts have become available from other far ultraviolet emission lines like C iii] 1907, 1909Å (Stark et al. 2015(Stark et al. , 2017;;Mainali et al. 2017) and infrared lines like [C ii]158m (Pentericci et al. 2016;Carniani et al. 2017;Endsley et al. 2022b).More recently, JWST has delivered the restframe optical emission lines to provide  sys and also find Ly emission at  > 7.2 (Tang et al. 2023;Saxena et al. 2023) and even almost  = 11 (Bunker et al. 2023).Thus, part of the requirement of having systemic redshifts and Ly EWs and velocity offsets are now falling into place.
In this Letter we take advantage of the fact that Ly EWs and velocity offsets are now being measured in the EoR.We study a sample of 23 galaxies at  = 6 − 11 with systemic redshifts and Ly EWs and velocity shifts, which we present in Section 2. Their emergent Ly profiles (i.e.those that leave the galaxy after radiation transport in the ISM and CGM) are not known, but using a large sample of low- galaxies (Hayes et al. 2023) where IGM attenuation is negligible, we build realistic models for the Ly that emerges from galaxies.While there is no guarantee that the emergent Ly spectral profiles at low- match those in the EoR, we showed in Hayes et al. (2021) that we do not find evidence of their evolution in currently-available data.Using the expected damping wing from a neutral universe, we build a model for the emergent Ly observables, fitting the size of the ionized region and emergent Ly EW in a hierarchical Bayesian framework.Thus, we empirically derive the distribution of the sizes of ionized bubbles that surround galaxies across most of the reionization timeline.This method is described in Section 3 and the results in Section 4. We discuss the impact of various assumptions in Section 5, and present our concluding remarks in Section 6.Throughout we assume a cosmology of { 0 , Ω M , Ω Λ } = {70 km s −1 Mpc −1 , 0.3, 0.7}.

THE GALAXY SAMPLE
We obtain the spectral measurements for 23 galaxies at  > 6, by combining various literature samples.We include the fifteen galaxies compiled by Endsley et al. (2022b) The Ly observables of our sample.Δ red (Ly) is plotted on the ordinate axis, with EW on the abscissa.Galaxy redshift is color-coded, with the colorbar to the right.The shaded region shows the parameter space that can be matched using the emergent line profile of a low- galaxy (Hayes et al. 2023), when broadened synthetically to FWHM = 500 km s −1 .We assume the IGM density at  = 7, vary the velocity from 0 to 3000 km s −1 , and the EW between 0 and the very high value of 800 Å. the primary references to obtain uncertainties on the Ly EW and its velocity offset, Δ red (Ly).We add eight more galaxies with recently-obtained measurements from JWST.Six of these are taken from the The Cosmic Evolution Early Release Science Survey (CEERS, Finkelstein et al. 2023) that have Ly detections: four from Tang et al. (2023) at  = 7.4 − 8.7, and two from Jung et al. (2023) at  ≃ 7.5 after excluding one system for which  sys is discussed as uncertain.The final two targets are JADES targets GNz11 at  ≃ 10.6 (Bunker et al. 2023) and GS-z7-LA at  ≃ 7.3 (Saxena et al. 2023).The main properties of interest are Ly EW, velocity offset of the red peak, systemic redshifts and UV magnitudes.We show the distribution of these properties in Figure 1.
We note that this is a compilation of sources reported in different surveys, and requires spectroscopic detections in both Ly and non-resonant emission lines.The quantitative interpretation will naturally be prone to selection effects.For example, over such a broad redshift range Malmquist biases are possible, but within this small sample there is currently no systematic evolution in the average luminosity.Different emission lines also measure  sys at different redshifts, potentially impacting redshift precision where only weak lines are visible, but these uncertainties are treated within our method (Section 3.3).There is a general trend for galaxies with larger velocity shifts to show smaller EWs: this trend has been noted before, and is demonstrated in low- samples (e.g.Hayes et al. 2023) where the IGM has no influence.The relation probably arises because more massive galaxies have higher column densities of neutral gas: Ly photons must therefore undergo more scattering events in order to take longer frequency excursions to the wings of the line profile and see the gas as optically thin (e.g.Verhamme et al. 2008;Hashimoto et al. 2013).This results in larger velocity offsets and also smaller Ly escape fractions, because of the increased probability of dust absorption.
As the IGM becomes thicker with increasing redshift, and ionized regions are presumably smaller, a trend of increasing Δ red (Ly) with redshift could be expected.However, the apparent trend for galaxies with larger velocity shifts to lie at  ≳ 7.5 is not found to be significant by a two-sample Kolmogorov-Smirnoff test.

The 'Emergent' Lyman alpha profile
We use three different adjectives to describe Ly in this paper: the intrinsic properties that are produced by the H ii regions, and the observed properties that reach the telescope are commonly-used terms.Here we also use the term emergent, which refers to the Ly emitted by the galaxy (after the CGM) but before attenuation by the IGM.The emergent Ly profile is estimated from the sample of starburst galaxies at 0.05 <  < 0.44 observed with the Cosmic Origins Spectrograph on HST (Hayes et al. 2021(Hayes et al. , 2023)).At this redshift, the Ly profile is unaffected by IGM absorption.
In Hayes et al. (2023) we produced stacked average Ly spectra, where we binned the sample by various galaxy properties.Here we use the same software to estimate the emergent Ly profiles from the EoR galaxies, and show some examples in the left panels of Figure 2. We stack spectra based up measured quantities that match galaxies recently observed with JWST/NIRspec (e.g.Brinchmann 2022;Cameron et al. 2023;Tang et al. 2023), such as high [O iii]5007/[O ii]3727 line ratios (≡ O 32 ).We adopt the stack of galaxies with O 32 ratios above 10, to approximately match the measured values at  > 6 where these measurements have been made (Tang et al. 2023;Saxena et al. 2023).While the optical spectroscopic properties of the low-and high- galaxies match well, the study relies upon invariance of the Ly profiles with redshift, which we cannot test directly for EoR galaxies.However we showed in Hayes et al. (2021) that these profiles accurately reproduce those of Ly-emitters at 3 ≲  ≲ 6 when the effects of an absorbing IGM are applied.We regard this as cause for optimism, and assert that the same general profile shape can be applied at higher redshifts still.
All spectra are first continuum-subtracted, using the modeled continuum spectra described in Hayes et al. (2023) and Hayes (2023).We then shift all spectra into the restframe, using the systemic redshifts measured from optical emission lines -average stacks of these spectra are shown in the left panel of Figure 2, where each has been normalized by the luminosity in the red Ly peak before combination.However, because we are mostly interested in the spectral shape of the Ly profile at velocities redwards of line-centre, we also re-normalize our spectra onto a common frequency metric: we convert the spectra to velocity space, and rescale each to FWHM=1000 km s −1 in the red peak -these are shown in the central panel of Figure 2. The individual resampled spectra are renomalized by redshifted Ly luminosity before combination.This provides us with the average shape of the Ly emission, and in computing it we have written up a function to give the Ly profile for an arbitrary input FWHM.Thus, our software treats FWHM as a free parameter.

Absorption by the Intergalactic Medium
With a model for the emergent Ly profiles, we attenuate the spectra with a model IGM.We take the cosmic hydrogen number density scaled by a factor of (1 +  sys ) 3 at each redshift.From this we calculate the expected Ly optical depth profile  IGM () as a function of velocity, as radiation redshifts through a neutral universe.
The unknown quantity is by how much the Ly is cosmologically redshifted before it encounters the absorbing H i in the foreground.We estimate the absorption of Ly using Voigt profiles, implementing successive absorptions numerically over velocity shift, Δ IGM , which equates to a distance in an expanding universe.We multiply the emergent Ly profile by exp(− IGM ()) and calculate the zeroth and first moments of the 'observed' profile.For a given emergent EW we rescale the zeroth moment to an 'observed' EW that accounts for the IGM.I.e  obs =  emerge ( 0,obs / 0,emerge ). emerge therefore becomes a free parameter in our model, and can be compared with the EWs and velocity shifts shown in Figure 1.

Bayesian Inference on EW em , FWHM red , and 𝑅 B
For each galaxy, we perform a hierarchical Bayesian inference analysis to estimate the free parameters of our model:  em ,    red , and  B .We model the likelihood of observing the available data, i.e., the pair of { obs ±  obs , Δ( Ly ) obs ± Δ( Ly ) }, with a bivariate Normal distribution including known errors.Defining the vector of model parameters as  = { em ,    red ,  IGM }, we write the posterior for  as: where (|  ) is the prior on the model parameters, given the UV magnitude of each galaxy.Subscript 'tr' refers to Ly quantities after transmission through the IGM, which are compared to the observed quantities given a subscript 'obs'.We assume that (| UV ) =  ( em ) ( IGM ) (   red | UV ), i.e., the priors on the emergent EW, the Ly line width and the bubble radius are independent.The prior on Δ IGM is the simplest to treatwe have no empirical knowledge of this and adopt a uniform prior.For the emergent FWHM distribution we base our prior upon strong trends between  UV and the FWHM of the red peak (Hayes et al. 2023).We adopt measurements from the Lyman alpha Spectral Database (LASD; Runnholm et al. 2021)1 and the GALEX-measured UV magnitudes (see Hayes et al. 2023, Figure 21 of the online-only material) and fit a power-law between FWHM and  UV ; our prior is then Normal around this relation, and includes errors on the fit.For the emergent Ly EW we adopt an exponential distribution based upon very deep MUSE and HST observations at  = 6.We take the exponential scale length of  0 = 212 ± 186Å from Hashimoto et al. (2017).This value is already corrected for IGM absorption in that paper (following Inoue et al. 2014, same as Hayes et al. 2021), and thus should be comparable to the emergent EW required for our method.We then make the assumption that the EW distribution does not evolve strongly over the redshift of the sample.
We sample the posterior distribution using a Metropolis Hastings sampler built in pycm (Salvatier et al. 2016).Δ IGM is converted to  B using Hubble's law.We show an example fit in the right panel of Figure 2. In Figure 3 we show the posterior distribution functions for  em ,    red , and  B for each galaxy (light curves), divided into three samples according to redshift.In the discussion that follows, we use the maximum a posteriori and its associated 68% credibility interval as our best estimate and associated uncertainty on each parameter.

An Example: GNz11
We show an example model for GNz11 (Bunker et al. 2023) in the right panel of Figure 2.This is currently the highest redshift known Ly-emitting galaxy.The galaxy has an observed Ly EW of 18 Å, and Δ red (Ly) of 555 km s −1 .While it may be an AGN, our models show that its spectral properties in Ly can be recovered with relatively normal conditions.Naively, it may be considered difficult to explain a Ly emission line from  = 10.6, but in actuality it is not especially hard.
The maximum of the posterior on the Ly FWHM is 340 km s −1 -this is quite high, but within the range of widths measuerd in low- galaxies (the broadest we find in the  ∼ 0 COS sample has a FWHM=350 km s −1 in its red peak).This modeled broad Ly has a wing visible out to  ≃ 1200 km s −1 at low flux density.The maximum of the posterior probability function for the IGM velocity is almost 900 km s −1 -at this velocity offset, the damping wing draws a steep diagonal line across the emerging red Ly peak, absorbing ∼ 8 times more flux at  = 0 compared to at  = 1000 km s −1 .The result is a significant shift in the first moment of the line, to the 555 km s −1 that is observed.In doing so the Ly is also significantly suppressed, and only ∼ 7% of the emergent Ly survives IGM absorption.The emergent Ly EW is ≃ 260 Å with a credibility range of 153-800 Å at 16-84%.This value is again high for a star-forming galaxy, but not higher than observed in galaxies at effectively all redshifts.Based upon the observed H flux and Case B recombination theory, Bunker et al. (2023) estimate that around 4% of the intrinsic Ly escapes the galaxy.The is comparable to the 7 % we estimate from the emergent Ly flux; the further factor of ∼ 2 could indicate that ∼ 50 % of the intrinsic Ly escapes the galaxy after internal radiative transfer effects.

Sample-Averaged Galaxy Properties
We take the H i velocity offsets, emergent Ly EW and FWHM from the maximum values of the posterior sampling (Section 3.3) and show the pdfs in Figure 3. Faint lines show the pdfs for individual galaxies, color coded by redshift, while solid lines show the arithmetic mean pdf in each bin.Emergent EWs are typically peaked towards the lower end of the distribution near 100 Å -a value that is quite typical of star-forming galaxies and does not require extreme stellar populations or AGN.Only a handful of individual galaxies show posteriors that peak at EWs above 200 Å, one of which is obviously JADES-GS-z7-LA, with its observed EW of 400 Å.
The FWHM also peak in ranges that are typical of star forming galaxies at low-, of around 200 km s −1 , with a tail up to 600 km s −1 .Here there is a more obvious signature for galaxies at the higher redshifts to be intrinsically broader than those at  = 6−8.5,although we note that this highest redshift bin contains only three galaxies.This may be a signature of a selection bias in the requirement of Ly, where only broad emergent Ly profiles are able evade the damping wing of the IGM and be reported as a Ly detection.

The Evolution of Ionized Regions Through the EoR
After converting the IGM velocity offset to a foreground H i distance using the Hubble parameter, we show the resulting combined posterior probabilities in the right panel of Figure 3.Estimated H i foreground distances generally show flat combined posteriors for the  = 6 − 8 subsamples, which indicates a broad range of bubble sizes, from ≃ 0.5 pMpc, and extending to ≃ 2.5 Mpc.However the combined posterior of the  > 8 subsample is much more strongly peaked around distances below 1 pMpc, shows a much sharper decline with redshift, and drops to zero at 1 pMpc.
In Figure 4 we show the evolution of the distance to the foreground screen of H i (taking the peak of the posterior) as a function of redshift.Distances to the foreground H i screen range from around 0.5 to 2.5 proper Mpc at the lower redshift end.We note that our approach -and observations in general -may not be sensitive to smaller bubble radii, as a certain minimum offset will be required for Ly to be detected, and therefore be reported in the literature.There is in principle no reason why the method could not be applied to galaxies without Ly emission, but in order to avoid great degeneracies such an approach would require more informed priors on the intrinsic EW.This could perhaps be attained from optical line emission (e.g.Runnholm et al. 2020;Hayes et al. 2023).
The interpretation of the larger radii at  ≃ 6 is likely to be that the these galaxies reside in an ionized universe, and far out in the damping wing (after ∼ 2000 km s −1 ) there is little power in the IGM to modulate the Ly observables -the method breaks down for large bubbles and all Δ IGM are equally valid.The lower end of the distribution (Δ IGM ≃ 1 pMpc) are consistent estimates from Ly spectral profiles (Mason & Gronke 2020) and from ionization balance considerations (Bagley et al. 2017), both of which are estimated at comparable redshifts.
Figure 4 is striking in its absence of larger bubbles at higher redshifts.Using Kendall's  statistic, we find the trend to be significant at  = 1×10 −2 , and we note that similar credibility regions are determined at  > 8 compared with the galaxies at the lower redshift end of the sample.If taken at face value, the result implies that ionized regions increase in size from  ≃ 11 to 6 -this empirical find of ionized regions growing with time fulfills the main expectations of reioinization.The selection bias of needing Ly in emission may also enter here: smaller bubbles will not enable Ly transfer through the IGM, implying a lower limit on the bubble sizes we can capture.As this is more likely to be the case at earlier times, the strength of the relationship would only increase.
In the same figure we also show computational results from Giri & Mellema (2021).The simulated relation takes a similar form: starting from  ∼ 11, the bubble size of GNz11 falls just at the upper edge of the distribution.The simulated bubble sizes then increases towards  ≃ 7, forming a fully ionized universe (effectively infinite bubble radii) by  = 6.5.However, simulated bubbles grow somewhat faster than our data suggest.There are four observed galaxies at  ∼ 6 that seemingly have much smaller bubble radii than suggested by the simulations.If the universe is fully ionized at  = 6 on these sightlines, the inference for these galaxies could be the result of proximate Ly-absorbing systems (possibly galaxies) whose damping wings suppress and redshift Ly, and cause our model estimate smaller radii.This result may reflect the high end of the broader range of Ly optical depths observed in quasar spectra at  ≈ 6 (e.g.Bosman et al. 2022).

The Escape of Ionizing Radiation
In Section 4.3 we estimated the distance between the galaxy and the absorbing H i gas in the foreground.Here we test whether these main targeted galaxies are capable of ionizing their own H ii regions and, if so, what their properties must be.We proceed by simply investigating the size of the Strömgren sphere: Here,  0 is the intrinsic production rate of ionizing photons and  LyC esc is the ionizing escape fraction; hence the product of the two is the emitted ionizing photon rate. B is the (temperature dependent) total recombination rate coefficient under Case B, for which we assume 10 4 K.This includes three main assumptions: (1.) that we can estimate  0 from data; (2.) that galaxies remain ionizing for long enough to ionize their surrounding media; (3.) the IGM is homogeneous and of fixed density.
We assume  H to be the cosmic value at the redshift of each galaxy, as we described also for the inference in Section 3.2.Each galaxy must be strongly star-forming, as evidenced by the emergent Ly EWs of 100 Å (intrinsic values are likely higher still).Addressing assumption 1 above, we assume the ionizing photon production efficiency ( ion ) measured for GNz11 of 5.2×10 25 Hz erg −1 (Bunker et al. 2023) and convert the UV luminosity to  0 .By first setting  LyC esc to 1, we compute  max S , finding values between 0.5 and 6 Mpc.We contrast this with the distances to the foreground H i (0.5-2.5 pMpc) in the left panel of Figure 5.It is immediately obvious that, if the above assumptions hold, then 22 of the 23 galaxies have sufficient ionizing power to ionize their own bubble.
Obviously  LyC esc cannot be 1, or the galaxies would not produce the strong Ly or observed nebular line emission. LyC esc cannot be measured at these redshifts, but under the assumption that each galaxy is singularly responsible for ionizing its own H ii region, we can solve for  LyC esc in each galaxy by setting  S equal to the inferred bubble size:  LyC esc = ( Bub / S ) 3 .We show  LyC esc against the inferred Ly EW in the right panel of Figure 5.These two quantities are strongly correlated at low-and mid-redshifts and we overplot the individual data points of Flury et al. (2022b) and stacks and Steidel et al. (2018) for comparison.Our estimates span the same range as the lower redshift estimates.While this plot is suggestive that relationship could also be identifiable for galaxies in the EoR, we do not yet have sufficient samples to make similar statements.The median value of  is 5.1 %, with a 16-84 percentile range of 1-21%, which is broadly consistent with requirements for cosmic reionization (e.g.Finkelstein et al. 2019).For GNz11 we estimate  LyC esc to be just 2.4 % with a 1 confidence intervals between 1 and 6.1%.This is remarkably coincident with the value of 0.03 +0.05 −0.02 , estimated by Bunker et al. (2023) using completely different methods.One galaxy significantly outlies the distribution, with  B almost 3 times its supposed Strömgren radius, leading to an apparent escape fraction of 13.This is the extreme galaxy GS-z7-LA (Saxena et al. 2023), with Ly EW of 400Å, and demonstrating that our assumptions do not hold everywhere.

SOURCES OF UNCERTAINTY
We now discuss each of the assumptions above as a source of uncertainty, beginning with the ionizing photon production efficiency.We have assumed log( ion /Hz erg −1 ) = 25.7 throughout.In our formulation  LyC esc ∝  −1 ion , and variations of a factor of 2 in  ion will change  LyC esc by the corresponding amount.Significant reductions in  ion are hard to envisage, as they must be able to reproduce the strong emergent Ly and nebular line emission observed in the available JWST spectra.Increasing  ion by a factor of 2 may perhaps be possible (e.g.Maseda et al. 2020), but values higher than this would become inconsistent with predictions from normal stellar populations.
Secondly we have assumed the Strömgren sphere can fully form, which would require ≃ 0.6 − 7 Myr based upon the inferred range of  B , and is largest at the lower redshift end.This timescale is short compared to the typical ionizing timescales of galaxies, even at this epoch.By modeling the SEDs of galaxies at a similar epoch, (Endsley et al. 2022a) find typical stellar ages greater than this for ≃ 80 % of their sample (and also exactly the same median value of  ion we adopt above).Moreover, starburst events should not coordinate on times shorter than the freefall timescale,  ff : to bring  ff below 7 Myr for a mass of 10 8 M ⊙ would require all gas to fall from a radius of just 250 pc, which we deem unrealistic.
Next we address the assumption of fixed  H :  LyC esc is proportional to the square of this density.It is worth noting that recombination timescale of gas at cosmic  H is about twice the Hubble time at  = 6: global recombinations do not significantly affect our calculations.If  H varies substantially, it will also impact our estimates of  B , and to test this we have re-run our inference with  H rescaled from the cosmic average.Reducing  H by a factor of 2 decreases the median  B by only 16%, but increases  S by 60%; consequently  LyC esc decreases to a median value of ≃ 1%.However it is very unlikely that the vicinity of galaxies is underdense at all compared to the cosmic average.The reverse situation is more likely: doubling  H marginally increases  B , while  S decreases, causing the median  LyC esc to increase to ≃ 18%, which remains a realistic value in the EoR.In the instance where gas is clumped on smaller scales within the bubble, dense regions will experience increased recombination rates while the regions between them will have lower densities, which allows ionization fronts to propagate faster.In this case, sightline effects could also become important: if the denser regions lie in front of a galaxy, the excess absorption will push the inferred  B to larger values to compensate.Assuming the denser clumps occupy a small volume, the decreased  H elsewhere will mean the true  S is slightly larger than our estimate.In this instance  LyC esc will be underestimated, although if the denser regions lie outside of our sightline, the effects will mostly cancel.
Finally, we address the question of multiple galaxies contributing to the ionization of a single H ii region, which is indeed likely.However, it is also probable that the UVselection of these targets has found the most luminous galaxy in the vicinity.Our median value of  B = 1.1 physical Mpc corresponds to a volume of ≃ 2850 comoving Mpc 3 at  = 7, and such a volume would on average contain only 0.04 galaxies brighter than  ★ according to the recent LF of Harikane et al. (2023).In an unclustered universe, we would have to integrate to luminosities ∼ 5 times fainter for the probability of finding another galaxy within  B to reach 1.With all other quantities held constant, this would add only 20% more ionizing luminosity and not substantially changing the results.Of course galaxies do cluster, but we expect these results to hold on average in cases where the most luminous local source has been identified.We expect, however, that this is not the case for the very faint galaxy GS-z7-LA (Saxena et al. 2023), which has an observed Ly EW of ≃ 400Å, and almost certainly requires assistance from nearby ionizing sources.

CONCLUDING REMARKS
We have built a hierarchical Bayesian model to estimate the intrinsic Ly observables (emergent equivalent width and kinematic offset from systemic velocity) of galaxies in the reionization epoch, as well as the size of the ionized regions ('bubbles') in which they must reside.The model is built upon empirical Ly spectral templates observed in lower redshift galaxies where the IGM has no impact, and estimates the IGM absorption at a given redshift that best matches observation.We have applied this framework to a sample of 23 ostensibly star-forming galaxies at redshift  = 6 − 11 where systemic redshifts are available, including very recent observations from JWST.We find the following main results.
• The observed galaxies occupy ionized regions with sizes between 0.5 and 2.5 proper Mpc.The poste-rior probability distribution of the bubble size is not invariant with redshift, and is skewed towards smaller bubbles at higher redshifts.We detect an upwards evolution of the bubble size with redshift that is significant at better than 3 and demonstrates that ionized regions grow with time.The recovered bubble radii are broadly consistent with results from numerical simulations of reionization.
• From the observed UV luminosity and reported ionizing photon production efficiencies, we compute the size of the Strömgren radius of each galaxy.The Strömgren radius does not correlate with the bubble radius.We use the ratio of these two radii to estimate the escape fraction of ionizing photons, recovering a median value of 5% -this is marginally consistent with the requirements for galaxies to reionize the universe.If the IGM density within the ionized region is overdense by a factor of 2, Figure1.The Ly observables of our sample.Δ red (Ly) is plotted on the ordinate axis, with EW on the abscissa.Galaxy redshift is color-coded, with the colorbar to the right.The shaded region shows the parameter space that can be matched using the emergent line profile of a low- galaxy(Hayes et al. 2023), when broadened synthetically to FWHM = 500 km s −1 .We assume the IGM density at  = 7, vary the velocity from 0 to 3000 km s −1 , and the EW between 0 and the very high value of 800 Å.

Figure 2 .
Figure2.Left: model Ly profiles.The far left spectrum shows the stacked profiles, where each spectrum is simply renormalized in luminosity at Δ > 0 before stacking.As illustration, each colored line shows the stack of a subsample where all galaxies exceed the labeled O 32 ratio: greater than 5 in pink to greater than 20 in black.The central spectrum shows the same as the left, but the spectra are also rescaled to the same FWHM of the red peak prior to stacking.The reference FWHM is set to 1000 km s −1 (note the different abscissa range of the two plots).Solid lines are mean stacks and dotted lines are median.The far right panel shows an example profile for GNz11 before and after IGM attenuation.We show the emergent Ly profile as the grey line and the shape of the best-recovered Gunn Peterson absorption in orange.The observed Ly profile is the product of the two and is shown in black.Various characteristic velocities are labeled.

Figure 3 .
Figure 3. Binned probability distribution functions from traces of the Monte Carlo sampler.Left shows the emergent Ly EW, and center shows the FWHM of the Ly red peak.Right shows the bubble size in proper Mpc.Each posterior is color-coded by redshift, as labeled in the left plot.Hard lines show the means of the posterior distributions in each redshift bin, while faint lines show those for individual galaxies.The number of galaxies per bin is shown in the legend of the central plot.

FUVFigure 4 .
Figure 4. Evolution of the bubble radius,  B (ordinate) with redshift (abscissa).Uncertainties refer to the 16th and 84th percentiles of the posterior.Absolute UV magnitude is color-coded.The shaded regions shows the 10-90-percentile range of bubble sizes from the simulation of Giri & Mellema (2021).

Figure 5 .
Figure 5. Left: the radii of the bubble inferred from Ly modeling ( B ), contrasted with that of the Strömgren sphere ( S ).Right: the comparison of the LyC escape fraction (  LyC esc ; inferred from the left plot) with the model-inferred EW of Ly.We over-plot the values measured at  ≃ 3 by Steidel et al. (2018) and at  ≃ 0.3 by Flury et al. (2022a).Redshift is color-coded in both cases.
LyC esc increases to 18%.Our estimates of  LyC esc and Ly EW are comparable to those derived at lower-, but we do not yet recover the correlation between these quantities.