The morphology of the redshifted 21-cm signal from the Cosmic Dawn

The spatial fluctuations in the tomographic maps of the redshifted 21-cm signal from the Cosmic Dawn (CD) crucially depend on the size and distribution of the regions with gas temperatures larger than the radio background temperature. In this article, we study the morphological characteristics of such emission regions and their absorption counterparts using the shape diagnostic tool SURFGEN2. Using simulated CD brightness temperature cubes of the 21-cm signal, we find that the emission regions percolate at stages with the filling factor of the emission regions $FF_{\rm emi}\gtrsim 0.15$. Percolation of the absorption regions occurs for $FF_{\rm abs}\gtrsim 0.05$. The largest emission and absorption regions are topologically complex and highly filamentary for most parts of the CD. The number density of these regions as a function of the volume shows the power-law nature with the power-law indexes $\approx -2$ and $-1.6$ for the emission and absorption regions, respectively. Overall, the planarity, filamentarity and genus increase with the increase of the volume of both emission and absorption regions.


INTRODUCTION
Our Universe became neutral during the epoch of recombination about four hundred thousand years after the Big Bang.At a later stage, radiation emitted from the first stars, galaxies, Quasars, Highmass X-ray binaries (HMXBs), etc. heat up and ionize the cold and neutral gas in the inter-galactic medium (IGM).Theoretical studies such as Pritchard & Furlanetto (2007); Mesinger et al. (2011); Ghara et al. (2016); Islam et al. (2019); Ross et al. (2019) suggest that the first X-ray sources such as HMXBs, and mini-Quasars heat up the gas in the IGM during a period earlier than the stage when the IGM became highly ionized.This period is known as the Cosmic Dawn (CD) of our Universe.The CD was succeeded by a period, known as the Epoch of Reionization (EoR), in which the IGM's primordial neutral gas became ionized by the UV radiation from the sources.The physical processes of these eras have major roles in shaping the current state of our Universes through various feedback mechanisms.Unfortunately, these epochs remain among the least understood epochs of Cosmic history (e.g.Fan et al. 2006;Bañados et al. 2018;Planck Collaboration et al. 2020;Mitra et al. 2015).
★ E-mail: ghara.raghunath@gmail.comMany details about the timing, feedback mechanisms, morphology of the ionized and heated regions, etc. are still poorly known.
Observation of the redshifted 21-cm line radiation emitted from the IGM's neutral hydrogen (Hi ) atoms is the most promising probe of the physical processes during the CD and EoR (see e.g., Madau et al. 1997;Furlanetto et al. 2006;Zaroubi 2013;Murmu et al. 2022;Shaw et al. 2023a;Bera et al. 2023).Many of the world's largest radio telescopes have dedicated their valuable resources to measuring the Hi 21-cm signal from these periods.One type of such observations uses radio interferometers such as the Low-Frequency Array (LOFAR) 1 (Mertens et al. 2020), the Precision Array for Probing the Epoch of Reionization (PAPER) 2 (Kolopanis et al. 2019), the Murchison Widefield Array (MWA) 3 (e.g.Wayth et al. 2018) and the Hydrogen Epoch of Reionization Array (HERA) 4 (DeBoer et al. 2017), the New Extension in Nançay Upgrading LOFAR (NenuFAR) 5 (Munshi, S. et al. 2024), the Am-sterdam ASTRON Radio Transients Facility And Analysis Center (AARTFAAC) (Gehlot et al. 2022).Due to limited sensitivity, the ongoing radio interferometer-based 21-cm signal observations of the CD and EoR aim to measure the time evolution of the spatial fluctuations of the signal at different scales (e.g., Gehlot et al. 2020;Mertens et al. 2020;Trott et al. 2020;Abdurashidova et al. 2023;Munshi, S. et al. 2024).The upcoming Square Kilometre Array (SKA) 6 will have a higher sensitivity and will further produce tomographic images of the distribution of the Hi signal on the sky (Mellema et al. 2015;Ghara et al. 2017).The second type of radio observations such as EDGES (Bowman & Rogers 2010), EDGES2 (Monsalve et al. 2017;Bowman et al. 2018), SARAS (Patra et al. 2015), SARAS2 (Singh et al. 2017), LEDA (Price et al. 2018) use single radio antenna and aim to measure the redshift evolution of the sky averaged Hi 21-cm signal from these epochs.
The redshifted Hi 21-cm signal from the CD and EoR is enriched with information about the first sources of our Universe.The measurements of the redshifted Hi signal from these epochs are often used to infer the properties of the sources formed during these epochs (e.g., Ghara et al. 2020;Greig et al. 2021;Mondal et al. 2020;Abdurashidova et al. 2022), nature of dark matter (e.g., Barkana 2018;Fialkov et al. 2018;Muñoz & Loeb 2018;Nebrin et al. 2019;Ghara & Mellema 2020;Ghara et al. 2022), etc.Furthermore, the same observations can be used to infer the properties of the IGM such as the thermal and ionization states and morphology of the 21-cm signal maps (e.g., Ghara et al. 2020;Ghara & Choudhury 2020;Ghara et al. 2021).The EoR 21-cm signal is expected to be enriched with information about the ionized regions while the CD 21-cm signal is expected to carry information about the characteristics of the heated/emission regions.The heated regions around the early X-ray emitting sources such as the HMXBs, mini-Quasars are expected to appear as emissions in the Cosmic microwave background (CMB) subtracted 21-cm signal maps while the cold regions will appear as absorption signal (see e.g., Ghara et al. 2016Ghara et al. , 2017)).
The radio interferometer-based observations probe multi-scale properties of the signal and thus, are expected to contain more information about the states of the IGM.The fluctuation in the observed 21-cm signal by the radio interferometers depends on the morphology and distribution of the ionized/emission regions.Thus it is important to study the morphology of the ionized/emission regions of the signal to understand the connection between the measured statistical quantities of the signal such as the power spectrum or the sky-averaged signal.Recently several studies using methods such as granulometry method (Kakiichi et al. 2017), Mean-free path (Mesinger & Furlanetto 2007), image segmentation method (Giri et al. 2018), percolation analyses (Iliev et al. 2006b(Iliev et al. , 2014a;;Furlanetto & Oh 2016) and the Minkowski functionals (Friedrich et al. 2011a;Yoshiura et al. 2017;Kapahtia et al. 2021), etc. have made a significant effort in understanding the morphology of the ionized and neutral regions during the EoR.These studies indicate that the characteristics of the ionized/neutral regions e.g., the probability density functions (PDF) of the bubble size distribution can be used to distinguish different reionization scenarios and thus to constrain the EoR source parameters, etc.
In this work, we focus on the morphological analysis of the Hi 21-cm signal from the CD.Similar to the ionized/neutral regions which are the focus of morphology studies of the 21-cm signal from the EoR (e.g., Mesinger & Furlanetto 2007;Kakiichi et al. 2017;Bag et al. 2018;Kapahtia et al. 2021), we consider the emis-6 http://www.skatelescope.org/sion/absorption regions and study the morphology and distribution of these regions.Our study is based on simulated 21-cm signal maps and surfgen2 algorithm (Bag et al. 2019(Bag et al. , 2023) ) which assesses the geometry, morphology and topology of a field above/below certain thresholds.This algorithm models the surface of a region through the Marchingcube 33 triangulation (Chernyaev 1987;Lorensen & Cline 1995).The accuracy in measuring the Minkowski functionals and Shapefinders is much better compared to the other existing methods, such as the cell counting-based Crofton's formula (Crofton 1868).Previously, Bag et al. (2018Bag et al. ( , 2019)); Pathak et al. (2022); Dasgupta et al. (2023) studied the morphology of the ionized regions during the EoR using the same algorithm.These studies found that the time evolution of the Largest Cluster Statistics (LCS) of the ionized regions is a robust means to study the percolation of the ionized regions and can be used to distinguish between different reionization scenarios.Here, we study the evolution of the emission regions LCS and the percolation state for the emission regions of the 21-cm signal during the CD.In fact, we would like to point out that this paper is the first effort to use LCS for the morphological analysis of the CD 21-cm signal, all previous work with LCS has been limited to 21-cm signal coming from the EoR.
This paper is structured in the following way.We first describe the simulation used in this study in Section 2, then we discuss the methodology to characterise the emission regions in Section 3. Section 4 shows the results from this study.We discuss and conclude the findings of this study in Section 5.This study uses cosmological parameters Ω m = 0.27, Ω Λ = 0.73, Ω B = 0.044, ℎ = 0.7 (Wilkinson Microwave Anisotropy Probe (WMAP); Hinshaw et al. 2013) which are the same as used in the N-body simulation used in this study.

SIMULATION
In this section, we briefly present the simulation that leads towards 21-cm signal brightness temperature maps used for the morphological study.Here, we use grizzly (Ghara et al. 2015(Ghara et al. , 2018) ) code to simulate the redshifted 21-cm signal from the CD and EoR.The algorithm is based on one-dimensional radiative transfer (RT) and is an independent implementation of the previous 1D RT-based code bears (Thomas et al. 2009;Krause et al. 2018).
The differential brightness temperature ( b ) of the redshifted Hi 21-cm signal emitted from a region with angular position x and redshift  can be expressed as (see e.g., Madau et al. 1997) where the Hi 21-cm signal optical depth  b can be denoted as (Bharadwaj & Ali 2005;Datta et al. 2022), (2) Here,   is the brightness temperature of the radio background.In the absence of exotic sources of radio emission (e.g., see Mondal et al. 2020;Ghara et al. 2022),   is expected to be the CMB brightness temperature   () = 2.73 ×(1 + ) K.The quantities  HI and  B are the neutral fraction and the density contrast of hydrogen in that region, respectively.The quantity  S stands for the spin temperature of the neutral hydrogen which quantifies the relative From left to right panels correspond to redshifts 14.7, 13.5, 12.9 and 12.3 with heated fractions 0.05, 0.3, 0.6, and 0.9, respectively.These maps correspond to   =  CMB .population of hydrogen in the two hyperfine states of the hydrogen atom's ground state.
The grizzly algorithm uses dark matter halo catalogues and density and velocity fields on grids as derived from N-body simulations.Given these inputs, grizzly generates  b maps for a source model.Here, we use the outcome of a dark-matteronly N-body simulation from the results of the PRACE7 project PRACE4LOFAR.The simulation was performed using cubep 3 m (Harnois-Déraps et al. 2013) code with 6912 3 particles in a box of volume (500 ℎ −1 Mpc) 3 which set the particle mass resolution to 4.05 × 10 7 M ⊙ .Later, the density and velocity fields are smoothed into grids.We use 300 3 grids in this paper.The dark matter halo catalogues were generated on the fly with a spherical overdensity halo finder (Watson et al. 2013).These contain all halos with masses larger than 10 9 M ⊙ which is ≈ 25 times the dark matter particles resolution.The simulation saved the gridded fields and halo catalogues at different redshifts while the time gap between two successive redshifts is 11.4 mega-years.Details of the N-body simulation can be found in Dixon et al. (2016); Giri et al. (2019b).Here we use 59 redshift snapshots between  ∼ 20 − 6.5 to simulate the CD and EoR for the following source model.
The astrophysical sources dependence in Equation 1 comes through  HI and  S .Both these quantities depend on the emissivity of UV, X-ray, radio and Ly photons, as well as the population number density of the sources (see e.g., Ross et al. 2019;Reis et al. 2020).The source model adopted in this study assumes that each of the dark matter halos on the halo catalogues emits UV, X-ray and Ly photons.The emission rates of these photons are proportional to the stellar mass ( ★ ) of the sources which we assume to be linearly related to the dark matter halo mass  halo , i.e.,  ★ = Ω m  halo .We fix the star formation efficiency  ★ = 0.02 in this study.The emission rate of the UV photons per stellar mass   is controlled by the Ionization efficiency parameter  as The cosmic heating and reionization depends also on the minimum mass of dark matter halos ( min ) that contain radiating sources.Here, we have assumed that all dark matter halos with masses larger than  min = 10 9 M ⊙ host sources of UV, X-rays and Ly photons.We use  = 0.1 for which the EoR ends around  ∼ 6.5 8 .The X-ray spectral energy distribution of the source is assumed to be a power-law with energy  as  X () ∝  −  with spectral index  = 1.2.The emission rate of the X-ray photons per stellar mass  X =  X × 10 42 s −1 M −1 ⊙ where  X is the X-ray heating efficiency parameter.The X-ray band spans from 100 eV to 10 keV.We choose  X = 100 for the heating model considered in this study.The details of the spectral energy distribution and source model can be found in Ghara et al. (2015); Islam et al. (2019).These references also provide the details of the method to simulate the coeval cubes of neutral fraction  HI and  K at different redshifts.
The solid line of the left panel of Figure 1 shows the redshift evolution of the ionization fraction while the dashed curve shows the evolution of the heated fraction  heat i.e., the volume fraction of regions with gas temperature larger than the CMB brightness temperature.The Figure shows that the IGM is significantly heated by redshift ∼ 13 while the ionization fraction remains less than 0.1.Figure 1 also shows slices of the gas temperature distribution at three different stages of heating.From left to right  K maps represent redshifts 14.7, 14.3 and 13.5 with  heat 0.05, 0.1 and 0.3, respectively.The slices show that significant overlap between heated regions already started as early as a stage with redshift 14.3.The Ly photon flux fields at different redshifts are generated assuming a 1/ 2 fall of the number density of Ly photons with radial distance  from the source.These Ly photon flux fields are later used to generate the  S fields.These  HI ,  S cubes and the density fields are then used to generate the  b maps using Equation 1. Finally, we incorporate the effect of the peculiar velocities of the gas or socalled redshift-space distortion using a cell movement method (see e.g., Mao et al. 2012;Ross et al. 2021).Figure 2 shows the  b slices at four different stages of the CD.The size of the heated/emission regions (with  b ≥ 0) increases as heating progresses.The panels also show that, above  heat ≳ 0.5, the distribution of the absorption regions (with  b < 0) is more relevant compared to the emission regions.

MORPHOLOGY ANALYSIS
There exists no unique way to characterize the morphology and size distribution of a complex spatial structure such as the distribution of  b emission regions in the IGM.Here we use surfgen2 algorithm (Bag et al. 2019(Bag et al. , 2023) ) to assess the morphology of individual heated regions in the 3D coeval  b boxes generated using grizzly.Given a coeval cube of the  b , we first generate a binary field  b bin (x, ) where  b bin (x, ) = 1 if  b (x, ) ≥ 0, and zero otherwise.These binary maps are then used in surfgen2 for the morphological analysis of the emission regions.The code implements periodic boundary conditions on the simulated  bin b (x, ) boxes and uses the friends-of-friends (FoF) algorithm to identify individual emission and absorption regions.Figure 3 shows slices of such binary fields at different stages of our simulated CD as described in the previous section.The white regions in the slices represent emission regions for   =  CMB .For  heat > 0.5, the size distribution of the absorption regions is more meaningful than the emission regions.This is due to the inside-out nature of the heating around the X-ray sources.Clearly, the individual emission regions grow in size over time and overlap to create a large region that dominates the total emission segment volume, while the morphology of the absorption regions towards the end of the CD appears different in nature.Visually the morphology of these binary emission regions is similar to the morphology of ionized regions during an inside-out reionization scenario (see e.g., Bag et al. 2018;Pathak et al. 2022).
The calculation of the Minkowski functionals and thereafter the Shapefinders of a connected (emission/absorption) region in surfgen2 is based on modelling the surface of that region through the advanced Marching Cube 33 triangulation scheme (Chernyaev 1987;Lorensen & Cline 1995).We refer the readers to Sheth et al. (2003); Bag et al. (2019) papers for the details of the surfgen2 algorithm.The algorithm first determines the four basic Minkowski functionals (Mecke et al. 1994) that well describe the morphology of an individual region.These are (i) volume: , (ii) surface area: , (iii) integrated mean curvature (IMC): (iv) integrated Gaussian curvature or Euler characteristic: Here  1 and  2 represent the two principle curvatures at the surface element  on the surface .The integrated Gaussian curvature  is related to the genus () of a closed surface as  = 1 − /2.Physically  ≡ (no. of tunnels) − (no. of isolated surfaces) + 1.Note that while  and  of the emission and absorption regions will always be positive, IMC and  might also be negative in certain cases.
To measure the shape of a region directly, we calculate the following three ratios of the Minkowski functionals, (5) Note that each of ,  and  are of the dimension of length and are spherically normalized, i.e.  = (4/3) .These ratios were first introduced in Sahni et al. (1998) as 'Shapefinders' since they assess the extent of a region in three dimensions.Here, ,  and  represent the length, breadth and thickness of a region, respectively.A sphere of radius  will have  =  =  = .These quantities should satisfy natural ordering  ≥  ≥ .However, this ordering might not be maintained naturally in some rare cases.In those cases, we define the largest value among the three in Equation 5as  and the smallest value as , to maintain the order  ≥  ≥ .Further, the planarity () and filamentarity () of a region are defined as the following dimensionless quantities based on the Shapefinders (see Sahni et al. 1998), As their names suggest, they measure the degree of planarity and filamentarity of a region.By definition,  and  values are bounded by 0 and 1.Note that, for a spherical object,  =  =  leads to  = 0 = .For a planar object, such as a sheet,  ∼  ≫  results in  ≫ .On the other hand,  ∼  ≪  for a filament and thus  ≫ .
In this work, we employ the above-mentioned Minkowski functionals and Shapefinders to statistically study the morphology of the emission and absorption regions of our simulated CD 21-cm signal.In addition, we consider the morphological properties of the largest emission and absorption regions which provides crucial information about the percolation state of these regions.

RESULTS
In this paper, we explore the morphological distribution of the emission and absorption regions of the 21-cm signal during the CD.We    (2022), the fall of the signal amplitude at the boundary of the emission regions where  b = 0 is not very sharp.Furthermore, the nature of the signal, either for emission or absorption, depends crucially on the radio background temperature   which is uncertain at the frequencies of our interest (e.g., Fixsen et al. 2011).Thus, we choose a number of different thresholds  b,thre on the  b (x, ) maps to study the morphology of the regions with  b larger than  b,thre which we refer to as emission segment.In addition, we also consider the complementary regions, i.e. the absorption regions with  b <  b,thre .We have seen from Figures 2 and 3 that the distribution of these absorption regions is more relevant for  heat ≳ 0.5.Different  b,thre values also represent different strengths of radio background at 1420 MHz at the redshift of our interest.We parameterize the  b,thre in terms of radio background temperature   =  T × 2.725(1 + ) K with the threshold parameter  T .Note that, similar to the majority of the 21-cm signal study,  T = 1 is the default threshold for our study and it will be the same if not mentioned otherwise.
The number density and size of the emission regions vary with the radio background temperature   and thus also with threshold parameter  T .The left panel of Figure 4 shows the evolution of the number density (Φ) of isolated emission (dashed curves) and absorption (solid curves) regions with  heat for four different  T values.The magenta circular points in the left panel of Figure 4 mark different CD stages where the morphology of the emission and absorption regions have been estimated.Φ of the emission regions first increases with  heat until  heat reaches ∼ 0.1 due to the formation of new heated regions.For  heat ≳ 0.1, Φ of the emission regions gradually falls before sharply falling to 1 when  heat → 1.This happens due to overlaps between individual emission regions.Both the maximum value of emission Φ and the corresponding  heat (at which the maximum occurs) remain smaller for a smaller value of  T .This is expected as by increasing  T , a large emission region at a later stage of heating can break into multiple smaller emission regions.The fall of Φ for emission regions is faster for a smaller value of  T .In fact, Φ changes very little in the range 0.1 ≲  heat ≲ 0.9 for  T = 10.This implies that the overlap between the emissions regions happens at a later stage for a larger  T value.
The number density of absorption regions, on the other hand, gradually increases with  heat since larger absorption regions break into smaller ones as heating progresses.Contrary to the Φ of the  emission regions, the Φ of the absorption regions (solid curves of Figure 4) decreases with the increase of  T .However, Φ of absorption regions sharply drops to zero at the end of the heating of the IGM, as expected.
The filling factor (which is defined as the volume fraction of emission or absorption segment), FF, of emission or absorption regions also depends on  T .The right panel of Figure 4 shows how FFs of the emission (dashed) and absorption (solid) regions evolve with  heat for different  T values.As expected, FF of the emission regions (FF emi ) is roughly the same as the  heat for  T = 1.For a certain  heat , FF emi decreases for a larger value of  T .On contrary to FF emi , FF of the absorption regions (FF abs ) decreases as heating progresses since FF abs = 1 − FF emi .

Largest cluster statistics
The 'largest cluster statistics' (LCS) for a segment (emission or absorption) is defined as the volume fraction enclosed by the largest region of that segment, LCS = volume of the largest emission/absorption region total volume of all the emission/absorption regions .
It has been established that the LCS carries robust information about the percolation process in different fields at cosmological scales, e.g. the large-scale matter distribution in the Universe (Shandarin 1983;Klypin & Shandarin 1993;Yess & Shandarin 1996) and the neutral hydrogen field during the EoR (Iliev et al. 2006b(Iliev et al. , 2014a;;Furlanetto & Oh 2016;Bag et al. 2018).For the ionization bubbles, the percolation happens early with volume averaged ionization fraction  HII ∼ 0.1 (see e.g., Furlanetto & Oh 2016;Bag et al. 2018).Here, we study the characteristics of the largest emission and absorption regions in the IGM at different stages of CD and the percolation transitions in the respective segment.
In the left panel of Figure 5, we present the evolution of LCS of the emission and absorption regions as a function of  heat for the four different choices of  T .At the initial stages of the heating (e.g., with FF emi ≲ 0.1), the emission regions are fairly isolated with nominal overlaps with each other.As heating progresses, many new emission regions appear around the X-ray emitting sources and the emission regions themselves grow in size resulting in a gradual increase in the LCS of the emission regions with  heat (see the dashed curves in the left panel of Figure 5).Soon the emission regions start to merge significantly.Due to the vigorous merging of smaller regions, an extremely large connected emission region that extends throughout the IGM appears abruptly at this stage.In our case, it stretches throughout the simulation box.We refer to this period as the 'percolation transition' in the emission regions.The LCS of the emission regions increases sharply during the percolation transition since the largest region abruptly accumulates most of the emission segment volume.The critical  heat at which percolation transition takes place depends on  T values.However, percolation transition roughly occurs at similar values of the filling factor FF emi ∼ 0.15 (marked by the vertical dotted line in the right panel of Figure 5) for all the threshold values chosen (see the dashed curves of the right panel of Figure 5).
The LCS of the absorption regions, on the other hand, remained ≈ 1 for a long time before sharply dropping to zero around  heat ∼ 1.This suggests that the absorption segment in the IGM remains percolated for a large extent of the CD and only breaks into isolated absorption regions towards the end period of the CD.When plotted against FF, the LCS of the absorption regions shows a sharp decrease around FF emi ∼ 0.95 for all the thresholds chosen.The evolution of the LCS of the emission regions with FF emi is almost identical for the different choices of  T .

Evolution of the largest emission and absorption regions
Next, we study the evolution of shape and morphology of the largest region, separately for the emission and absorption segments, in terms of Shapefinders.Note that we stick to  T = 1 for the rest of the study.This threshold corresponds to the CMB as the radio background.The left panel of Figure 6 shows the evolution of the first three Minkowski functionals ,  and IMC of the largest emission (dashed) and absorption (solid) regions as functions of FF emi for  T = 1.All these three quantities evolve rapidly around the percolation transition of the corresponding segment (at FF emi ≈ 0.15 and 0.95 for the emission and absorption segments respectively).In fact, this panel demonstrates how the largest emission/absorption region grows rapidly at the onset of their respective percolation transitions and this is manifested in the evolution of the LCS in Figure 5.
Let us now focus on the largest emission region.As expected, its  always increases with the increase of FF emi .Its , on the other hand, increases with FF emi until FF emi ∼ 0.6 and slowly decrease thereafter.This shows at stages with FF emi ≳ 0.6, the largest emission region in the simulation box becomes less complex in shape with less amount of absorption regions surface inside it, although its volume increases with time.While  and  are always positive quantities, IMC can be negative also.The change in the sign of the IMC occurs around FF emi ∼ 0.65 for this model CD.The negative part of the IMC is shown in the light-blue colour in the left panel of Figure 6.For FF emi ≳ 0.7, one can think of the IGM as a distribution of several absorption regions in an emission region background.At these stages, the sum of IMCs of all the absorption regions is positive which makes the IMC of the largest emission region negative.The evolution of the largest absorption region exhibits similar behaviour but in the opposite direction along the FF emi axis.
The middle panel of Figure 6 shows the evolution of the three Shapefinders -'thickness' (), 'breadth' () and 'length' () -of the largest emission and absorption regions as functions of FF emi , again for  T = 1.The 'length' () of the largest emission or absorption region increases rapidly at the onset of the percolation transition (in the respective segment, i.e.FF emi ≈ 0.15 and 0.95 for the emission and absorption segment respectively) as compared to its  and .
The evolution of  of the largest emission region merely follows the evolution of its |IMC| (see Equation 5).
The right panel of Figure 6 shows the evolution of morphological quantities,  and , and the topology, , of the largest emission (dashed curves) and absorption (solid curves) regions. of the largest emission region is more than 2 orders of magnitude larger than its  around the percolation transition stage.This leads to the largest emission region being highly filamentary with  ∼ 1 (see Equation 6) at percolation.The same is also applicable to the largest absorption region.Also, the genus values of both the largest emission and absorption regions increase at the onset of the percolation transition in the respective segments (see the blue curves), and remain stable at post-percolation stages.This implies that the largest regions become more multiply connected with complex topology as they grow and percolate.Note that, at some stage, well after percolation, the genus value of the largest region (emission or absorption) abruptly becomes negative.In this phase where the genus is negative, we show its absolute value by the light-blue coloured curves (solid and dashed).The sign flip, which indicates the change in the topology of the largest region, manifests in negative IMC values (see left panel of Figure 6) which in turn leads to the sharp dip/rise in different Shapefinders (as evident from the middle panel and explained above) at this stage.Indeed, well after percolation in the emission (absorption) segment, many isolated small pockets of absorption (emission) regions emerge inside the percolated largest emission (absorption) region that brings down the latter's genus value and eventually rendering it negative.

LCS comparison between the emission and ionized regions
Next, we compare the LCS of the CD emission regions with the LCS of the EoR ionized regions.We use outputs from the same with grizzly simulation to study the LCS of the ionized regions.
The ionization history is presented in the left panel of Figure 1.We consider binary field of ionization fraction  bin HII (x, ) between redshifts 20 and 6.5 where  bin HII (x, ) = 1 if  HII (x, ) ≥ 0.5.These binary ionization fields are then used as inputs to surfgen2 for estimating the LCS of the ionized regions.We present the LCS as a function of the filling factor of the ionized regions (FF HII ) in the left panel of Figure 7.The percolation of the ionized regions occurs at a stage with  HII ∼ 0.15 which is consistent with studies such as Iliev et al. (2006b); Chardin et al. (2012);Furlanetto & Oh (2016); Bag et al. (2018).This suggests that the sharp increase of LCS of the emission regions at FF emi ∼ 0.15 as seen in the left panel of Figure 5 is similar to the percolation case of the ionized regions.
The second left panel of Figure 7 shows the evolution of ,  and IMC of the largest ionized region as functions of FF HII .All these three quantities evolve rapidly around FF HII ∼ 0.15.These are consistent with previous studies such as Bag et al. (2018).Similar to the emission region's LCS case (see the left panel of Figure 6) these quantities evolve sharply around the percolation transition of the ionized segment.Similar to the sign change for the IMC of the largest emission region (see the left panel of Figure 6), we also see a sign change for the IMC of the ionized regions as shown in the second left panel of Figure 7. Overall, the evolution of ,  and IMC of the largest emission region around the percolation transition is similar to that of the same quantities of the largest ionized region in the ionization field.
The third left panel of Figure 7 shows the evolution of ,  and  of the largest ionized region as functions of FF HII while the right panel shows the evolution of ,  and .The evolution of ,  and  of the largest ionized region in our simulation box as functions of FF HII is consistent with previous studies such as Bag et al. (2018); Pathak et al. (2022); Dasgupta et al. (2023).These also highlight the filamentary nature of the largest ionized region during ionization percolation and agree well with the evolution of the largest emission region during the CD as we have seen in section 4.2.

Size Distribution of the emission and absorption regions
The left panels of Figure 8  The size distributions in the top panels are estimated using the mean-free path method, while the bottom panels are for the granulometry method.Here,  is the number of rays used in the MFP method whose sizes are between  and  + .On the other hand, the quantity  (< ) from the granulometry method is the volume fraction of the emission (or absorption) regions whose size is smaller than a radius .
inside the emission (or absorption) regions in random directions and record the lengths of the rays until those reach the edge of the regions.The top panels of Figure 9 show the probability density function (PDF) / of the recorded ray lengths at different stages of the CD.Here,  is the number of rays whose sizes are between  and  + .The top left panel shows the PDFs for the emission regions while the top right panel represents the absorption regions.The overlap of the smaller emission regions at a later stage of heating is prominent by the characteristic scales or the peaks of the PDFs.The peak of the PDFs of the emission regions moves toward the larger values as the heating progresses.The PDFs show log-normal type behaviour which is similar to the distribution of the ionized regions as well (see e.g., Mesinger & Furlanetto 2007).The size distribution of the absorption regions also shows a log-normal type nature and is similar to the distribution of the neutral regions during the EoR (see e.g., Giri et al. 2019a).
The bottom panels of Figure 9 show the PDFs  (< )/ of the size distribution as obtained by using a granulometry method.Here, the quantity  (< ) is the volume fraction of the emission (or absorption) regions whose size is smaller than a radius .We follow the method of Kakiichi et al. (2017) which is based on Minkowski subtraction and addition steps to estimate the PDFs.Note that, this method is based on applying spherical filters, and is sensitive to the smallest dimension of a complex-shaped region.
Thus, the algorithm does not find any structure with a radius as large as the largest ray length found in the MFP method.In our case, e.g., for FF ≲ 0.6 the granulometry method does not find any spherical structure of the emission or absorption regions with a radius more than ∼ 20 ℎ −1 Mpc.

Shape distribution of emission and absorption regions
In this subsection, we study how the shape (together with topology and morphology) of the emission and absorption regions is distributed over their volume and how these distributions evolve over time.The top panels of Figure 10 show the distributions of the ,  and  of the emission (dashed curves) and absorption (solid) regions as functions of their volume at different states of the CD.For both types of regions,  is significantly larger compared to their  and  while the difference increases with the increase of volume of the individual regions.E.g., for FF ∼ 0.1, / or / is ∼ 1 − 2 at  ∼ 10 2 (ℎ −1 Mpc) 3 , while the ratios increases to ∼ 20 for  ∼ 10 4 (ℎ −1 Mpc) 3 .This implies that small-volume emission/absorption regions are more spherical while the larger regions are more filament-type.The same conclusions can be made from the bottom panels of Figure 10 as well which show the distribution of ,  and  of the emission (dashed curves) and absorption (solid) regions as functions of their volume for different FF.While  of the emission regions is ∼ 0.2 for  ∼ 10 2 (ℎ −1 Mpc) 3 , it increases to ∼ 1 at  ∼ 2 × 10 3 (ℎ −1 Mpc) 3 .For the same range of , the change of  is smaller compared to that of .From the genus distribution (blue dashed and solid curves) it is evident that larger emission/absorption regions tend to be more multiply connected with more complex topology.Strikingly, the shape (and morphology) distributions of the absorption regions remain similar to that of the emission regions when both segments are compared at the same filling factor.

SUMMARY & DISCUSSION
The spatial fluctuation of the redshifted Hi 21-cm signal from 400 million years to the first billion years of our Universe's history carries information about the first radiating sources, the ionization and the thermal states of the IGM.The observations of this signal from the Epoch of Reionization (EoR) and the Cosmic Dawn (CD) using the world's largest radio interferometers have the ability to reveal many unknown facts about these epochs.Theoretical studies suggest that the spatial fluctuations in the signal from the Cosmic Dawn depend on the distribution and properties of the heating sources.An alternative interest focuses on the distribution of the emission regions around the heating sources.Such a morphological study is important to understand the size and distribution of such emission regions or the remaining part of the sky which is seen in absorption.
In this study, we focus on the shape and morphology of the emission and absorption regions during the Cosmic Dawn.This study is based on a simulation of the Cosmic Dawn 21-cm signal using grizzly (Ghara et al. 2015) code and analyzing it using surfgen2 algorithm (Bag et al. 2019) to produce the geometrical properties of the regions.The main findings of our study are the following.
• If the radio background is the CMB, the number density of the emission regions gets a maximum at a stage when the volume fraction of regions with gas temperature larger than the CMB temperature (heated fraction  heat ) is ∼ 0.1.In case the radio background is () of the emission (dashed) and the absorption (solid) regions functions of volume () for the threshold parameter  T = 1 at different stages of the CD.The bottom panels show the 'planarity' (), 'filamentarity' () and 'genus' () of the emission and the absorption regions as functions of  for  = 1 at different stages of the CD.FF stands for the filling factor corresponding to the type of region, i.e.FF stands for FF emi (FF abs ) when we consider the emission (absorption) regions.
stronger than CMB, the peak of the number density vs  heat shifts towards a larger  heat value and can occur between  heat = 0.1 and 0.9 depending on the radio background strength.The peak of the number density of the absorption regions as a function of  heat occurs at a stage when  heat ≳ 0.9.
• The percolation transition in the emission regions roughly occurs at a filling factor (defined as the volume of the emission segments) of FF emi ≈∼ 0.15.For the absorption regions, the percolation transition occurs at the filling factor FF abs ≈ 0.05 which corresponds to FF emi ≈∼ 0.95.The evolution of the 'largest cluster statistics' LCS (defined as the ratio of the largest region volume and the total volume of all regions with the same type) of the emission, as well as the absorption regions, as a function of the filling factor, are roughly identical for the different radio background levels chosen in this study.
• The overall evolution of the Minkowski functionals -volume (), surface area (), integrated mean curvature (IMC) and genus () -and Shapefinders -thickness (), breadth (), length (), planarity () and filamentarity () -of the largest emission region (when studied functions of the filling factor of the emission regions) during the CD is found to be similar to the evolution of the same characteristics of the largest ionized region (when studied against ionized filling factor) during the EoR.
• The number density of emission regions / as a function of their volume , estimated from the surfgen2 algorithm, shows a rough power-law dependence on the volume, i.e., / () ∝   with  ≈ −2.We also see similar behaviour for the number density of the absorption regions while  ≈ −1.6 for that.After the percolation transition, the largest region appears as a disjoint part of the size distribution of these regions.
• The size distribution of the emission regions when estimated using the mean free path method shows a log-normal type behaviour as a function of the length.A similar type of size distribution can be seen for the ionized regions during the EoR.On the other hand, the size distribution of the absorption regions is similar to the size distribution of the neutral regions during the EoR.A similar conclusion can also be made when the size distribution is estimated using the granulometry method.
• We find that the distributions of the Shapefinders (and genus) of the emission regions, with respect to their volumes, are quite similar to those of the absorption regions when both segments are compared at the same filling factor.In particular,  of both emission and absorption type regions is significantly larger compared to  and .This difference increases with the increase of the volume of these regions resulting in higher values .Therefore, we find that larger emission/absorption regions tend to be more filamentary with higher values of  as well.
The morphology and distribution of both the EoR ionized regions and CD emission regions are crucially dependent on the source properties and the number density of sources during these epochs.For example, a larger value of  min (i.e., the minimum mass of dark matter halos that hosts sources) will result in a smaller number density of X-ray emitting sources and a more patchy heating scenario.The conclusions of the work are based on a single scenario of CD.Thus, it will be important to investigate the robustness of our conclusion for changing CD scenarios.For the reionization case, as shown in Pathak et al. (2022), the percolation transition in the ionized regions (as probed by the evolution of LCS of the ionized regions) and the redshift evolution of the Shapefinders can constrain reionization scenarios and shed light on the ionizing sources.Similarly, the redshift evolution of the LCS and Shapefinders in emission regions might distinguish the fundamentally different heating scenarios and put constraints on the properties of the heating sources.We plan to investigate this in a future study.
The simulation volume ∼ (500 ℎ −1 Mpc) 3 as used in this work is much larger than the cosmological scale (∼ 100 − 200 ℎ −1 Mpc) beyond which the cosmological principles of isotropy and homogeneity hold for the EoR and CD 21-cm signal (Iliev et al. 2014b;Kaur et al. 2020).This ensures that the evolution of the LCS, which detects the percolation transition in the field unambiguously, is independent of the choice of the simulation volume (see Pathak et al. 2022, in the context of EoR).Note that the results of geometrical tools including ours might depend on the resolution of the 21-cm signal fields.Pathak et al. (2022) investigated the effect of resolution on LCS of ionized regions and found that the LCS results are consistent for different choices of simulation resolutions.We plan to study the effect of simulation resolution on the LCS of the CD 21-cm signal in a follow-up work.
This study has been done in a very idealistic situation.In reality, the 21-cm signal as observed by the radio telescopes will be contaminated by the astrophysical foreground, system noise, radio frequency interference (RFI), etc.Even for the most optimistic scenario with a full  coverage of the instrument, an efficient method to mitigate the foregrounds, RFI, etc., a perfect calibration technique, the presence of the system noise from a realistic observation (e.g., 1000 hours with the SKA) will affect the recovery of such high-resolution images as we have used in this study.Thus, the recovered size distribution will be affected.Recently, Dasgupta et al. (2023) considered mock SKA1-low observation and studied LCS of EoR 21-cm signal in the presence of telescope noise and other telescope-related effects such as image distortions due to the synthesized beam.The study found that these effects introduce a bias in the redshift evolution of the LCS, while the percolation transition in the field can be affirmatively recovered even with some calibration error.We expect similar results in the case of CD 21-cm observations as well.Studying the impact of the noise of the system on the recovery of the size distribution of the CD emission regions is beyond the scope of this work.We will address this issue in a follow-up paper.

Figure 1 .Figure 2 .
Figure 1.Ionization and heating history of the simulated Cosmic Dawn and reionization scenario.The left panel shows the evolution of the volume fraction of the ionized region (i.e., volume averaged ionization fraction xHII ) and heated regions (i.e., heated fraction  heat ) of our simulated Cosmic Dawn (CD) and Epoch of Reionization (EoR) as functions of redshift.We define a 'heated region' as a region with a gas temperature larger than the CMB temperature.The left to right two-dimensional maps show the gas temperature at three different stages of the CD with heated fractions 0.05, 0.1 and 0.3, respectively.The corresponding redshifts are 14.7, 14.3 and 13.5, respectively.The colour bar shows the temperature in Kelvin.

Figure 4 .
Figure 4. Evolution of the number density of emission and absorption regions (left panel) and their filling factors FF (right panel) as functions of the heated fraction, i.e., the volume fraction of regions with gas temperature larger than  CMB .Different types of lines represent different threshold parameter (  T ) values for the radio background temperature as chosen in this study.The magenta circular dots mark different stages of the CD where the morphology of the emission regions has been estimated.The emission and absorption regions at different stages of the CD are determined using surfgen2 algorithm.

Figure 5 .
Figure 5. Evolution of the Largest cluster statistics (LCS) as a function of the heated fraction (left panel) and filling factor (right panel) of the emission regions.The dashed curves represent emission regions while the solid curves represent the absorption regions of the 21-cm signal  b coeval cube.The magenta circular points mark different stages of CD when the morphology of the emission and absorption regions has been estimated.The largest emission and absorption regions at different redshifts of the CD are determined using surfgen2 algorithm.

Figure 6 .Figure 7 .
Figure6.Evolution of different shape defining quantities of the largest region of the emission (dashed) and absorption (solid) segments in the IGM as functions of filling factor of the emission regions correspond to  T = 1.The x-axis is equivalent to the heated fraction of the simulation.The left panel shows the volume (), surface area () and integrated mean curvature (IMC) of the largest emission and absorption region.The middle panel shows the thickness (), breadth () and length () of the largest regions, while the rightmost panel shows their planarity (), filamentarity () and genus ().

Figure 9 .
Figure 9. Size distribution of the emission (left panels) and absorption (right panels) regions during the CD.The size distributions in the top panels are estimated using the mean-free path method, while the bottom panels are for the granulometry method.Here,  is the number of rays used in the MFP method whose sizes are between  and  + .On the other hand, the quantity  (< ) from the granulometry method is the volume fraction of the emission (or absorption) regions whose size is smaller than a radius .

Figure 10 .
Figure10.Distribution of the shapefinders at different stages of Cosmic Dawn.The top panels show the Shapefinders 'length' (), 'breadth' () and 'thickness' () of the emission (dashed) and the absorption (solid) regions functions of volume () for the threshold parameter  T = 1 at different stages of the CD.The bottom panels show the 'planarity' (), 'filamentarity' () and 'genus' () of the emission and the absorption regions as functions of  for  = 1 at different stages of the CD.FF stands for the filling factor corresponding to the type of region, i.e.FF stands for FF emi (FF abs ) when we consider the emission (absorption) regions.