Radiative transfer of 21-cm line through ionised cavities in an expanding universe

The optical depth parameterisation is typically used to study the 21-cm signals associated with the properties of the neutral hydrogen (HI) gas and the ionisation morphology during the Epoch of Reionisation (EoR), without solving the radiative transfer equation. To assess the uncertainties resulting from this simplification, we conduct explicit radiative transfer calculations using the cosmological 21-cm radiative transfer (C21LRT) code and examine the imprints of ionisation structures on the 21-cm spectrum. We consider a globally averaged reionisation history and implement fully ionised cavities (HII bubbles) of diameters $d$ ranging from 0.01 Mpc to 10 Mpc at epochs within the emission and the absorption regimes of the 21-cm global signal. The single-ray C21LRT calculations show that the shape of the imprinted spectral features are primarily determined by $d$ and the 21-cm line profile, which is parametrised by the turbulent velocity of the HI gas. It reveals the spectral features tied to the transition from ionised to neutral regions that calculations based on the optical depth parametrisation were unable to capture. We also present analytical approximations of the calculated spectral features of the HII bubbles. The multiple-ray calculations show that the apparent shape of a HII bubble (of $d=5$ Mpc at $z=8$), because of the finite speed of light, differs depending on whether the bubble's ionisation front is stationary or expanding. Our study shows the necessity of properly accounting for the effects of line-continuum interaction, line broadening and cosmological expansion to correctly predict the EoR 21-cm signals.


INTRODUCTION
At about 0.3 Myr after the Big Bang, electrons and protons began to combine to form neutral hydrogen (H I) atoms (Planck Collaboration et al. 2014).This allows the photons to decouple from the charged particles and form the relic radiations which then evolved to become the Cosmological Microwave Background (CMB) observed today.The Universe entered a dark age before the first luminous objects began to form.The UV photons from the first stars and the X-rays from the first active galactic nuclei (AGN) are capable of ionising H I gas, carving out ionised hydrogen (H II) cavities (see e.g.Loeb & Barkana 2001).As these ionised cavities expand (Furlanetto et al. 2004(Furlanetto et al. , 2006a) ) and merge (Iliev et al. 2006;Robertson et al. 2010), the Universe eventually became almost completely reionised, except in some dense self-shielded regions (see e.g.Chardin et al. 2018).This Epoch of Reionisation (EoR) is believed to start at  ∼14 and ★ E-mail: kinwah.wu@ucl.ac.uk (KW), qin.han.21@ucl.ac.uk (QH), jyhchan@cita.utoronto.ca(JYHC) complete at  ∼6 based on observations of high redshift luminous quasars (Fan et al. 2006;Robertson et al. 2010) and large-scale polarisation of the CMB (Planck Collaboration et al. 2016b,c).
These observational constraints, however, suffer from various shortcomings, making it difficult for us to obtain a detailed history of reionisation (see e.g.Natarajan 2014;Mesinger 2016).It can be compensated by observing the 21-cm signals produced during the EoR (Morales & Wyithe 2010;Baek et al. 2010;Eide et al. 2018;Pritchard & Loeb 2012;Eide et al. 2020).The 21-cm signals are generated by the spin-flip of H I atoms in the ground state, and corresponds to a frequency of ν 21cm = 1.42 GHz (Hellwig et al. 1970;Essen et al. 1971).Due to the abundance of H I gas of the Universe, we will be able to access the full history of the reionisation via the 21-cm signals when we overcome the difficulties in observation and data analysis.Several projects including the Low Frequency Array (LOFAR, Zarka et al. 2012;Mertens et al. 2020), Murchison Widefield Array (MWA, Tingay et al. 2013;Trott et al. 2020) the Hydrogen Epoch of Reionization Array (HERA, DeBoer et al. 2017;HERA Collaboration 2022) and the Square Kilometer Array (SKA, Koopmans et al. 2015) are underway.
Currently, most of the studies focus on simulating luminous sources and their contributions to reionisation (Santos et al. 2010;Mesinger et al. 2011;Xu et al. 2016;Park et al. 2019;Mangena et al. 2020;Gillet et al. 2021;Doussot & Semelin 2022).The calculation of 21-cm signals are merely analytical post-processing of simulation results based on the optical depth parametrisation (Furlanetto et al. 2006b;Pritchard & Loeb 2012).Adopting this parametrisation implies that the observed 21-cm signal at  = 0 at a specific frequency ν 0 unambiguously reflects the properties of H I gas at redshift  = ν 21cm /ν 0 − 1 along the line-of-sight (LoS).Line broadening and radiative transfer effects are not correctly accounted for in this parametrisation (see also Chapman & Santos 2019).
To quantify the inaccuracies of the 21 cm signals computed with the optical depth parameterisation, we adopt a covariant radiative transfer formalism which does not require various assumptions which were used in the derivation of the 21-cm optical depth.We investigate the spectral features imprinted by ionised cavities on the 21-cm spectra at  = 0.As most simulations study reionisation on scales larger than galaxies (∼ 1 kpc) and present the statistical properties of reionisation, we do not use their detailed ionisation structures.Instead, we construct more generic descriptions for evolving ionised cavities in an expanding Universe which is applicable for all length scales and various ionising sources, in Section 2. To facilitate the comparison between the covariant method and the optical depth paramertrisation method, we use the most simplified fully ionised spherical cavities as input ionisation structure and calculate their imprints on the 21-cm spectra at  = 0. We use C21LRT (cosmological 21-cm line radiative transfer) code which tracks the 21 cm signal in both frequency and redshift space, and takes full account of cosmological expansion, radiative transfer effects, and possible dynamics and kinetic effects, such as those caused by macroscopic bulk flow and microscopic turbulence (Chan et al. 2023).The C21LRT formulation and the computational set-up of radiative transfer calculaitions are presented in Section 3. We then analyse the spectral features of fully ionised cavities in the absorption and emission regimes of 21-cm global signals and discuss the implications in Section 4. Finally, we summarise the findings in Section 5.

IONISED CAVITIES AND TIME EVOLUTION
The physics of ionisation processes and the dynamic models of ionisation-front expansion can be found in Axford (1961) (see also, Newman & Axford (1968); Yorke (1986); Franco et al. (2000)).For the purposes of this study, a generic description for the ionisation bubbles and their expansion is sufficient.However, we summarise all the physical processes pedagogically for completeness of this paper.We then discuss whether the corresponding physics are incorporated in reionisation simulations and the implications for the accuracy of the predicted 21 cm signals.

Model ionised cavities
We only consider H atoms (H I and H II) in our cavity models and C21LRT calculations without any other species.In the following texts, we use 'H I zone' and 'H II zone' for regions filled with only H I or H II gas, respectively, 'partially ionised zone' for regions with both H I and H II gas, and 'ionised cavities' refers to more general/realistic ionisation structures.The H II zones are assumed to be spherical.Although they would not be spherical in reality, adopting this idealised geometry is sufficient for the purpose of this study, that is, to demonstrate the development of the 21-cm line propagating through evolving ionised cavities.which are driven by continuing ionisation together with cosmological expansion.
For a H II zone surrounded by an infinitely extended H I zone, its size increase are due to (i) photoionisation, which is caused by the radiation emitted from the sources embedded inside the H II zone, and (ii) length stretching as a consequence of cosmological expansion.We consider first a two-zone model, which has an insignificant thickness of the transitional interface between the H II zone and the H I zone.We next consider a three-zone model.It has a H II zone surrounded by infinitely extended H I zone, and a layer of partially ionised gas between the H I zone and H II zone.This partially ionised zone has sufficient contribution to the 21-cm line opacity.In this study, we consider that the temperature within each zone is uniform at any instant.Thus, this simplification does not ignore the temperature evolution of the gases in each of the zones.
For both the two-zone and three-zone models, all zones expand only radially.The thickness of the partially ionised zone also increases radially in the three-zone model.To avoid using excessive number of parameters, which may cause unnecessary complications that can mask the physics and the 21-cm line structure, we consider that the partially ionised zone is also uniform at any instant.A schematic illustration of the two-zone and the three-zone models is shown in Fig. 1.

Ionisation induced expansion
We first consider the situation that the time for light traversing the H II zone is negligible compared to the Hubble time.In this situation, cosmological effects can be ignored, and the radiative transfer of the 21-cm line is subject only to the dynamical evolution of the H II zone, which is defined by the propagation of the ionisation front.Suppose that the H II zone has a radius  0 initially.It expands radially at a speed v.The radius of the H II zone at time  is then  () =  0 +   () where  = v/,  () = , and  is the speed of light.Now consider a plane light front propagating in the direction (in the Cartesian coordinates) from an initial reference location,  0 .We may set  0 = − 0 .Then the light front will reach () = − 0 +  () at time .It follows that, on the -plane, the light front will reach the H II zone boundary at the location (, ), given by Clearly, at  = 0, (, ) = (− 0 , 0).Also, we have implying that the path that the light front traverses through the H II zone is (3) The time that the light front meets the H II zone boundary for a given  therefore has two values, and they are For a non-expanding H II zone ( = 0), This means that expansion leads to LoS length distortion, in terms of a length ratio R: Note that Δ  = 0 occurs at instead of at  = ±  0 , the values expected for a non-expanding H II zone -it has grown bigger since the light front first intersects the H II zone boundary.Note also that the light path Δ and hence the light traversing time Δ diverge when  → 1 (see Equation 3).This corresponds to the situation where the H II zone is expanding so rapidly, in the speed of light, that light propagating from the boundary of one side of the H II zone will never reach the boundary of the other side of the H II zone, as seen by a distant observer.
In the two-zone model, the opacity of the H II zone is provided by the free-electron through Compton scattering.The scattering optical depth of radiation transported across the H II zone at  is where  e is the electron number density and  Th is the Thompson cross section.This can be generalised for ionised cavities with any arbitrary opacity  ν , at frequency ν, which gives a specific optical depth which is in contrast to the specific optical depth expected for a non-expanding cavity: The remaining question is now how fast the H II zone would expand.We may obtain a rough estimate for the expansion speed v from the following consideration.The boundary of the H II zone advances only when the amount of ionising radiation is sufficient to convert H I atoms to ions.That is, where   is the number density of ionising photons, and  HI is the number density of H I in the surrounding.The efficiency of converting H I atoms into ions may parameterised by a variable Υ, which is determined by the local atomic and thermodynamic properties (i.e.not the global geometry of the system) and the ratio ( HI /  ).The expression in Eqn. ( 10) is essentially a measure of how saturated the ionisation process is.
The number density of ionising photon   () attenuates with  due to the dilution of radiation over distance and the consumption of the ionising photons to facilitate the expansion.Thus, the expansion of the H II zone even slows down as ionisation approaches saturation, where processes, such as recombination, counteract ionisation when the photon supply becomes insufficient.Nonetheless, in the very initial stage when there is a sudden burst of ionising sources, the supply of   is abundant, and hence the ionisation is far from saturation.For   ≫ Υ HI , the expansion velocity v would be expected to almost reach the speed of light, i.e.  ≲ 1.For modest unsaturated ionisation, which gives rise a (constant) normalised expansion velocity  = 0.5, gives | * | = √ 3  0 , and  ν = 4  ν  0 = 2 ′ ν at  = 0 (i.e.radiative transfer along the symmetry axis of the H II zone).Increasing the expansion velocity to  = 0.75 gives | * | = √ 7  0 , and Here, it shows that ignoring the expansion of H II zone could lead to incorrect inferences for its size and the optical depth.
The situation is slightly more complicated, if there is a partially ionised zone enveloping the H II zone.The partially ionised zone determines the radiative transfer process together with the H II zone.The temperature in the partially ionised zone could be an important variable, especially in the context of radiative transfer of the 21-cm line, as its relative contrast with the temperatures of the background H I gas and the temperature of the H II gas would determine whether a line would appear as emission or as absorption relative to the continuum radiation.Apart from the microscopic aspects, such as radiation processes and radiative transfer, there are subtle issues in the macropscopic perspectives.The inner and outer boundary of the partially ionised zone could expand asynchronously at different speeds.There are also additional issues, which is generally insignificant for the development of the H II zone but might need to be take into account for the development of the partially ionised zone.This is due to the fact that the boundary of the H II zone is an ionisation front, not a hydrodynamics shock front, which needs to satisfy certain shock-jump condition, while the interface between the partially ionised zone and the H I zone could be determined by hydrodynamics and also thermodynamic processes.Nonetheless, we leave the investigation of these subtle yet important physics in a future paper.In this section, we employ a simple parametric model to illustrate how the presence of an additional partially ionised zone would alter the radiative transfer process.
Assume that the H II zone and the partially ionised zone expand in uniform normalised speeds: the outer boundary of the partially ionised zone has an initial radius  a0 and expand with  a .The boundary of the H II zone has an initial radius  b0 and expands with and  b (with  b ≤  a , which avoids the unphysical situation where the expansion of the H II zone overruns the expansion of the partially ionised zone).At  = 0, when the light front, propagating in the direction (as in the two-zone model) reaches the outer boundary of the partially ionised zone, which is located at  = − a0 .The results obtained above for the two-zone model can be adopted using the substitutions  0 →  a0 and  →  a .The advancing of the light front is now () = − a0 +  (), and the radius of the outer boundary of the partially ionised zone is  () =  a0 +  a  ().Thus, The corresponding maximum size of the partially ionised zone perceived by the distant observer is For the inner boundary of the partially ionised zone, some additional transformations for the corresponding variables are required, before applying the results obtained for the two-zone model.First, the light front reaches this boundary at a time t, at a location between  = − a0 and  = − b0 .Denote this location as  = − rb0 .We may then determine rb0 , t, and hence ξ (= t) by setting − a0 + ξ = − b0 −  b ξ.This gives and where  =  b0 / a0 ≤ 1.With the substitutions of  a0 → rb0 ,   →   and  → ( − ξ) in Eqns.( 11), ( 12) and ( 13), we obtain for the interception of the light front with the boundary of the H II zone.It follows that implying that where and is always positive.Also, In where  νa is the specific opacity of the gas in the partially ionised zone.We denote the corresponding specific optical depth for these three passages as ν and [k] ν respectively, and the total specific optical depth is the sum of them, i.e.  ν = Using similar procedures as in the derivation of Eqn. ( 21), we obtain for the H II zone1 .To derive the specific optical depths of the two passages through the partially ionised zone, we may consider the expression to simplify the algebraic steps.The specific optical depth corresponding to the first passage is the specific optical depth corresponding to the second passage is With non-zero  νa and  νb , summing the specific optical depths in Eqns.( 22), ( 24) and ( 25) yields The structures of real ionised cavities are expected to be more complex than that of the two-zone and three-zone model bubbles.Ionised cavities produced by quasars and galaxies have been studied in detail with analytical (e.g.Shapiro & Giroux 1987;Wyithe & Loeb 2004), semi-numerical simulations (e.g.Geil & Wyithe 2008;Mesinger et al. 2011) and hydrodynamical simulations (e.g.Thomas et al. 2009;Geil et al. 2017;Hutter et al. 2021;Kannan et al. 2022).While these studies focused on the global progression of reionisation, some also targeted individual bubbles (e.g.Ghara et al. 2017).The statistics properties of the cavities were also investigated (e.g.Zahn et al. 2007;Shin et al. 2008;Lin et al. 2016;Muñoz et al. 2022;Schaeffer et al. 2023;Lu et al. 2024).The 21-cm signal associated with complex ionisation structures are not ideal for comparing our covariant formulation with the optical depth parametersation.Hence we adopt the simplified the two-zone and three-zone models.

Cosmological expansion
Cosmological expansion affects the observational properties different to the local expansions, such as the geometrical expansion caused by the advance of an ionisation front.When the Universe expands, an ionised cavity expands accordingly, even when the ionisation front is stationary in the local reference frame.Cosmological expansion also alters the thermodynamics and hydrodynamic properties of the ionised cavities.It sets a new balance between the radiative processes, hence modifying the radiation and the observational characteristics of the ionised cavity.
Two effects are the most noticeable among the others.(i) The size of the cavity perceived by an observer at a lower redshift is larger than the size of the cavity measured in its local reference frame, where the radiative and hydrodynamic processes operates.(ii) The wavelength of the radiation is stretched as it propagates across the cavity.A 21-cm line originated from the far side of the cavity will have a wavelength longer than 21 cm when it reaches another side of the cavity.This change in the waveleghth will alter the transport of the radiation.For the 21-cm line, the interaction with its neighbouring continuum becomes prominent, while the resonance absorption can become insignificant.This has not taken account of  When the radiation enters the first H I slab, its specific intensity is  ν,0 ; when the radiation exits, it is  ν,1 .When the radiation enters the second H I slab, its specific intensity is  ′ ν,1 ; when the radiation exits, it is  ν,2 .The sequence of redshifts at which the radiation passes the slab boundaries in its propagation is further complications by effects associated with thermal and hydrodynamics evolution of the cavity and of the medium surrounding the cavity.Fig. 2 shows an illustration to elaborate these effects.A photon of a wavelength  1 starts its journey from the far side of a cavity boundary, at redshift  1 (corresponding to a cosmological time  1 ) and arrives at the near side of the cavity boundary, at a lower redshift  2 (corresponding to the a cosmological time  2 ).When the photon reached the other side of the cavity boundary its wavelength became  2 .The cavity was bathed in the CMB of temperature  CMB ( 1 ) at the moment the photon started its journey, but the temperature of the ambient CMB had dropped to a lower temperature  CMB ( 2 ) by the time the photon arrived at the other side of the cavity.
In a Friedmann-Lemaître-Robertson-Walker (FLRW) universe with a zero curvature, the expansion of the universe can be parameterised by a scale parameter (), which evolves with the cosmological redshift as This gives the stretch of the wavelength of radiation and the evolution of the CMB temperature The number density of particles evolves with the cosmological redshift as Consider that a photon propagates from a location  1 to a location at  2 .It starts from  1 at time  1 and reaches  2 at time  2 .As photons travel along null geodesics, in a flat FLRW universe  d  = () d along the LoS.It follows that Here  1 and  2 are the redshifts, with respect to a present observer located at  = 0 at redshift  = 0, respectively, and correspond to the time  1 (when the photon starts its journey from  1 ) and  2 (the time when the photon completes its journey reaching  2 ), with () = 1 at  = 0. Hence,  2 satisfies In the ΛCDM framework, (see e.g.Peacock 1999) where  0 is the Hubble parameter at present, and Ω r,0 , Ω m,0 and Ω Λ,0 are the density parameters for radiation, matter and dark energy, respectively.For the spherical H II zones in the two-zone model without expansion induced by ionisation (Fig. 1, see also Sec. 2.2), ( 1 −  2 ) = −2 0 along the symmetry axis.With this specified,  2 is readily determined for a given  1 , and this can be generalised for the two-zone and threezone models.Assigning the location of the zone boundaries as  2 , which specifies the value for ( 1 −  2 ).The corresponding  2 (and hence  2 ) can be obtained by a direction integration of Eqn.(32).
During the Dark Ages and the EoR, the cosmological evolution is matter dominated.We may ignore Ω r,0 and Ω Λ,0 , which results in an analytic expression for  2 : We now use a three-slab model to illustrate the difference between the spectra, consisting of a line component and a continuum component, in the presence and in the absence of cosmological expansion.
Consider a ray of radiation with initial specific intensity  ν,0 , passing through a slab of H II gas sandwiched by two thermal absorptive slabs of H I gas, both of co-moving thickness  (see Figure 3).Without losing generality, the total optical depth of the H II slab is assumed to be negligible, i.e. the slab practically does not contribute to absorption and emission of the radiation.The two H I slabs are transparent to the continuum but have sufficient optical depths at the line centre frequency, ν 21cm , and its neighbouring frequencies.The local absorption in the H I slabs is specified by an absorption coefficient, which takes the form: , where  a is the number density of the absorbers,  a is the normalised absorption cross section, ν 21cm is the centre frequency measured in the local rest frame, and A thermal absorptive media also emits and the local emissivity can be derived from the absorption coefficient and the source function, which is the Planck function  ν () in a thermal medium.(Note that, for simplicity, we have omitted the contribution from stimulated emission in this illustration.For the radiative transfer calculations of the hyperfine 21-cm H I line used for diagnosing cosmological reionisation, stimulated emission must be included.)At low-frequencies (i.e. in the Rayleigh-Jeans regime), the Planck function  ν →  RJ ν = ν 2 , where  = 2 B  −2 and  B is the Boltzmann constant.For Δ 1 ≪ (1 +  1 ) and Δ 2 ≪ (1 +  2 ), we may approximate that the emission and absorption are constant with time.Suppose also that the background is a thermal black-body of a temperature  0 .Then, in this situation, where  a,1 is the number density of absorber in H I slab 1 and  1 is the thermal temperature associated with the line-formation process in H I slab 1.As the H II slab does not contribute to absorption and emission, ( ν /ν 3 ) is invariant.Thus, where where  a,2 is the number density of absorber in H I slab 2 and  2 is the thermal temperature associated with the line-formation process in H I slab 2. Generally,  0 ≠  1 ≠  2 and  a,1 ≠  a,2 .The spectrum of the radiation that emerges has a continuum component and two line components, one associated the absorption in H I slab 1 and another associated with the absorption in H I slab 2. The continuum has a specific intensity which is independent of the thermal conditions in the two H I slabs.This is consistent with our assumption that the two H I slabs are optically thin to the continuum and the H II slab has no contribution to the radiation, and that ν −3  ν is an invariant quantity.For a sufficiently large difference between  1 and  2 such that the two line-profile functions do not overlap, the line associated with the first H I slab is centred at ν * 21cm = ν 21cm (1 +  2 )/(1 +  1 ), and the specific intensity at the line centre is The line associated with the second H I slab is centred at ν 21cm , and the specific intensity at the line centre is Here, it shows the two line components are not identical with the correction of the shift in the line centre energies, even for identical thermal conditions in the two H I slabs at time  0 , when the radiation enters the first H I slab.Even when the micro-astrophysical processes, such as heating and photon-pumping are absent, cosmological expansion would alter the temperature and the number density of the absorbers in the slabs.If cosmological expansion is insignificant over the interval when the radiation completes its journey through the two H I slabs and the H II slab, the spectrum has a continuum and only one line.The specific intensity of the continuum is simply The line is centred at frequency ν 21cm , and the specific intensity at the line centre is The line-centre specific intensity becomes the same as the specific intensity of the continuum neighbouring to the line, i.e. the line vanishes, when  2 =  1 =  0 .An emission line will result for  0 >  1 >  2 and an absorption line will result for  2 >  1 >  0 .These results are as expected from the line-formation criteria.
The effect of cosmological expansion in the transport of ionising photons and the development of ionised cavities have been studied analytically and implemented in reionisation simulations (see e.g.Shapiro & Giroux 1987;Bisbas et al. 2015;Weber et al. 2013;Fedchenko & Krasnobaev 2018).The subsequent effects on 21-cm signals are also discussed (Yu 2005).This cosmological effect is sometimes named as 'light cone effect' and may lead to anisoptropy in the 21-cm power spectrum (Barkana & Loeb 2006;Datta et al. 2012;La Plante et al. 2014).However, this effect is generally ignored in recent studies of 21-cm signals which are based on reionisation simulations.

COMPUTATIONAL SET-UP FOR 21-CM RADIATIVE TRANSFER CALCULATION
We briefly recapitulate the C21LRT equation and the input default reionisation history model used throughout this paper first.The computational set-up for the radiative transfer calculations as the same as in Chan et al. (2023), except for the inserted H II zones.The C21LRT equation in covariant form with our adopted FLRW cosmological model, when there is negligible scattering, is as follows, d d where the subscript "L" and "C" denote the 21-cm line and its neighbouring continuum (i.e.CMB in our calculations).Following their definitions in the previous section,  L,ν ,  L,ν and  L,ν are the specific intensity, line absorption coefficient and line emission coefficient at 21-cm line centre, respectively.Ξ is the factor for the stimulated emission. C,ν and  C,ν are the continuum absorption coefficient and continuum emission coefficient and  ν is the normalised line-profile functions. represents the photon's path length, and d/d =  d/d can be calculated based on our assumed cosmology (see Eqn. ( 33)).
The line coefficients are calculated in a local rest frame (the same practice as in Fuerst & Wu (2004); Younsi et al. (2012); Chan et al. ( 2019)) as follows, where the subscript "u" denotes the upper energy state and "l" denotes the lower energy state of the H I hyperfine transition.Here, ν ul = ν 21cm ,  u and  l are the multiplicities (degeneracies) of the upper and lower energy states, respectively,  u and  l are the number density of particles in the upper and lower energy states, respectively, and  ul and  ul are the Einstein coefficients.We assume that the normalised line profile functions  ν,x ≡  x (ν − ν line,0 ), with x ∈ {abs, emi, sti} corresponding to absorption, spontaneous emission and stimulated emission, respectively.Then the factor for the stimulated emission is Ξ = ( u  l )/( l  u ).
With this C21LRT equation, the input H I gas properties in EoR are the globally averaged spin temperature  s () and ionised fraction  i () from the 'EOS 2021 all galaxies simulation result' (which used the 21CMFAST code) (Muñoz et al. 2022).The density of H I gas is calculated based on  i and the adopted cosmological model.Then  u and  l are calculated with  s .We only consider CMB and leave other types of continuum emission to future studies.One 21-cm line profile is adopted for all redshifts throughout one calculation and it is parameterised as follows, where is the Doppler parameter.
As little do we know about the turbulent velocity of HI gas, from observations or from reliable theoretical modelling, in the ionisation front created by luminous sources, we choose the turbulent velocity of gas in intra-cluster medium (ICM), IGM and inter stellar medium (ISM) as benchmarks.The characteristic v turb for ICM is ∼ 100 km s −1 , created by large scale shocks (Subramanian et al. 2006 et al. 2018) and can the turbulence can be driven by cold gas infall, gravitational instability or star formation activities (Patrício et al. 2018).We therefore adopt v turb ∼ 1, 10, 100 km s −1 in our calculations.As will be demonstrated, the exact value of v turb does not effect our conclusions.We then assume that  k is always 0 in this paper 2 .More detailed description of the C21LRT formulation and the adopted default reionisation history can be found in Chan et al. (2023).
With the default reionisation history specified, we then insert H II zones of spherical shapes (referred as 'bubbles' in the following texts).Various methods have been employed in the literature to extract the topology and size distributions of ionised cavities from large scale simulation studies (e.g.Furlanetto et al. 2006a;Shin et al. 2008;Friedrich et al. 2011;Lin et al. 2016).These studies predict the distribution of sizes of ionised cavities and study the evolution of the distribution function throughout EoR.The characteristic sizes varies from ≲ 1 Mpc at the beginning of EoR to ≳ 10 Mpc at the end of EoR (Iliev et al. 2006;Shin et al. 2008;Friedrich et al. 2011;Lin et al. 2016).To cover the possible sizes, we choose bubble diameters  of 0.01, 0.1, 1.0, 10.0 Mpc.Only one bubble is inserted in each scenario.For most of the scenarios, radiative transfer calculation is carried out along a single ray through the bubble centre.The comoving distance intercepted by the ray is .We always fix the boundary of the bubble at higher redshift e.g. begin = 12.0 first.With the size d of the bubble specified by , we then calculate the boundary of the bubble at lower redshift  end =  begin − d.Along this ray,  HI = 0 in the redshift cells between  begin and  end .We insert bubbles in this way to facilitate the comparison of bubble features.We do not modify  s in these redshift cells (as this should cause no difference).Each bubble in our calculation is resolved with at least 10 redshift cells.
2 To match v turb = 10 kms −1 ,  k = 6060.67K which is already higher than the expected globally averaged temperature of H I (cf. Figure 1 of Chan et al. (2023)).Also, such a high temperature is usually accompanied with high ionisation fraction  i , diminishing the 21 cm signal from HI gas.end  begin = 12.0, which corresponds to ν = 109.2620MHz, is marked with red vertical line in all three panels.In the top panel, the globally averaged 21-cm spectrum without any bubble is plotted in thick grey line.The bubble features created by bubbles with diameters  = 0.01, 0.1, 1.0, 10.0 Mpc are also plotted with think black lines in contrast.The bubble features are not distinguishable from each other in this panel and are shown in more details in the middle and bottom panels. = 12.0 is inside the absorption regime and all the bubble features can be understood as the decrease of the absorption signal.In the middle panel, the features created by bubbles with  = 10.0, 1.0, 0.1 Mpc are plotted in black with solid, dashed and dash-dot lines.The original 21-cm spectra (in grey) and the features created by smaller bubbles (in blue) are also preserved in this panel for comparison.
The transition between a bump and a valley happened when the bubble is ≳ 10 Mpc, this transition bubble size is mostly determined by v turb as explained in the main text.In the bottom panel, the features created by bubbles with  = 0.01, 0.1 Mpc are plotted in blue.The original spectra is plotted in grey but not visible as the curve representing  = 0.01 Mpc is almost identical in this panel.

What determines the properties of bubble features
We can first estimate the frequency range where the 21-cm spectra is affected by a bubble specified with  begin and  end .Its effective size when measured in frequency (at  = 0) dν bubble is It is natural to assume that the spectrum will be affected in frequencies from ν 21cm (1 +  end ) −1 to ν 21cm (1 +  begin ) −1 .We also expect the 'neighbouring' frequencies to be affected due to the finite end  begin = 8.0, which corresponds to ν = 157.8229MHz, is marked with red vertical line in all three panels. begin = 8.0 is inside the emission regime and all the bubble features can be understood as decreasing of the emission signal.In the top panel, the globally averaged 21-cm spectrum without any bubble is plotted in thick grey line.The bubble features created by bubbles with diameters  = 0.01, 0.1, 1.0, 10.0 Mpc are also plotted with thin black lines in contrast.The bubbles features by larger bubbles and smaller bubbles are shown in more details in the middle and bottom panels.In the middle panel, the features created by bubbles with  = 10.0, 1.0, 0.1 Mpc are plotted in black with solid, dashed and dash-dot lines.The original 21cm spectra (in grey) and the features created by smaller bubbles (in blue) are also preserved in this panel for comparison.The transition between a dip and a valley happened when the bubble is ≳ 10 Mpc, this transition bubble size is mostly determined by v turb as explained in the main text.In the bottom panel, the features created by bubbles with  = 0.01, 0.1 Mpc are plotted in blue.The original spectra is plotted in grey but not visible as the curve representing  = 0.01 Mpc is almost identical in this panel.
width of 21-cm line.The full width half maximum of the 21-cm line profile when measured with frequency in the local rest frame is Due to cosmological expansion, a 21-cm line with FWHM ν at  will become narrower, i.e. ν  = FWHM ν 1 1+ when propagated to  = 0. We therefore expect the spectrum to be affected at frequencies ν 21cm (1 +  end ) −1 ± dν  and ν 21cm (1 +  begin ) −1 ± dν  .Note that we defined both dν bubble and dν  at  = 0 to facilitate the analysis of bubble features in the observed spectra, which are at  = 0.
This estimation is consistent with the C21LRT results.We first analysed scenarios at  = 12.0 calculated with v turb = 100 km s −1 as shown in Fig. 4. The original spectra without any bubble is plotted with thick grey line and the spectra with various bubble sizes in .Gaussian fit analysis of the features created by a  = 0.01 Mpc bubble at  = 12.0.In all three panels, the difference between the original spectra ( nobubble ) and the spectra when there is a  = 0.01 Mpc bubble at  = 12.0 ( bubble ) is plotted with black solid lines.The difference in spectra ( bubble −  nobubble -ν) is then fitted with a Gaussian function, shown in yellow dashed curves.From top to bottom panel, the bubble features were calculated with v turb = 100, 10 km s −1 , respectively.As explained in the main text, the widths of Gaussian features are dominated by the line profiles because a 0.01 Mpc bubble at  = 12.0 is narrower compared to the widths of the line profiles, when measured in frequency.The amplitudes of the Gaussian features can be approximated as  bubble −  nobubble ∼  •  L,ν 21cm , see main text for details.black in the top panel.In the middle panel, we zoomed into the relevant frequency range for the inserted bubbles.The smaller bubbles with  = 0.01, 0.1, 1 Mpc are plotted in blue.Their effective widths are dν bubble = 4.9572 × 10 −4 , 4.9574 × 10 −3 , 4.9670 × 10 −2 MHz, respectively, at  = 0.These bubbles are narrower than the width of 21-cm line profile (dν bubble ≪ dν  = 0.06069 MHz), and the spectral width of the bubble feature (at  = 0) dν feature is dominated by dν  .The effective size of the 10 Mpc bubble (5.0018 × 10 −1 MHz) is larger than dν  .The corresponding bubble feature is plotted in black, with width dominated by dν bubble .It is significantly broadened to be plateaus instead of bumps.The absolute value of specific intensity (| L −  C |) is reduced to 0 at the plateau.In the bottom panel, we zoom in further to show the shape of the features created by smaller bubbles.Comparing the three curves calculated with  = 0.01, 0.1 Mpc, we can see that whilst the widths of the bubble features are similar (dominated by dν  ), the heights of these features increase as the sizes of the bubble increase.
The analysis of bubbles in the emission regime of the redshifted 21-cm global signal is similar.These bubbles have  = 0.01 − 10 Mpc at  begin = 8.0.Their features are also calculated with v turb = 100 km s −1 as shown in Fig. 5. Bubbles manifest as dips or valleys (reduced emission) as shown in the top panel.We zoomed into the bubble region in the middle panel.The bubble features created by small bubbles ( = 0.01, 0.1, 1 Mpc, dν bubble = 5.9570 × 10 −4 , 5.9572 × 10 −3 , 5.9697 × 10 −2 MHz) have a dip like shape (in blue).The 10 Mpc bubble (5.9989 × 10 −1 MHz) has effective size larger than dν  = 0.8766 MHz and produce a valley-like shape.The transition between these two types of shapes is also determined by dν  .The features of smaller bubbles are enlarged in the bottom panel.Their depths also increase with the size of the bubble dν bubble .
In the results above, we chose a large turbulent velocity (v turb = 100 kms −1 ) to clearly demonstrate features created by the bubbles of  = 0.01 − 10 Mpc, also the transitions from dip (peak) to valley (plateau) when bubble sizes increases in emission (absorption) regime.As the transition is determined by comparing dν bubble to dν  , the critical bubble size where the transition happens is determined by dν  (which is determined by v turb in our calculation).To demonstrate this dependence, we also calculated the bubble features at  = 12 with v turb = 10 km s −1 , as shown in 6.The line profile width for this v turb is dν  = 0.006069 MHz.In the top panel, the transition now happens between  = 0.1 Mpc and  = 1 Mpc (the transition happens between  = 1 Mpc and  = 10 Mpc for v turb = 100 km s −1 in the middle panel of Fig. 4).

Analytical approximations for bubble features
For a  = 0.01 Mpc bubble at  = 12.0, its effective size dν bubble = 4.9572 × 10 −4 MHz is smaller than the line profile width dν  when adopting turbulent velocities v turb = 100, 10 km s −1 .The width of the bubble features calculated with these v turb are therefore all dominated by dν  .We then fit the features created with  = 0.01 Mpc with the line profile function (i.e.Gaussian functions in our calculations) for a more quantitative analysis.
The fitting results in Fig. 7 are for v turb = 100 and 10 km s −1 in the upper and lower panel.The bubble features calculated with C21LRT (black solid lines) and the fitted Gauasian functions (yellow dashed lines) are consistent over the frequency ranges from where the fitting was carried out.The frequency corresponding to the high redshift boundary  begin = 12.0 of all these bubbles are marked with red vertical lines.The fitted width dν fit and height  fit of the Gaussian functions are listed in Table 1.The fitted width dν fit for the two bubble features are very close to dν  as determined by v turb .Also,  fit is very close to the maximum specific intensity of the bubble feature  max .We calculated the relative difference between the bubble feature and the fitted Gaussian function and averaged it over from (ν 21cm − 1.5FWHM ν )/(1 + ) to (ν 21cm + 1.5FWHM ν )/(1 + ) 3 .The averaged relative difference  for the two bubble features are all less than 10 −4 .The small values of dν  − dν fit ,  max −  fit and  shows that the bubble features are very close to Gaussian functions and the bubble width is well approximated with dν  .
We then analyse the height of the feature  fit of these bubble features.As discussed above,  fit grows with dν bubble .We also expect it to be proportional to the 21-cm emission coefficient 4 .We found that  fit can be approximated as where  redshift = (1 + ) 3 results from cosmological expansion (we may view this the factor as the ν −3 factor in the expression of the invariant intensity). dν = 1 +  is caused by the width of the relevant integration frequency range scales as 1+ when the features propagate from  to the observer on earth.The  convolution factor resulted from convolution effect of line broadening in the local rest frames and the propagation of radiation.The values of  analytical for the three bubble features are also listed in Table 1.By comparing  fit with  analytical , the value of  convolution is approximately 0.75.
These three examples are calculated with fixed dν bubble and different dν  .The results show that when the bubble feature is dominated by the line profile (dν bubble ≪ dν  ), the width of its spectral signature is well approximated by dν  .The height of the bubble's spectral feature is proportional to dν bubble /dν  and can be approximated with Eqn.51 .This analysis also applies to other small bubbles.For example, at a given redshift the bubble feature produced with  = 0.1 Mpc calculated with v turb = 100 km s −1 is comparable to that produced with  = 0.01 Mpc calculated with v turb = 10 km s −1 .
The large bubbles can be approximated in a similar way.The features created by large bubbles (which satisfy dν bubble ≫ dν  ) can be separated into three part.For example, the 10 Mpc bubble feature in the middle panel of Fig. 4 starts and ends approximately at frequencies and This features can be approximately separated by and into left, middle and right part.The left part (from ν begin to ν left ) and Table 1.Results of Gaussian fit analysis.We fitted the bubble feature with a Gaussian function (see main text for the details of bubble feature and the fitting processes).From left to right, the columns are, 1 turbulent velocity v turb adopted when calculating the bubble feature for  = 0.01 Mpc at  begin = 12.0, 2. fitted width at half maximum dν fit of the bubble feature, 3. adopted line profile width at half maximum dν  as determined by v turb at  begin = 12.0, propagated to  = 0, 4. difference between dν  and dν fit , 5. the maximum specific intensity (difference in intensity between the result with and without bubble)  max , 6. the fitted height of Gaussian  fit , 7. approximated height of the bubble feature calculated using Eqn. 51,8. relative  For bubbles of intermediate size (dν bubble ∼ dν  ), they produce spectral features between ν left and ν right .Their spectral features are similar to a Gaussian function but can not be fitted perfectly with a Gaussian function.This is because their spectral features are determined by the convolution between the line profile and other radiative transfer effects over a relatively wide redshift range 5 .
We note that these approximations are only valid because we adopted a constant v turb , hence a constant normalised line profile, for each calculation and all the bubbles we inserted are fully ionised.The features will be too complicated to be approximated analytically if we adopt realistic  HI and v turb values, which could have large spatial variation within one ionised cavity created by an astrophysical source.

Resolve a bubble with multiple rays
When ionised bubbles are large enough, it is possible to resolve each of them with multiple rays and study the changes in 21-cm features due to their spatial variations.
We used ten rays to resolve a bubble with  = 5.0 Mpc ( = 2.5 Mpc) at  = 8.Different to the previous calculations, here we set the bubble centre at  = 8 and calculated  begin and  end for each ray.The bubble and the ten rays are illustrated in Fig. 8.The first ray is determined by the bubble centre and the observer on earth (d = 0), which also determines the LoS.All the other 9 rays are parallel to the first ray, with distance d = /10, 2/10, ......, 9/10.We use d to denote distance to bubble centre in the perpendicular direction with respect to LoS.The edge of the bubble is marked by the thick grey line (d = ), where the ray does not intersect the H II bubble.We adopted a turbulent velocity of 100 km s −1 6 .
The 21-cm spectra at  = 0 for the bubble region are shown in Fig. 9.The spectra are coded with the same color and line styles as 5 When we fit these bubble features with Gaussian function, the fitted Gaussian function is flatter than the bubble feature, i.e. the height of bubble feature |  max | is noticeable larger than the fitted height |  fit |. 6 We note that this 100 km s −1 turbulent velocity may not be achieved for a realistic ionisation bubble.It is chosen such that the distortion of the apparent shape of the ionisation front, due to the finite speed of light, can be clearly demonstrated.the corresponding rays in Fig. 8.The bubble features grows in both depth and width from d = 9/10 to d = 0.Each bubble feature appears symmetrical to the frequency that corresponds to the bubble centre at  = 8.0.When d = 9/10, the distance intercepted by the ray is small, and the bubble feature is still dominated by line profile width.As d decreases, the intercepted distance grows and dominates over the line profile width.So far, we have been measuring the bubble size or the intercepted distance in the rest frame of the (assumed) ionising source at the centre of the bubbles.This is no longer accurate if the bubble is expanding with velocities v comparable to the speed of light.When v is comparable to , the size of the bubble changes significantly within the light crossing time (e.g. of the ray which intersect the bubble centre) as discussed in Sec.2.3.The apparent shape probed with 21-cm line is different from the shape observed from the local rest frame of the ionising source.We followed Yu (2005) and calculated the apparent shape of a bubble with a constant v = 0.9 , as shown with the black ellipse in Fig. 10.Instead of assigning rays by d = /10, 2/10, ..., we first calculate the maximum distance to the centre ray  max and assign rays with d =  max /10, 2 max /10, ...9 max /10.We can see that due to a large v, the apparent shape significantly changed7 .The intercepted distances decreased for all the ten rays compared to Fig. 8.The resulting 21-cm spectra  = 0 are shown in Fig. 11.The rays close to the edge of bubbles intercept much shorter distances (where there is no H I) and the width of the corresponding bubble features are dominated by line profile width.All the spectra features are narrower compared to those in Fig. 9.These spectra no longer appear symmetrical with respect to the frequency marked with red vertical line ( = 8.0).

Differences between C21LRT and optical depth parametrisation
Using identical HI gas properties as inputs, the absolute difference in the spectra at  = 0 calculated with C21LRT and optical depth parmetrisation method, denoted by , primarily arises from two factors.First, the optical depth method tracks the signal evolution only in redshift space, while C21LRT tracks 21-cm line in a two-dimensional (redshift and frequency) space.This results in an approximate 5 − 10 % difference across all redshifts.The second factor stems from the simplifications used in deriving the optical depth parametrisation, at the expense of neglecting smallscale variations of HI gas.Both factors are discussed in Chan et al. (2023).The second factor is evident when HI gas properties undergo abrupt changes, such as the transition between neutral and ionised zones.For a non-zero FWHM ν ,  reaches its maximum at approximately ν left and ν right as defined in Eqns.54 and 55, where  L −  C immediately drops to 0 based on optical depth method, whereas || ≈ ( L −  C )/2 based on C21LRT calculations 8 .For smaller bubbles where dν  > dν bubble , || remains greater than 0 throughout the relevant frequency range (from ν begin to ν end , as defined in Eqns.52 and 53).For larger bubbles, || decreases to 0 between ν left and ν right .These differences mostly result from the rest frame 21-cm broadening, with a minor contribution from radiative transfer effects.
The bubble features can also be understood as follows: For the optical depth method (with v turb = 0 implicitly assumed), the 'loss' of the 21 cm signal due the presence of an ionised bubble is distributed between ν 21cm /(1+ )| = begin and ν 21cm /(1+ )| = end .For the cases presented in Figs. 4, 5, 6, 7, the value of v turb determines the FWHM ν (eq.50).The 'loss' of 21 cm signal is distributed over a broader frequency range between ν begin and ν end .The shape of the spectral feature is determined mainly by the 21-cm line profile.
Additionally, when using the optical depth method to compute the 21-cm signals seen by a present-day observer, it often involves a separate prescription to incorporate the light from past epochs that can reach us now.This is often done by constructing a light cone from sequential snapshots of the simulated box, where the same part of the universe has evolved with cosmological time.However, the optical depth method is inaccurate when 21-cm line is broadened or not optically thin (Chapman & Santos 2019).It can also overlook certain cosmological effects, such as the alterations in the apparent shape of the expanding ionisation fronts.

Implications of our results
For more sophisticated ionised cavities in cosmological simulations, their impacts on the 21-cm spectra and power spectra are currently all calculated with the optical depth parameterisation.This means that the radiative processes (e.g.line broadening due to turbulent velocity) which happens on length scales of or smaller than the sizes of the ionised cavities are missing in their calculation.When observational data become available, the sizes of ionised cavities in EoR inferred from their calculated 21-cm power spectra can have large uncertainties.For example, if ionised cavities all have sharp boundaries (the ionised fraction changes quickly from 1 to 0) and have a characteristic diameter , the optical depth parametrisation would predict a sharp change in the power spectra near  = 1/, while full radiative transfer calculations would predict gradual changes near  = 1/.

Some remarks on the ionised cavity models
In this work we have used the two-zone model to show the features that it imprints on the 21-cm spectra.The two-zone model is a plausible representation of the ionised cavities carved out by the very luminous astrophysical sources when the partially ionised zone has negligible size, such as quasars 9 .The propagation of the ionising front created by powerful ionising sources, as captured in the expanding bubble scenario (presented in Sec 4.1.3)has demonstrated that the time evolution of ionisation structures needs to be corrected.
In some situations, the two-zone model may not be able to fully captures thermodynamics and geometrical properties of the ionised cavities.For instance, UV radiation from individual sources such as Population III stars or Population II stars would ionise the surrounding H I gas gradually over time.This would allow the photoionisation and recombination processes to reach a quasi-equilibrium state, 9 We tested the three-zone model by adding a thin partially zone near the boundary (similar to the illustration in Figure 1) and found that the features in the 21-cm spectra of the two-zone and three-zone models are actually very similar.We therefore do not show the results of three-zone models in this paper.
in which the size of the cavity would be approximately the Strömgren radius (see Strömgren 1939).In the more complicated situations, such as that the ionisation source is a group of star-forming galaxies or a galaxy hosting AGNs, the radiation flux and spectrum would vary with time.The interplay between microscopic processes, such as photoionisation, recombination and radiative cooling loss, and macroscopic processes, such as hydrodynamical bulk flows, will give rise to instabilities.The distribution of H I gas near the ionising sources is no longer uniform and smooth.The media both inside and outside the ionised cavities can be clumpy, and the 21cm spectral signature of these cavities will be very different to those obtained from the model bubbles constructed using a globally averaged  HI .The development of the ionisation structure created by these sources are jointly determined by the ionisation, recombination processes as well as the thermodynamic processes (Axford 1961;Newman & Axford 1968;Yorke 1986;Franco et al. 2000).Transition layer of partially ionised gas can vary in thickness, ionisation state and thermodynamics properties.Furthermore, in the presence of X-ray emissions, which penetrate more deeply than UV photons and are more efficient at heating gas, ionised cavities are expected to much larger with a hazier boundary.Over time, ionised bubbles also came to overlap.All these factors are expected to individually and collectively affect the 21-cm power spectra (Finlator et al. 2012;Kaurov & Gnedin 2014;Hassan et al. 2016;Mao et al. 2020).Our future studies will look into these complex ionisation structures in more detail, as C21LRT can interface with the simulated data with small scale (≲ 10 kpc) structures in the future (see Appendix A for details).

CONCLUSION
We investigated the 21-cm spectral features imprinted by individual spherical ionised cavities enveloped by H I medium, adopting a covariant formulation for tracking 21-cm signals in both redshift and frequency space.We studied the improvement in accuracies of the 21-cm signals, compared to adpoting the optical depth parametersation, which tracks 21-cm in the redshift space only.
We use a cosmological radiative transfer code C21LRT for the numerically calculation of 21-cm signals.We showed that the evolution of ionised cavities as driven by ionisation, thermodynamical and cosmological effects, would affect the apparent shape of the cavities when probed with 21-cm line.The apparent shape of an evolving cavity is different from its shape in the rest frame of its stationary cavity centre, as there is a time lap between the radiation from the farside and from the nearside of the cavity.
We employed a set of single-ray calculations to show how the spectral features evolve with various bubble diameters  = 0.01, 0.1, 1.0, 10.0 Mpc.We found that the widths of spectral features imprinted by bubbles are jointly determined by the 21-cm FWHM ν , the bubble diameter  and the redshift of the bubble .The spectral feature width at redshift  can be approximated with max(dν  , dν bubble ), where dν  =FWHM ν /(1 + ) and dν bubble is the bubble diameter when measured in frequency space at  = 0.
When the 21-cm line FWHM ν dominates (dν bubble ≪ dν  ), as we adopted a Gaussian function shape 21-cm line profile, the bubble spectral feature has a Gaussian function like shape.For bubbles with smallest  in our paper (0.01 Mpc), their spectral features can be fitted with Gaussian function with negligible residuals.As  increases (dν bubble ∼ dν  ), the convolution between line profile and other radiative transfer and cosmological effects becomes nonnegligible and their spectral shape deviates from Gaussian func-tions.For these bubbles (dν bubble ≲ dν  ), the 21-cm specific intensity  L −  C does not decrease to 0 in the corresponding frequency ranges.
For large bubbles (dν bubble ≫ dν  ), the widths of their spectral features are dominated by dν bubble .These features can be divided into three parts.The shape of the left and right parts of the spectral feature is similar to half Gaussian functions.They are produced by the transition between H I and H II zone and their shape is mainly determined by the 21-cm line profile.The middle part of their spectral feature has  L −  C = 0.In comparison, the optical depth parametersation predicts that 21-cm signal will diminish to 0 for a fully ionised bubble, regardless of its size.We showed that the 21-cm intensities computed with it can have large discrepancies in the transition zones (ionisation fronts) from those computed with C21LRT.
The 21-cm signals associated with length scales equal to or smaller than the sizes of the ionized cavities need to be tracked in both redshift and frequency space for accuracy.For these scales, physical processes, such as line-broadening due to turbulent motion of the gas imprints on the 21-cm signals, can not be accurately computed with the optical depth parametrisation.Explicit covariant radiative transfer, such as the C21LRT, is necessary for correctly and self-consistently accounting for the convolution of local (thermodynamics and atomic processes and bubble dynamics) and global (cosmological expansion) effects onto the radiation that we receive from EoR.

Figure 1 .
Figure 1.A schematic illustration of the two-zone and the three-zone models.The vertical lines are the locations of the light fronts.In the two-zone model (top row), they are are marked by  (−) and  (+) , at the moments when the front reaches and leaves the H II zone respectively (at  = 0).In the three-zone model (middle and bottom rows), they are marked by  a(−) and  b(−) ,  b(+) and  a(+) , at the moments when the front reaches the outer boundary of the partially ionised zone, finishes the first passage through the partially ionised zone, starts the second passage through the partially ionised zone, and leaves the partially ionised zone, respectively (at  = 0).
the context of radiative transfer (and ray-tracing), a fraction of the radiation (those at | * a | ≥ || > | * b |) will pass through the partially ionised zone, and another fraction of the radiation (those at | * b | ≥ ||) will pass through the partially ionised zone before entering the H II zone and then pass through the partially ionised zone later afterwards.For the rays at | * a | ≥ || > | * b |, the specific optical depth, at frequency ν, across the H II and partially ionised zones is For rays with | * b | ≥ ||, the specific optical depth is the sum of three components: two from the partially ionised zone and one form the H II zone.Along the ray, they are the segments of the light path specified by ( a(−) ,  b(−) ) for the first passage in the partially ionised zone, ( b(−) ,  b(+) ) for the passage in the H II zone, and ( b(+) ,  a(+) ) for the second passage in the partially ionised zone.
) for | * b | > ||.Note that the expression for the specific optical depth above becomes the same as that in Eqn.(21), by equating  νb and  νa , regardless of what values  b takes.However, if we set  → 1 and  b →  a in Eqn.(26), then Eqn.(21) becomes a special case of it where  νb equals  νa .

Figure 2 .
Figure2.An illustration showing the changes in the properties of the radiation, e.g. its wavelength , and its environment, e.g. the CMB temperature  CMB , along the LoS between the moment that it started propagating from the far side of the cavity boundary (with variables denoted by the subscript "1") and the moment that it reached the near side of the cavity boundary (with variables denoted by the subscript "2"), as perceived by a distant observer at present.A useful reference for the changes is the cosmological redshift , which has a 1-1 correspondence to a cosmological time .The rate of change in the scale factor  in the FLRW metric is a measure of the cosmological expansion, which may be parameterised by  or .The same LoS is drawn twice in the figure to illustrate separately the conditions when the radiation started and finished its propagation across the cavity.

Figure 3 .
Figure 3.A schematic illustration of the three-slab model to show radiative transfer effects on line and continuum of radiation through an H II slab sandwiched by two initially identical thin H I slabs in an expanding universe.When the radiation enters the first H I slab, its specific intensity is  ν,0 ; when the radiation exits, it is  ν,1 .When the radiation enters the second H I slab, its specific intensity is  ′ ν,1 ; when the radiation exits, it is  ν,2 .The sequence of redshifts at which the radiation passes the slab boundaries in its propagation is {  1 ,  1 − Δ 1 ,  2 + Δ 2 ,  2 }.The corresponding sequence of time is { 0 ,  1 ,  1 ′ ,  2 }.

Figure 4 .
Figure 4. Features in the 21-cm spectra at  = 0 created by bubbles at  begin = 12.All the spectra in this figure are calculated with v turb = 100 km s −1 .The boundary of all the bubbles at the high redshift

Figure 5 .
Figure 5. Features in the 21-cm spectra at  = 0 created by bles at  begin = 8.0.All the spectra in this figure are calculated with v turb = 100 km s −1 .The boundary of all the bubbles at the high redshift

Figure 6 .
Figure 6.Features in the 21-cm spectra at  = 0 created by bubbles at  begin = 12.All the spectra in this figure are calculated with v turb = 10 km s −1 .In the top panel, the globally averaged 21-cm spectrum without any bubble is plotted in thick grey line (identical with some of the blues curved hence not visible).The bubble features created by bubbles with diameters  = 0.01, 0.1, 1.0 Mpc are plotted in this panel.Due to a smaller v turb , the transition between a bump and a valley happened when the bubble is ≳ 1 Mpc.In the bottom panel, the features created by bubbles with  = 0.01, 0.1 Mpc are plotted in blue.The original spectra is plotted in grey, although only visible in the bottom panel.The frequency ν = 109.2620MHz corresponding to  begin = 12.0 is marked with red vertical lines in both panels.

Figure 8 .Figure 9 .
Figure 8.Using ten rays to resolve a bubble with  = 2.5 Mpc at  = 8.The first ray is determined by the bubble centre and the observer on earth (d = 0), plotted in black.All the other 9 rays are parallel to the first ray, with distance d = /10, 2/10, ......, 9/10.We use d to denote distance to the first ray.The edge of the bubble is marked by the thick grey line (d = ), where the ray does not intersect the H II zone.

Figure 10 .Figure 11 .
Figure10.Using ten rays to resolve an expanding bubble at  = 8 whose ionisation front is expanding with a velocity of v = 0.9  in the local frame of an assumed ionising source, marked by the red cross.When the bubble has expanded to  = 2.5 Mpc in the rest frame (marked by the blue circle), we calculate the apparent shape due to finite speed of light (see main text for details).The apparent shape is indicated with the black ellipse.The ten rays are equally spaced, with the first ray intercept the bubble centre (d = 0), plotted in black.The other 9 rays have distance d =  max /10, 2 max /10, ......, 9 max /10.The edge of the bubble is marked by the thick grey line (d = ), where the ray does not intersect the H II zone.
difference between fitted Gaussian function and the calculated bubble feature averaged over from (ν 21cm − 1.5FWHM ν )/(1 + ) to (ν 21cm + 1.5FWHM ν )/(1 + ).right part (from ν right to ν end ) are both approximately half of a Gaussian shape and their width are approximately 0.5FWHM ν /(1 + ) = 0.5dν  , the intensity of the middle part (from ν left to ν right ) is very close to 0. These estimated frequencies are of ≲ 0.1 MHz accuracy.For example, ν left = 109.2923MHzforthe 10 Mpc bubble feature in Fig.4; the bubble feature ( bubble −  nobubble ) computed from C21LRT calculation has the (left) maximum value at 109.3358 MHz.