with compiled pdf

.gitignore

@@ -1,3 +1,7 @@
 *.pdf
 # but keep the pdfs in the "assets" folder
 !assets/*
+
+# finally - allow the main.pdf to be commited for completeness
+!main.pdf
+

abstract.typ

@@ -1,4 +1,8 @@
 #import "helpers.typ": *
-= Abstract
-We present an improved version of the #beorn framework, a semi-numerical simulation suite that generates maps of the cosmic dawn and the epoch of reionization. The refinements include a self-consistent treatment of the evolution of individual galaxies, a parametrization of stochasticity of the mass accretion rate, and a general optimization that allows for speedier simulation runs.
-We validate the improved version of the suite against ??. We employ the Thesan-Dark simulation to inder halo mass history and demonstrate the effect of this more detailed treatment.
+
+We present #beorn, a semi-numerical simulation suite that uses a simplified description of radiation sources to efficiently generate maps of the cosmic dawn and the epoch of reionization.
+In its original formulation, the framework assumed exponential halo growth with a fixed accretion rate. We show how this assumption introduces systematic inaccuracies in the predicted observables of reionization. In this work, we implement a more realistic description of halo mass growth that follows the mass accretion history of individual halos.
+We achieve this consistent treatment of halo growth by incorporating individual accretion histories extracted from the merger trees provided by the #thesan simulations.
+This required an extensive refactoring of the #beorn code base, redesigning the computation of the radiative source properties to accommodate variable growth while ensuring scalability and computational efficiency.
+We describe these modifications in detail and validate the improved suite against the original implementation.
+We demonstrate that the refined treatment alters the timing and morphology properties of the reionization. Finally, we discuss how this new version of #beorn provides a more robust platform for exploring reionization scenarios, and we outline planned extensions and applications.

acknowledgements.typ

@@ -2,3 +2,7 @@
 #heading(numbering: none, level: 1, outlined: false)[Acknowledgements]
 
 We would like to thank Sambit Giri and Yu-Hsiu Huang for their valuable input and helpful discussions during the development of this code. Their expertise and insights have significantly contributed to the robustness and accuracy of the simulation suite.
+
+
+
+This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID `uzh27 - 32205` on Alps.

appendix.typ

@@ -13,11 +13,12 @@
 #figure(
   image_cell(notebook, cell_id: "halo_mass_functions"),
   caption: [
+    The halo mass functions of the #thesandark 1 (blue) and #thesandark 2 (red) plotted at different redshifts. The error bars indicate the Poisson error in each mass bin.
   ]
 ) <fig:halo_mass_functions>
 
-// TODO - comment
+We compare the halo mass functions of the two dark-matter-only simulations #thesandark 1 and #thesandark 2. As shown in @fig:halo_mass_functions, the higher resolution of the former allows to resolve smaller halos more accurately. The two simulations only agree above a mass of about $5 dot 10^9 M_dot.circle$ and then yield the same halo mass function. This is consistent with our observations in the validation (@validation) of both simulations and poses a limitation on the minimum halo mass that can be resolved in #beorn when using merger trees from #thesandark 2.
 
 
 == B - Generation of the cover image
-The cover image of this report has been generated using #beorn. From a simulation run with a grid resolution of $512^3$ cells, a slice of the brightness temperature map has been extracted at $z = 8.07$. The slice shows an emission due to the spin temperature being higher than the CMB temperature. The first ionization bubbles appear as dark patches where the ionized hydrogen does not contribute to the 21 cm signal.
+The cover image of this report has been generated using #beorn. From a simulation run with a grid resolution of $512^3$ cells, a slice of the brightness temperature map has been extracted at $z = 8.07$. The slice shows an emission due to the spin temperature being higher than the CMB temperature. The first ionization bubbles appear as dark patches where the ionized hydrogen does not contribute to the 21-cm signal.

conclusion.typ

@@ -3,16 +3,18 @@
 
 = Conclusion <conclusion>
 
-With the highly anticipated detection of the reionization signal and upcoming observations of the conditions of the intergalactic medium during the cosmic dawn, the interpretation of these observations requires accurate predictions from theoretical models and simulations. #beorn
+With the highly anticipated detection of the reionization signal and upcoming observations of the conditions of the IGM during the cosmic dawn, the interpretation of these observations requires accurate predictions from theoretical models and simulations. #beorn
 by @Schaeffer_2023
 is a semi-numerical simulation framework that implements the halo model of reionization
 #cite(<schneider2023cosmologicalforecast21cmpower>, form: "normal") #cite(<Schneider_2021>, form: "normal").
-It uses flux profiles to express the emission of radiation by sources in terms of their host halo to simulate the reionization on large volumes and to obtain predictions for the 21-cm signal. It excels in its computational efficiency and flexibility, allowing for fast and flexible execution.
+It uses flux profiles to express the emission of radiation by sources in terms of their host halo to efficiently simulate the reionization on large volumes and to generate predictions for the 21-cm signal. It excels in its computational efficiency and flexibility, allowing for fast and flexible execution.
 
-We have presented an extension to #beorn that improves the physical accuracy by implementing a more consistent growth of galaxies based on the individual mass accretion histories of their host dark matter halo. We use the fact that the input data from the underlying #nbody simulation already includes constraints on the growth from the halo properties at different snapshots. Disregarding this information and instead assuming a fixed accretion rate for all halos is an oversimplification. The proof-of-concept implementation presented here leverages the halo history encoded in the merger trees of the #thesan simulation. More broadly, the updated framework is now better suited to incorporate more detailed growth simulations and can be easily extended to other simulations.
+We have presented an extension to #beorn that improves the physical accuracy by implementing a more consistent growth of galaxies based on the individual mass accretion histories of their host dark matter halo. We use the fact that the input data from the underlying #nbody simulation already includes constraints on the growth from the halo properties at different snapshots. Disregarding this information and instead assuming a fixed accretion rate for all halos is an oversimplification. The proof-of-concept implementation presented here leverages the halo history encoded in the merger trees of the #thesan simulation. More broadly, the updated framework is now better suited to incorporate more detailed growth simulations and can be easily extended to other simulations. We also refactored the simulation procedure and achieved a faster execution time by a factor of two. This further enhances the usability of #beorn for large parameter studies.
 
-After validating the new procedure we have shown that the consistent modeling of halo growth produces simulation outputs which have distinct features compared to simpler models. We compared map outputs direcly and also analyzed global quantities and their derived signal. The results are sensitive to the distribution of accretion rates, highlighting the importance of careful modeling of the halo growth.
+After validating the new procedure we have shown that the consistent modeling of halo growth produces simulation outputs which have distinct features compared to simpler models. We compared map outputs directly and also analyzed global quantities and their derived signal. The results are sensitive to the distribution of accretion rates, highlighting the importance of careful modeling of the halo growth.
+
+
+Works going beyond this proof-of-concept implementation should utilize more sophisticated history tracking that ensures the consistency of halo properties across mutlitple timesteps (e.g. the `rockstar` halo finder by @Behroozi_2012). We also highlighted the limitations incurred by the coarse mass resolution of the #thesan simulation, which is why subsequent research should be based on higher resolution simulations in order to benefit from accurate accretion rate matching down to the lowest halo masses.
+
+Furthermore, our analysis of halo growth shows that a simple modeling relying on the halo mass alone is insufficient. Many of the radiative properties in the halo model of reionization have similarly been expressed as a function of halo mass. Given the increasing evidence against simple mass-only models such as the stellar-to-halo-mass relation, the refinement of these models using stochasticity or additional halo properties is a promising avenue for future research. The impact of revised models, including our refinement to the halo growth will be subject to a future publication.
 
-Works going beyond this proof-of-concept implementation should utilize more sophisticated history tracking that ensures the consistency of halo properties accross mutlitple timesteps (e.g. the `rockstar` halo finder by @Behroozi_2012). Furthermore, our analysis of halo growth shows that a simple modeling with a mass dependent accretion is insufficient. Similar cases can be made for the other parameters that govern the emission of radiation from galaxies.
-// which ones??
-Investigations of the effect of stochasticity in the stellar-to-halo mass relation and the escape fraction of ionizing photons are promising direction for future research.

halo_mass_history.typ

@@ -3,8 +3,7 @@
 
 = Halo mass history <halo_mass_history>
 
-This section shows the impact of the halo growth on the resulting radiation profiles and motivates the need for a more precise treatment of the halo mass history. We show how to leverage simulation data for a refined simulation.
-// Don't like refined simulation
+This section shows the impact of the halo growth on the resulting radiation profiles and motivates the need for a more precise treatment of the halo mass history. We show how to leverage simulation data for an improved consistency during the simulation.
 
 
 == Modeling mass accretion
@@ -12,97 +11,64 @@ This section shows the impact of the halo growth on the resulting radiation prof
 As described in @hmreio the fundamental assumption of #beorn is the halo model of reionization by @schneider2023cosmologicalforecast21cmpower.
 // no need to recite?
 It describes how observables of reionization can be parametrized in terms of the halo mass and more specifically its rate of change since they are derived from the star formation rate expressed in @eq:star_formation_rate.
-In this simplified model, for a given star formation efficiency
-#footnote[
-  Note that the assumption of a fixed star formation efficiency or even an analytic expression as a function of halo mass is a simplification.
-  // Citation
-  The investigation of stochasticity has been subject to separate research (e.g. "missing").
-],
-the halo mass history is the single most impactful property besides the mass itself.
+In this simplified model, for a given star formation efficiency the halo mass history is the single most impactful property besides the mass itself.
 
 
-#beorn's goal is to provide simulations of the map-level contributions to the 21-cm signal, meaning that we cannot rely on a distribution of halo masses and accretion rates alone. Instead, #beorn leverages large scale N-body simulations to provide a spatial distribution of halos. For the first iteration halo growth was modeled through an exponential growth model
+#beorn's goal is to provide simulations of the map-level contributions to the 21-cm signal, meaning that we cannot rely on a distribution of halo masses and accretion rates alone. Instead, #beorn leverages large scale N-body simulations to provide a spatial distribution of halos. Halo growth is considered to follow the exponential growth model
 $
   M_"h" (z) = M_"h" (z_0) dot exp[-alpha (z - z_0)]
 $ <eq:exponential_growth>
-where $alpha = - dot(M_"h") / M_"h"$ is a free parameter describing the specific mass accretion rate. Following `@???`
-// TODODO
-a value of $alpha = 0.79$ was used as a fiducial value for all halos, independent of their mass or redshift. This meant that the requirements on the simulation data were minimal: Only a single halo catalog at a given redshift was required to generate a map at that redshift.
+where $alpha = - dot(M_"h") / M_"h"$ is a free parameter describing the specific mass accretion rate. Following @10.1093-mnras-stt1338
+a value of $alpha = 0.79$ was assumed for all halos, independent of their mass or redshift in the initial version. This means that the requirements for the simulation data were minimal: Only a single halo catalog at a given redshift was required to generate a map at that redshift.
 
-Using a simple exponential growth model is a significant simplification of the complex process of halo growth
-// maybe a citation
-but the most obvious
-// maybe a better word
-limitation is the assumption of a constant accretion rate $alpha$ for all halos, independently of their position, mass or redshift. In a realistic scenario we expect to observe a correlation with both the halo mass and redshift, in addition to the stochasticity of the accretion process. From a statistical perspective, this has been investigated by @Schneider_2021 who also consider a halo growth following the extended Press-Schechter formalism. This more detailed treatment shows that in particular small scales deviate from the simple exponential growth model. From a simulation perspective an even more precise treatment is possible since the growth history of each halo is already encoded in the successive snapshots of the N-body simulation. Ignoring this information introduces inconsistencies by painting halos using profiles that might not reflect their actual growth history.
+Using a simple exponential growth model is a significant simplification of the complex and time-sensitive process of halo growth (e.g. @McBride_2009). Another limitation is the assumption of a constant accretion rate $alpha$ for all halos, independently of their position, mass or redshift. In a realistic scenario we expect to observe significant stochasticity of the accretion process as well as systematic effects from the halo mass and the redshift. From a statistical perspective, this has been investigated by @Schneider_2021 who also consider a halo growth following the extended Press-Schechter formalism. This more detailed treatment shows that in particular small scales deviate from the simple exponential growth model. From a simulation perspective, an even more precise treatment is possible since the growth history of each halo is already encoded in the successive snapshots of the N-body simulation. Ignoring this information introduces inconsistencies by painting halos using profiles that might not reflect their actual growth history.
 
 
-// In a purely formal investigation where a qualitative prediction is derived from a well-defined halo mass distribution, the mass history is simply obtained as a direct derivation from the mass distribution. The simulations made by #beorn aim to provide 3D data that allows for quantitative conclusions. To this end a spatial distribution of the halo mass history is required, as provided by large scale simulations
-// #footnote[
-//   As described previously the halo model allows us to restrict the simulation to dark matter only, allowing for a more efficient simulation of the large scale structure.
-// ]
-// // Cite pkdgrav, Illustris, THESAN
-// .
-
-// @Schneider_2021 already compared exp growth to other models and found that following a more rigorous EPS approach show less growth at small masses.
-// // e.g. papers like "2309...." suggest a revised halo mass growth.
-
 
 
 == Effect on radiation profiles
 #let notebook = json("../workdir/11_visualization/alpha_dependence_of_profiles.ipynb")
 
-In order to illustrate the necessity of a more precise treatment of the halo mass history, we first investigate the effect of different mass accretion rates on the resulting radiation profiles. To this end we consider halos at fixed masses and vary their accretion rates around the fiducial value of $alpha = 0.79$.
+To illustrate the necessity of a more precise treatment of the halo mass history, we first investigate the effect of different mass accretion rates on the resulting radiation profiles. To this end, we consider halos at fixed masses and vary their accretion rates around the fiducial value of $alpha = 0.79$.
 
 #figure(
   image_cell(notebook, cell_id: "profile_plot_alpha_dependence"),
   caption: [
     Flux profiles around halos with varying accretion rates.
-    _Left:_ Profile of the Lyman-$alpha$ coupling coefficient.
-    _Center:_ Profile of the kinetic temperature $T_k$.
-    _Right:_ Ionization fraction profile.
+    _Left:_ Profile of the Lyman-$alpha$ coupling coefficient $rho_alpha$.
+    _Center:_ Profile of the heating profile $rho_h$.
+    _Right:_ Ionization fraction profile $x_"HII"$.
     The effect of different mass accretion rates is visualized by the color gradient where bluer colors correspond to lower accretion rates and redder colors to higher accretion rates.]
 
 ) <fig:profile_plot_alpha_dependence>
 
-@fig:profile_plot_alpha_dependence shows the three relevant profiles, computed for $M_"h1" = ??$ and $M_"h2" = ??$. The variation of the accretion rate leads to noticeable differences in all three profiles, even at high radial distances. There is a clear and consistent trend for all three profiles: Higher accretion rates lead to higher fluxes, i.e. an effect that is more outreaching. This is expected as a higher accretion rate leads to a higher star formation rate and thus to the production of more photons.
+@fig:profile_plot_alpha_dependence shows the three relevant profiles, computed for $M_"h" = 6.08 dot 10^11 M_dot.circle$. The variation of the accretion rate leads to noticeable differences in all three profiles, even at high radial distances. There is a clear and consistent trend for all three profiles: Higher accretion rates lead to higher fluxes, i.e. an effect that is more outreaching. This is expected as a higher accretion rate leads to a higher star formation rate and thus to the production of more photons. Note that the spread of the profiles seems symmetric. This is due to the logarithmic scaling of the plot and the higher range of alpha values above the fiducial value. In other words, the magnitude of the effect seems to be the same when shifting the accretion above or below the fiducial value.
 
 // TODO - how far should I comment on that?
-This picture is more complex once we consider a distribution of accretion rates instead of a single value.
-// Do I need to show a plot of that as well?
-We note that the dominating factor when considering a distribution is the contribution from the mean accretion rate. The scatter around the mean has a significantly smaller effect. We do not elaborate on the stochasticity of the accretion rate since the usage of #nbody simulations allows for a more sophisticated investigation. Instead of assuming pure stochasticity we can extract the actual growth history of each halo and use it to assign a more meaningful accretion rate.
+This picture is more complex once we consider a distribution of accretion rates instead of a single value across all halos. We note that the dominating factor when considering a distribution is the contribution from the mean accretion rate. The scatter around the mean has a significantly smaller effect. We do not elaborate on the stochasticity of the accretion rate since the usage of #nbody simulations allows for a more sophisticated investigation. Instead of assuming pure stochasticity we can extract the actual growth history of each halo and use it to assign a more meaningful accretion rate.
 
 
 
 == Merger trees
 
-=== Using THESAN
+=== Using #thesan simulation data
 
 In order to generate precise map-level predictions of the 21-cm signal, #beorn combines the halo model of reionization with large-scale #nbody simulations which provide realistic snapshots of the dark matter distribution. They give a spatial context to the generated profiles.
 
 As described in @procedure #beorn was initially used to postprocess the #pkdgrav
-// cite!
-simulation suite and to obtain a meaningful signal capable of constraining astrophysical parameters related to star formation. The aim of this thesis is not to merely increase the precision but to leverage the mass history that can be extracted directly from the simulation to refine the underlying model.
+#cite(<potter2016pkdgrav3trillionparticlecosmological>, form: "normal")
+simulation suite and to obtain a meaningful signal capable of constraining astrophysical parameters related to star formation. The aim of this thesis is not to merely increase the precision but to develop a proof of concept that leverages the mass history which can be extracted directly from the simulation to refine the underlying model.
 
 To this end, we use the publicly available data from the #thesan simulation suite
 #cite(<Kannan_2021>, form: "normal")
 #cite(<Garaldi_2022>, form: "normal")
 #cite(<Smith_2022>, form: "normal")
-. The #thesandark simulation in particular provides a dark-matter-only simulation and already provides halo catalogs and merger trees generated by the `LHaloTree` tree builder by @Springel2005. This will allow us to extract the growth of each halo accross different snapshots without signifcant preprocessing.
+. The #thesandark simulation in particular provides a dark-matter-only simulation that conveniently already includes halo catalogs and merger trees generated by the `LHaloTree` tree builder by @Springel2005. This will allow us to extract the growth of each halo across different snapshots without significant preprocessing.
 
 With a box length of $95.5 "cMpc"$ the simulation provides a sufficient volume to avoid box size effects
-// CITATION
+(e.g. @Iliev_2014)
 while still allowing us to iterate quickly and test the refined model without excessive computational cost. The simulation has two variants with different mass resolutions:
-#thesandark 1 with $2100^3$ particles for a mass resolution of $3.70 dot 10^6 M_dot.circle$ per particle and #thesandark 2 with $1050^3$ particles for a mass resolution of $2.96 dot 10^7 M_dot.circle$ per particle. Unless specified otherwise we use #thesandark 2 since it provides a good compromise between resolution and computational cost. We make use of #thesandark 1 to perform convergence tests as described in @validation.
-
-
-
-// TODO - below
-// @Kannan_2021 also shows that reionization history is different for different gas densitites, i.e. halo masses. We also show from a profile perspective that treating halo accretion as a free parameter can lead to significant differences in the resulting profiles.
-
-
-// Thesan halo catalog and the motivation to increase the cutoff.
-
-// At the same time THESAN low mass halos seem overabundant which is why we use boosted models of star formation efficiency.
+#thesandark 1 with $2100^3$ particles for a mass resolution of $3.70 dot 10^6 M_dot.circle$ per particle and #thesandark 2 with $1050^3$ particles for a mass resolution of $2.96 dot 10^7 M_dot.circle$ per particle. Unless specified otherwise we use #thesandark 2 since it provides a good compromise between good resolution and computational cost. We make use of #thesandark 1 to perform convergence tests as described in @validation.
 
 
 
@@ -112,31 +78,31 @@ while still allowing us to iterate quickly and test the refined model without ex
 Growth of structure in #lambdacdm is hierarchical: Small structures form first and merge to form larger structures. The growth of halos can be represented using merger trees. These tree-like structures describe the halo history in terms of the mergers of its smaller progenitors. A merger tree is constructed by linking halos in consecutive snapshots of the simulation where each halo as a single descendant but potentially multiple progenitors.
 // As described in ... THESAN
 
-The main progenitor serves as a tracer of the halo mass history if we assume that the halo mass growth is dominated by mergers.
+The main progenitor is the most massive progenitor and serves as a tracer of the halo mass history if we assume that the halo mass growth is dominated by mergers.
 // Has this been explicitlyshown somewhere?
 Beyond that, we expect the main progenitor to be most representative of the baryonic conditions inside and outside the halo as the merger occurs.
 // Might need to reformulate
 For the identification of accretion rates for #beorn we therefore focus solely on the main progenitor branch of each halo.
 
-Reducing the breadth of the merger tree reduces the data volume significantly and allows us to implement the tree handling in memory without excessive computational cost. To this end, we provide a simple implementation of a tree walker that copies the simplified trees to a single file for easier access. Other preprocessing is not required which allows #beorn to keep all parameters related to the mass history as free parameters to be specified at runtime.
+Reducing the breadth of the merger tree reduces the data volume significantly and allows us to implement the tree handling in memory without excessive computational cost. For this purpose we provide a simple implementation of a tree walker that copies the simplified trees to a single file for easier access. Other preprocessing is not required, which allows #beorn to keep all parameters related to the mass history as free parameters to be specified at runtime.
 
 
 === Fitting procedure
-The restriction to the main progenitor corresponds to a reduction the dimensionality of the mass history to a one-dimensional function of redshift compatible with the orginal assumption of an exponential growth model as in @eq:exponential_growth.
+The restriction to the main progenitor corresponds to a reduction of the dimensionality of the mass history to a one-dimensional function of redshift compatible with the orginal assumption of an exponential growth model as in @eq:exponential_growth.
 
 #figure(
   image_cell(notebook, cell_id: "merger_tree_and_fitting"),
   caption: [
     Usage of merger tree fitting to obtain accretion rate estimates.
     _Left:_ Collection of normalized main progenitor branches with mass $M_"mp"$ starting at $z = 10.3$ and looking back over $n=10$ snapshots. Select histories and their corresponding exponential fits are highlighted.
-    _Right:_ Distribution of best-fit accretion rates $alpha$ for all halos at $z=10.3$.
+    _Right:_ Distribution of best-fit accretion rates $alpha$ for all halos at $z=8.29$.
   ]
 
 ) <fig:merger_tree_and_fitting>
 
-We use a linear regression in log-space to obtain estimates of the accretion rate $alpha$ for each halo. This is implemented in a vectorized fashion to allow for efficient processing of the full dataset. For this fit we enforce the current halo mass as a boundary condition. This prevents inconsistent fits where the latest fitted mass deviates from the actual current halo mass. As a visualization of the fitting procedure @fig:merger_tree_and_fitting shows a collection of normalized main progenitor branches starting at $z=10.3$ and looking back over $n=10$ snapshots. After fitting we overlay the estimated exponential growth history for a selection of halos. The right panel shows the distribution of best-fit accretion rates $alpha$ for all halos at $z=10.3$.
+We use a linear regression in log-space to obtain estimates of the accretion rate $alpha$ for each halo. This is implemented in a vectorized fashion to allow for efficient processing of the full dataset. For this fit we enforce the current halo mass as a boundary condition. This prevents inconsistent fits where the latest fitted mass deviates from the actual current halo mass. As a visualization of the fitting procedure @fig:merger_tree_and_fitting shows a collection of normalized main progenitor branches starting at $z=8.29$ and looking back over $n=10$ snapshots. After fitting we overlay the estimated exponential growth history for a selection of halos. The right panel shows the distribution of best-fit accretion rates $alpha$ for all halos at $z=8.29$. Given the relative low mass of the halos we observe a strong clustering of accretion rates around a value of $alpha approx 0.5$. Outliers with significantly deviating values appear nevertheless and are not linked to a specific mass range.
 
-Similarly to the halo mass itself the accretion rate can then be taken into account during the painting procedure by selecting a profile corresponding to the halo mass and accretion rate of each halo. Hence the accretion rate is binned as well and the range that is covered during the painting is finite. We leave this as a free parameter to be specified at runtime.
+Similarly to the halo mass itself, the accretion rate can then be taken into account during the painting procedure by selecting a profile corresponding to the halo mass and accretion rate of each halo. Consequently, the accretion rate is binned as well and the range that is covered during the painting is finite. We leave this as a free parameter to be specified at runtime.
 
 
 
@@ -146,15 +112,15 @@ Similarly to the halo mass itself the accretion rate can then be taken into acco
 #figure(
   image_cell(notebook, cell_id: "alpha_evolution_vs_redshift"),
   caption: [
-    Evolution of the mean of the fitted accretion rates and the $1 sigma$ standard deviation (shaded area). For a given snapshot we consider different numbers of snapshots $n$.
+    Evolution of the mean of the fitted accretion rates and the $1 sigma$ standard deviation (shaded area). For a given snapshot we consider different numbers of lookback snapshots $n$.
   ]
 ) <fig:alpha_evolution_vs_redshift>
 
-In order to obtain a sensible range of $alpha$ values to cover during the painting procedure, we investigate the global result of the fitting procedure. Our method of fitting trades speed and convenience for absolute precision: Not all halos are well represented in the merger tree. Additionally, we need to account for unphysical or incomplete histories due to limitations of the halo finder. We discuss this step in @implementation. For the current investigation we disregard these halos and only consider well-behaved, fully reolved trees. @fig:alpha_evolution_vs_redshift shows how the fitted accretion rate evolves when starting from the different snapshots. We plot the mean and $1 sigma$ standard deviation of the resulting distribution of $alpha$ values. We consider different lookback lengths with the goal of assessing the stability of the fitting procedure.
+In order to obtain a sensible range of $alpha$ values to cover during the painting procedure, we investigate the global result of the fitting procedure. Our method of fitting trades speed and convenience for absolute precision: Not all halos are well represented in the merger tree and no further processing is done beyond the consideration of the tree. Additionally, we need to account for unphysical or incomplete histories due to limitations of the halo finder. We discuss this step in @implementation. For the current investigation we disregard these halos and only consider well-behaved, fully resolved trees. @fig:alpha_evolution_vs_redshift shows how the fitted accretion rate evolves when starting from the different snapshots. We plot the mean and $1 sigma$ standard deviation of the resulting distribution of $alpha$ values. We consider different lookback lengths with the goal of assessing the stability of the fitting procedure.
 
-We observe a clear stabilization of the mean accretion for longer lookbacks. Not only does it make sense to consider longer lookbacks because of their causal connection, but also because it helps to absorb short-term fluctuations most likely introduced by the halo finder. This is especially noticeable in the first few snapshots where the $1 sigma$ uncertainty is significantly higher. This is likely due to the overabundance of low mass halos whose mass history is more erratic and harder to reconstruct, accentuated by displacements of the halos.
+We observe a clear stabilization of the mean accretion for longer lookbacks. Not only does it make sense to consider longer lookbacks because of their causal connection, but also because it helps to absorb short-term fluctuations most likely introduced by the halo finder. This is especially noticeable in the first few snapshots where mean and $1 sigma$ uncertainty are significantly higher. This is probably due to the overabundance of low mass halos whose mass history is more erratic and harder to reconstruct.
 
 Numerically, the advantage of longer lookbacks is the stabilization of the fit leading to reduced scatter in the resulting distribution. We note that these behaviors stabilize once we consider around $n = 10$ snapshots of lookback. Both the mean and standard deviation follow a stable trend and the mean settles at $alpha approx 0.6$.
 We attribute the slight offset of the means to the fact that discarding incomplete trees favors more massive halos at higher lookbacks. These halos are more stable in terms of detection by the halo finder and are expected to have fewer fluctuations.
 
-Physically, the lookback time is motivated from the flux profiles of the halos themselves. Due to the size up to the $"Mpc"$ range (see e.g. @fig:profile_plot_alpha_dependence) we attribute to each profile a timescale that causally affects the region defined by that profile. For a profile of radius $ #sym.tilde.op 100 "Mpc"$ this time is on the order of $Delta t = 300 "Myr" #sym.arrow.l.r.double.long Delta z = 4$ (when looking back from a redshift of $z=8$). Given the spacing of snapshots in #thesan using $n=10$ snapshots still lies below the causal range. Since the fitted behavior seems to stabilize we suggest to not go beyond that since the consideration of additional snapshots slows down the simulation considerably.
+Physically, the lookback time is motivated from the flux profiles of the halos themselves. Due to their size in the $"cMpc"$ range (see e.g. @fig:profile_plot_alpha_dependence) we attribute to each profile a timescale during which there is a causal effect on the region defined by the extent of that profile. For a profile of radius $#sym.tilde.op 100 "cMpc"$ this time is of the order of $Delta t = 300 "Myr"$, corresponding to $Delta z = 4$ (when looking back from a redshift of $z=8$). Given the spacing of snapshots in #thesan using $n=10$ snapshots still lies below the causal range. Since the fitted behavior seems to stabilize we suggest to not go beyond that since the consideration of additional snapshots slows down the simulation considerably.

implementation.typ

@@ -9,34 +9,34 @@ For each halo we require a flux profile that matches the halo properties which n
 // Maybe reformulate
 Since the dynamic range of accretion rates is large, the resulting parameter space rapidly expands. The computation of the profiles therefore utilizes vectorized operations to achieve reasonable runtimes.
 
-// TODO - maybe put somewhere else
-// at least explain why it isn't a problem
-// Reformulate
-Note that this introduces another "second degree" inconsistency: The flux profile attributes a radiative behavior to the halo that is motivated by its history. This is repeated for each snapshot creating possibly conflicting histories. In the case of stable halo growth this is not a problem but in the case of erratic growth (e.g. major mergers) this can lead to unphysical behavior. A more consistent approach would be to assume a more flexible mass growth model that distinguishes different growth modes/regimes.
+Note that the precomputation of the profiles introduces another "second degree" inconsistency: The flux profile attributes a radiative behavior to the halo that is motivated by its history. This is repeated for each snapshot creating possibly conflicting histories. In the case of stable halo growth this is not a problem but in the case of erratic growth (e.g. major mergers) this can lead to unphysical behavior. A more consistent approach would be to assume a more flexible mass growth model that distinguishes different regimes of growth. This would require a much more complex handling of the precomputed profiles and is beyond the scope of this work. The current approach remains a good approximation for the majority of halos.
 
 
-== Parallel binned painting
-Similary to the computation of profiles, the painting step is affected by the increased parameter space. #beorn's fast simulation times revolve around the crucial simplification of the halo model: Halos with the same core properties are treated identically and can be mapped onto the grid in a single operation. Through the addition of the accretion rate as a parameter the degeneracy is reduced. The number of halos that can be treated simultaneously decreases, even though their mass is identical. To mitigate this effect we implement a parallelized version of the painting step that distributes the workload to multiple processes
+== Parallel painting of profile bins
+Similarly to the computation of profiles, the painting step is affected by the increased parameter space. #beorn's fast simulation times revolve around the crucial simplification of the halo model: Halos with the same core properties are treated identically and can be mapped onto the grid in a single operation (see @painting). Through the addition of the accretion rate as a parameter the degeneracy of identical halos is reduced. The number of halos that can be treated simultaneously decreases even though they have the same mass. To mitigate this effect we implement a parallelized version of the painting step that distributes the workload to multiple processes
 #footnote[
-  A rudimentary parallel implementation using `MPI` already exists. It leverages the fact that each snapshot can processed independently and distributes the snapshots to multiple processes.
+  A rudimentary parallel implementation using the message passing interface (`MPI`) already exists. It leverages the fact that each snapshot can be processed independently and distributes the snapshots to multiple processes.
+].
+This implementation utilizes a shared memory approach and uses processes on a single node that share a common memory space to store the grid. This allows for a more efficient usage of node resources since the memory overhead of duplicating the grid for each process is avoided. The required pre- and postprocessing that ensure the correct execution of all processes are justified by the performance gain which is nearly linear with the number of processes used
+#footnote[
+  We test the scaling with a parallelization up to 70 processes and observe a continuous speedup. Part of this speedup is absorbed by the overhead of the much larger number of bins.
 ].
-This implementation utilizes a shared memory approach and uses processes on a single node that share a common memory space to store the grid. This allows for a more efficient usage of node resources since the memory overhead of duplicating the grid for each process is avoided.
 
-Part of the painting procedure remains inherently sequential: The final ionization map requires conservation of the total photon count. This is achieved by distributing duplicate ionizations to neighboring cells.
-// Reformulate
-a parallel approach cannot guarantee perfect consistency. We aim to keep the single process computations to a minimum.
+Part of the painting procedure remains inherently sequential. For instance the final ionization map requires conservation of the total photon count. This is achieved by distributing duplicate ionizations from overlapping bubbles to neighboring cells. A parallel spreading approach might create new overlaps that would require further iterations to resolve. We therefore perform this step in a single process. We aim to keep these inefficient computations to a minimum.
+
 
 == Merger tree processing
 The central improvement of the simulation procedure is the consideration of the individual halo mass accretion histories during the painting and not just the assumption of a predefined value. As described in @halo_mass_history we utilize the merger trees provided by the #thesan simulation. The inference of the accretion rate is performed at runtime. Further preprocessing of the simulation is not required, only a single step that merges the individual tree files into a single file.
 
-The generated alphas are binned as a result of the painting procedure and the permitted range is restricted as specified in the configuration. For our runs we find that an upper limit of $alpha = 5$ only affects a sub-percent fraction of halos. Many of these halos exhibit erratic growth suggesting that allowing for very high accretion rates is not physical.
+The generated $alpha$ values are binned as a result of the painting procedure and the permitted range is restricted as specified in the configuration. For our runs we find that an upper limit of $alpha = 5$ only affects a sub-percent fraction of halos. Many of these halos exhibit erratic growth suggesting that allowing for very high accretion rates is not physical.
 
-The #thesan data provides a convenient way to iterate and refine the above procedure but is not without shortcomings. The merger trees are constructed in post-processing and do not guarantee self-consistency of halo properties accross multiple snapshots. This manifests itself through negative growth rates that cannot be represented in the current model. Furhtermore the mass resolution of the #thesandark simulations is apparently too coarse to accurately resolve halos down to the atomic cooling limit of $M_"h" = 10^8 M_dot.circle$. This is an issue that becomes apparent in @validation where we compare the impact of the different mass resolutions. To account for this we follow the description of star formation efficiency employed by @Schaeffer_2023 picking a "boosted" model for the description of our halos. The resulting parameters for @eq:star_formation_efficiency are $f_(star,0) = 0.1$, $M_p = 2.8 times 10^(10) M_dot.circle$, $g_1 = 0.49$ and $g_2 = -0.61$.
+The #thesan data provides a convenient way to iterate and refine the above procedure but is not without shortcomings. The merger trees are constructed in postprocessing and do not guarantee self-consistency of halo properties across multiple snapshots. This manifests itself through negative growth rates that cannot be represented in the current model. Furthermore the mass resolution of the #thesandark simulations is apparently too coarse to accurately resolve halos down to the atomic cooling limit of $M_"h" = 10^8 M_dot.circle$. This is an issue that becomes apparent in @validation where we compare the impact of the different mass resolutions. To account for this we follow the description of star formation efficiency employed by @Schaeffer_2023, picking a "boosted" model for the description of our halos. The resulting parameters for @eq:star_formation_efficiency are $f_(star,0) = 0.1$, $M_p = 2.8 times 10^(10) M_dot.circle$, $gamma_1 = 0.49$ and $gamma_2 = -0.61$.
 
 
 == Secondary changes
-Additionally to the changes directly linked to the new accretion model we implement several improvements that allow for better usability and reproducability of the simulation outputs.
+In addition to the changes directly linked to the new accretion model we implement several improvements that lead to better usability and reproducibility of the simulation outputs.
 
 We improve the input/output handling by implementing proper `HDF5` support and caching of intermediate results. This allows for a more efficient usage of disk space and faster loading times. It also enables the resumption of interrupted simulations.
-The import of data from the original #nbody simulation has been generalized to a reference class to ensure modularity and allow for easier extension to other simulations. This has been part of a larger overhaul of the codebase to improve modularity and readability.
-A general speedup from the cumulated effect of the above changes and code optimizations allows for a faster painting procedure. A contribution to that speedup comes from the ussage of `Pylians` by @Pylians. It provides efficient implementations in `C` of of the grid mapping of the individual particles. This additionally allows for a rigorous implementation of redshift space distortions (RSD) by utilizing the exact velocity information of each dark matter particle individually. Previous implementations of RSD in #beorn were based on approximations of the velocity field derived from the density field. The impact of RSD on the 21-cm signal has been discussed e.g. by @Ross_2021 but is not the focus of this work.
+The import of data from the original #nbody simulation has been generalized to a reference class to ensure modularity and easy adaptation to other simulations. This has been part of a larger overhaul of the codebase to improve modularity and readability. #beorn aims to be a flexible framework that produces fast results that the end user can customize to reflect their parameter choices. Usability is therefore a key aspect of the code design.
+
+A general speedup from the cumulated effect of the above changes and code optimizations results in a faster painting procedure. A contribution to that speedup comes from the usage of the `pylians` package by @Pylians. It provides efficient implementations of the grid mapping of the individual particles. This additionally allows for a rigorous implementation of redshift space distortions (RSD) by utilizing the exact velocity information of each dark matter particle individually. Previous implementations of RSD in #beorn were based on approximations of the velocity field derived from the density field. The impact of RSD on the 21-cm signal has been discussed e.g. by @Ross_2021 but is not the focus of this work.

introduction.typ

@@ -1,65 +1,64 @@
 #import "helpers.typ": *
 = Introduction
 
-The earliest cosmological events (such as the formation of the first astrophysical objects - stars, galaxies, black holes...) have a profound influence on the evolution of the universe. Though poorly understood, these events have shaped the characteristics of our current universe, including the structure and distribution of matter itself.
+The formation of the earliest astrophysical objects, such as stars, galaxies, and black holes, have a profound influence on the evolution of the universe. Though poorly understood, these events have shaped the characteristics of our current universe, including the structure and distribution of matter itself.
 // Citation about an overview paper on ionization vs structure formation.
 Despite the milestones achieved in observational cosmology, many aspects of the early universe and its dark ages remain difficult to probe. While traditional measurements provide insights into relatively recent epochs, and the cosmic microwave background (CMB) serves as an early snapshot of the universe, the dark ages are incompatible with direct observations. They represent a critical link between the late-time universe and the primordial conditions that has remained largely unexplored.
 
 The epoch of reionization (EOR) spans the time period from the end of the dark ages until the time when the universe is fully ionized again. It is a period of complex interactions between matter and radiation but it is crucial to understand as it sets the stage for the subsequent evolution of the universe.
 // including the formation of galaxies and large-scale structures.
 // It simultaneously is affected by the fundamental mechanisms and also affects the subsequent evolution of the universe.
-Beyond its impact on the late universe, a detailed understanding of the reionization process has been shown to provide new and competitive constraints on the current cosmological model (e.g
+Beyond its impact on the late universe, a detailed understanding of the reionization process has been shown to provide new and competitive constraints on the current cosmological model (e.g.
 @Mao_2008
 @McQuinn_2006
 @schneider2023cosmologicalforecast21cmpower
 ).
-Understanding and being able to model the EOR is therefore crucial for a comprehensive picture of cosmology.
+Understanding and being able to model the EOR is therefore essential for a comprehensive picture of cosmology.
 
-The dark ages of the universe refer to the period after recombination where the primordial atoms remained neutral. They are characterized by the total lack of sources of radiation (beyond the radiation background). The dominant interactions during that period were either gravitational or due to the cooling of the primordial gas. The formation of the first stars was obstructed by the lack of efficient cooling mechanisms in the absence of heavier nuclei. With the simplest cooling channel being the deexcitation of atomic hydrogen, the gas inside a virialized structure can only collapse if the enclosed mass is high enough. This so called atomic cooling limit sets a minimum mass for the halos that can host star formation at around $10^8 M_(dot.circle)$. Other cooling channels such as the deexcitation of molecular hydrogen were suppressed by the emission of photons from the first stars.
+The dark ages of the universe refer to the period after recombination where the primordial atoms remain neutral. They are characterized by the total lack of sources of radiation (beyond the radiation background). The dominant interactions during that period are either gravitational or due to the cooling of the primordial gas. During that time the formation of the first stars is obstructed by the lack of efficient cooling mechanisms in the absence of heavier nuclei. With the simplest cooling channel being the deexcitation of atomic hydrogen, the gas inside a virialized structure can only collapse if the enclosed mass is high enough. This so-called atomic cooling limit sets a minimum mass for the halos that can host star formation at around $10^8 M_(dot.circle)$. Other cooling channels such as the deexcitation of molecular hydrogen are suppressed by the emission of photons from the first stars.
 // => argument that there is no "galaxy" in that sense below
 
-The first stars mark the end of the dark ages. These so called population III stars were metal-free and their short lifespan ended in supernovae that enriched and heated the surrounding gas in the intergalactic medium (IGM).
+The formation of the first stars marks the end of the dark ages. These so-called population III stars have zero metallicity and very distinct characteristics compared to later generations of stars. Their existence has not been confirmed observationally but they are thought to have shaped the subsequent formation of stars and galaxies and to have played a crucial role in the reionization of the universe (e.g. @Mebane_2020).
 
-// Citation about Pop III stars and their role in the cosmic dawn.
-which...
-During the cosmic dawn ...
+// Paragraph talking about the evolution of the IGM and the formation of ionized bubbles around sources.
+Driven by the newly formed stars and galaxies, reionization is explained as an inside-out process (e.g. @10.1111-j.1365-2966.2006.10502.x) expanding from within the halos that host the first galaxies. The ionizing radiation emitted by these sources reaches the intergalactic medium (IGM) and creates ionized bubbles that grow and eventually overlap to fully ionize the universe again.
 
+While reionization marks the gradual disappearance of neutral hydrogen, the preceeding abundance during the dark ages and cosmic dawn allows for an additional mode of observation: the 21-cm line. Due to the hyperfine transition of neutral hydrogen there is a characteristic emission or absorption of photons at a frequency of $1420 "MHz"$. The strength of this signal depends on the local conditions, encoded by the spin temperature. The redshifting of the photons allows to probe different epochs through the observed frequency.
 
+The detection of the 21-cm signal of reionization is a major goal of current and upcoming radio telescopes, for instance the
+Square Kilometer Array #cite(<SKAlow>, form: "normal", supplement: "SKA")
+or the Hydrogen Epoch of Reionization Array #cite(<HERA>, form: "normal", supplement: "HERA")
+. These instruments are expected to detect the power spectrum of the 21-cm signal, providing further insights into the dynamics of the early universe. In particular the low-frequency component SKA-Low is expected to have the sensitivity to image the 21-cm signal directly and to produce maps of the ionization field during the EOR.
 
-The large amounts of neutral hydrogen in the intergalactic medium during the dark ages and cosmic dawn allow for an additional mode of observation: the 21-cm line emission. Due to the hyperfine transition of neutral hydrogen there is a characteristic emission or absorption of photons at a frequency of $1420 "MHz"$. The strength of this signal depends on the local conditions, in particular the redshifting of the photons allows to probe different epochs through the observed frequency.
+Beyond observations, an additional pillar of understanding the EOR is the modeling and simulation of the universe during that time. The main purpose of simulations is to constrain the EOR observables. Combined with the first observations, simulations will generate a wealth of information about the early universe at a range of redshifts that has previously been inaccessible. With the highest sensitivity and resolution forecasted for these observations, the simulations must be able to capture the full dynamic range of the interactions from the small-scale physics of star formation and feedback to the large-scale structure of the universe.
 
-The main purpose of simulations is to constrain EOR observables, in particular the 21-cm signal.
-// Keep the below?
-Combined with upcoming observations from ... these simulations will generate a wealth of information about the early universe, at a range of redshifts that has previously been inaccessible. With the highest sensitivity and resolution forecasted for these observations, the simulations must be able to capture the full dynamic range of the interactions, from the small-scale physics of star formation and feedback to the large-scale structure of the universe.
-
-State of the art simulations need to implement a range of physical processes, including gravitational interactions, hydrodynamics, radiative transfer, and feedback mechanisms. Prominent examples include the THESAN simulations
+State-of-the-art simulations need to implement a range of physical processes, including gravitational interactions, hydrodynamics, radiative transfer, and feedback mechanisms. Prominent examples include for instance the #thesan simulations
 #cite(<Kannan_2021>, form: "normal")
 #cite(<Garaldi_2022>, form: "normal")
 #cite(<Smith_2022>, form: "normal")
-and ... .
-Another approach is to use ray-tracing algorithms which give detailed descriptions of the radiative transfer.
-// C2ray?
-These methods are computationally expensive, which limits their applicability for large-scale simulations.
+.
+Another approach is to use ray-tracing algorithms which give detailed descriptions of the radiative transfer (see e.g. @MELLEMA2006374).
+// Shortcomings of similar codes => justification for the development of #beorn (@Schaeffer_2023).
+These simulations are computationally expensive and cannot be used to to repeatedly explore the parameter space of reionization.
 
-// Shortcomings of similar codes (as noted in #beorn paper). => justification for the development of #beorn (@Schaeffer_2023).
 
-This work presents #beorn, the _Bubbles during the Epoch of Reionization Numerical simulator_ by @Schaeffer_2023, and the refinements we make to achieve self-consistency.
-// not clear!
-In its simplest description #beorn is the implementation of the "halo model of reionization" by @schneider2023cosmologicalforecast21cmpower. In this model the radiative interactions are treated as spherically symmetric around a halo-scale source. This effectively reduces the dimensionality of the radiative transfer problem. #beorn uses the one-dimensional (1-d) profiles generated by this model to paint the 3-d space around sources which are obtained from a large-scale #nbody simulation. A distinguishing feature of #beorn is the self-consistent treatment of the growth of individual sources over the course of the simulation. The first iteration of #beorn focused on the effect of emitted photons whereas this work focuses on the effects of gravitational mass accretion. We show that the radiation profiles are sensitive to the growth rate of the sources and that an accurate treatment of the source growth has an impact on the resulting 21-cm signal.
-// Mention that this is treated in more detail in @procedure
+This work presents #beorn, the _Bubbles during the Epoch of Reionization Numerical simulator_ by @Schaeffer_2023. In its simplest description #beorn is the implementation of the "halo model of reionization" by @schneider2023cosmologicalforecast21cmpower. In this model the radiative interactions are described through spherically symmetric profiles around sources embedded in dark matter halos. This effectively reduces the dimensionality of the radiative transfer problem. #beorn uses the one-dimensional (1-d) profiles generated by this model to paint the 3-d space around halos which are obtained from a large-scale #nbody simulation. A distinguishing feature of #beorn is the self-consistent treatment of the expansion of the affected regions around the sources. This approach allows for simulations at the largest scales while still taking into account the core processes of reionization. The computational efficiency of #beorn makes it suitable to explore the underlying parameters.
 
-This report is structured as follows: @procedure describes the details of the simulation procedure, including the underlying model. @halo_mass_history explains how mass evolution is modeled and its impact on the profiles.
-// not any profiles.
-In @implementation we give an overview of the implementation of the modeling assumed by #beorn and the steps required to produce a full 3-d lightcone simulation.
+// growth of individual sources over the course of the simulation.
+The first iteration of #beorn focused on the impact of parameters related to the emission of photons whereas this work focuses on the effects of gravitational mass accretion. We show that the radiation profiles are sensitive to the growth rate and that the mass accretion history provdied by #nbody simulations is too complex to be captured by simple parametrizations. Our improved version of #beorn permits a more consistent treatment by considering the individual mass accretion history of each source. We demonstrate the resulting measurable effects on the 21-cm observables when compared to the simpler models.
+
+This report is structured as follows: @procedure describes the details of the simulation procedure, including the underlying model. @halo_mass_history explains how mass evolution is modeled and its impact on the flux profiles used by #beorn.
+In @implementation we give an overview of the implementation and changes required by the refined modeling.
+// assumed by #beorn and the steps required to produce a full 3-d lightcone simulation.
 // self-consistent treatment of mass accretion.
 @validation details the validation we perform on the refined procedure and in @results we compare the resulting signals to quantify the impact of different models of mass accretion. @conclusion summarizes our findings and discusses potential future improvements.
 
+Note that throughout this report, physical distance units are specified with the prefix "p", while comoving distance units are specified with the prefix "c".
 
 
+// Other points to mention
+// - wouthuysen
+// - cold reionization
+// - comoving distances - check consistency
 
-Other points to mention
-- wouthuysen
-- cold reionization
-- comoving distances - check consistency
-
-how #beorn compares to traditional approaches
+// how #beorn compares to traditional approaches

main.typ

@@ -26,7 +26,7 @@
   title: "Simulating the EOR with self-consistent growth of galaxies",
   subtitle: "Master's Thesis",
   authors: ("Rémy Moll",),
-  supervisors: ("Prof. Aurel Schneider",),
+  supervisors: ("Prof. Aurel Schneider", "Prof. Alexandre Refregier"),
   affiliation: "ETH Zürich, Universität Zürich",
   abstract: include("abstract.typ"),
   background-color: color.rgb(32, 64, 123),
@@ -54,4 +54,5 @@
 
 #bibliography("references.bib", style: "assets/the-astrophysical-journal.csl")
 
+#pagebreak()
 #include "appendix.typ"

procedure.typ

@@ -11,11 +11,11 @@ This section presents the model describing the sources of radiation that drive r
 
 == The halo mass model of reionization <hmreio>
 
-The central feature
+The distinguishing feature
 // don't like that word
-of #beorn is the parametrization of sources of radiation through the properties of their host dark matter halos. This approach is based on the model presented by @schneider2023cosmologicalforecast21cmpower and gives a description
+of #beorn is the parametrization of sources of radiation through the properties of their host dark matter halo. This approach is based on the model presented by @schneider2023cosmologicalforecast21cmpower and gives a description
 // bad word
-of the 21-cm signal through the treatment of flux profiles around sources. Using these profiles and allowing them to overlap enables efficient computations of the ionization state of the IGM without the need for detailed radiative transfer simulations.
+of the 21-cm signal through the treatment of flux profiles around sources. Using these profiles and allowing them to overlap enables efficient computations of the ionization state of the IGM without the need for detailed radiative transfer simulations. The following description follows the derivation made by @schneider2023cosmologicalforecast21cmpower @Schneider_2021 and we refer to these works for components that are not defined here.
 
 The model describes the emission of ionizing radiation by galaxies. It assumes that the sources are hosted by dark matter halos and expresses the star formation and radiation properties as a function of the halo mass $M_"h"$ and mass accretion rate $dot(M_"h")$. The modeling
 // maybe treatment
@@ -30,70 +30,77 @@ Motivated by abundance matching, @schneider2023cosmologicalforecast21cmpower use
 $
   f_star (M_"h") = f_(star,0) dot (2 (Omega_b / Omega_m)) / ((M_"h"/M_"p")^(gamma_1) + (M_"h"/M_"p")^(gamma_2)) dot S(M_"h")
 $ <eq:star_formation_efficiency>
-where $M_"p"$ is the pivot mass where the efficiency peaks, $gamma_1$ and $gamma_2$ are the low and high mass slopes, and $f_(star,0)$ is the normalization chosen at approximately $0.1$. An additional suppression factor $S(M_"h")$ is introduced to account for reduced star formation in low mass halos, its effect is discussed by @Schaeffer_2023.
+with $M_"p"$ the pivot mass where the efficiency peaks, and $gamma_1$, $gamma_2$ are the low and high mass slopes, and $f_(star,0)$ is the normalization chosen at a value of $f_(star,0) #sym.tilde 0.1$. An additional suppression factor $S(M_"h")$ is introduced to account for reduced star formation in low mass halos whose effect is discussed by @Schaeffer_2023.
 
-// TODO - make clear that this follows @Schneider
 
 === Expression of the profiles
 Derived from the star formation rate the halo model predicts the production and distribution of photons in three distinct energy bands:
-Lyman-$alpha$ photons, ionizing UV photons, and X-ray photons. Each of these bands has a different effect on the IGM and is treated separately.
+Lyman-$alpha$ photons, x-ray photons, and ionizing UV photons. Each of these bands has a different effect on the IGM and is treated separately.
 // Not really sure that's true
-
+==== Lyman-$alpha$ flux profile
 Lyman-$alpha$ photons induce a coupling between the spin temperature and the kinetic temperature of the gas. This effect, known as the Wouthuysen-Field effect
-// TODO - check
 #cite(<Wouthuysen>, form: "normal")
 #cite(<Field>, form: "normal"),
-causes absorption of 21-cm photons before reionization. This is reflected in the absorption expected in the global signal before reionization.
+causes a characteristic absorption and emission feature in the 21-cm signal. Before the cosmic dawn the gas temperature is lower than the background temperature and the coupling leads to absorption. As the first stars heat the gas, the signal transitions to emission.
+
+The Lyman-$alpha$ flux profile around a halo of mass $M$ at redshift $z$ is given by
 $
   rho_alpha (r bar M, z) = (1 + z)^2 / (4 pi r^2) dot sum_(n=2)^(n_m)f_n dot epsilon_alpha (nu prime) dot f_star dot dot(M)(z prime bar M, z)
 $
+which depends on the lookback redshift $z prime$ at which the photon was emitted, hence the expression in terms of $nu prime = nu dot (1+z prime) slash (1+z)$.
 
+Finally the position dependent coefficient is expressed as
+$
+  x_alpha (r bar M, z) = (1.81 dot 10^11) / (1 + z) dot S_alpha (z) dot rho_alpha (r bar M, z)
+$
+using a suppression factor $S_alpha (z)$.
 
-
-The temperature around the sources is described
-// bad word
-by the cooling temperature of the adiabatically expanding universe and the heating due to X-ray photons emitted by the newly formed stars. The temperature profile follows
-
+==== Temperature profile
+The temperature of the IGM around the sources has a strong impact on the 21-cm signal. It is governed by two effects: the cooling temperature of the adiabatically expanding universe and the heating due to X-ray photons emitted by the newly formed stars. The heating profile $rho_h$ follows
 $
   3/2 dot derivative(rho_h (r bar M, z), z) = (3 rho_h (r bar M, z)) / (1 + z) - (rho_"xray" (r bar M, z)) /(k_B (1 + z) H(z))
 $
-which is based on ????
+which is based on the flux profile of x-ray photons $rho_"xray" (r bar M, z)$. The Boltzmann constant is given by $k_B$ and $H(z)$ is the Hubble parameter at redshift $z$.
 
 
+==== Reionization profile
 
-Ionizing photons, i.e. photons with energies above $13.6 "eV"$ experience a large optical depth which justifies the expression
-
-
+The comoving ionized volume around a source of ionizing photons satisfies the differential equation
 $
   derivative(V, t) = dot(N)_"ion"(t) / overline(n)_H^0 - alpha_B dot C / a^3 dot overline(n)_H^0 dot V
 $
+where $alpha_B$ is the recombination coefficient, $C$ is the clumping factor, $a$ is the scale factor, and $overline(n)_H^0$ is the mean density of hydrogen. We expressed this volume in terms of $dot(N)_"ion"$ the rate of change of the total number of ionizing photons per baryon. This description is not dependent on the frequency since we simply consider the contribution of all photons above the ionization threshold.
 
+Ionizing photons, i.e. photons with energies above $13.6 "eV"$, experience a large optical depth, which justifies the sharp cutoff of the ionization profile at the bubble radius $R_b = root(3, 3/ (4pi) V(M,z))$. The radial dependence of the ionized fraction is expressed through the Heaviside step function $theta_"H"$ and reads
 $
-x_("HII")(r bar M, z) = theta_"H" (R_b (M, z) - r) = theta_"H" (root(3, 3/ (4pi) V(M,z)) - r)
+  x_("HII")(r bar M, z) = theta_"H" lr([R_b (M, z) - r], size: #150%)
 $
 
-// introduced inaccuracies
-// e.g. bursty star formation as presented by Romain Teyssier
-
 === Expression of the reionization signal
 
-The observable signal of the 21-cm line is expressed as the differential brightness temperature $d T_"b"$ which describes the contrast to the foreground with the CMB background.
+The above profiles express the local effect of radiation around a single halo as a 1-d simplification. A representation of a typical profile can be seen in @fig:alpha_evolution_vs_redshift. Using an estimate of the spatial distribution of halos, these profiles can be applied to generate a full 3-d map if we assume spherical symmetry.
+
+The observable signal of the 21-cm line is now a spatially dependent quantity obtained from a combination of the mapped quantities. It is expressed as the differential brightness temperature $d T_"b"$ which describes the contrast between the foreground and the background radiation emitted by the CMB.
 // might want to rephrase that
 Following e.g. @Pritchard2012 an expression for $d T_"b"$ is given by
 $
   d T_"b" (bold(x), z) tilde.eq T_0 (z) dot x_"HI" (bold(x), z)  dot (1 + delta_b (bold(x), z)) dot (x_alpha (bold(x), z)) / (1 + x_alpha (bold(x), z) ) dot ((1 - T_"CMB" (z)) / (T_"gas" (bold(x), z)))
 $ <eq:dTb>
 // where $T_0 (z) = 27 "mK" sqrt((1 + z)/10 (0.15 / (Omega_m h^2))) (Omega_b h^2 / 0.023)$, $x_"HI"$ is the neutral hydrogen fraction, $delta_b$ is the baryonic overdensity, $x_alpha$ is the coupling coefficient introduced by the Wouthuysen-Field effect, $T_"CMB"$ is the temperature of the CMB, and $T_"gas"$ is the kinetic temperature of the gas. The expression assumes that the spin temperature is closely coupled to the kinetic temperature which is valid in most regimes of interest during reionization.
-where the background radiation originates from the CMB.
+where the background radiation originates from the CMB. The previous considerations allow us to infer the values of $x_"HI"$, $x_alpha$, and $T_"gas"$ from $x_"HII"$, $rho_alpha$, and $rho_h$ respectively. The baryonic overdensity $delta_b$ is assumed to trace the dark matter overdensity $delta_"dm"$ which is obtained from the underlying #nbody simulation. The amplitude of the signal is given by
+$
+  T_0 (z) = 27 dot (Omega_b h^2) / 0.023 dot sqrt((1 + z)/10 0.15 / (Omega_m h^2)) "mK"
+$
+where $Omega_m$ and $Omega_b$ are the matter and baryonic density parameters.
+
 
 
 == Simulation steps
-The simulation procedure revolves around the implementation of the spherical radiation profiles around halos. We give a brief overview of the main steps here. For a more detailed description of the implementation we refer to @Schaeffer_2023. We discuss our improvements and changes to the original implementation in @implementation.
+The simulation procedure revolves around the implementation of the above spherical radiation profiles around halos. We give a brief overview of the main steps here. For a more detailed description of the implementation we refer to @Schaeffer_2023. We discuss our improvements and changes to the original implementation in @implementation.
 
 === Halo catalog - #nbody simulations
-As a prerequisite, the generation of map data requires a spatial distribution of dark matter halos as well as the underlying density field. Each snapshot can be used to generate a map at the corresponding redshift. #beorn has been successfully validated against mock maps generated by `21cmFAST`
-// TODO cite
-and produced results by using the #pkdgrav #cite(<potter2016pkdgrav3trillionparticlecosmological>, form: "normal") simulation suite.
+As a prerequisite, the generation of map data requires a spatial distribution of dark matter halos as well as the underlying density field. Each snapshot can be used to generate a map at the corresponding redshift. #beorn has been successfully validated against mock maps generated by `21cmFAST` #cite(<21cmfast>, form: "normal")
+and produced meaningful signals by using the #pkdgrav #cite(<potter2016pkdgrav3trillionparticlecosmological>, form: "normal") simulation suite.
 
 === Computation of radiation profiles
 In accordance with the astrophysical parameters set by the user, radiation profiles are computed in order to be applied according to the halo catalog in a subsequent step. The computation of the profiles is deliberately separated from their application onto halos for a more efficient processing. The range of halo masses and redshifts covered by this precomputation is largely determined by the underlying halo catalog since it provides upper bounds on the halo masses
@@ -101,11 +108,15 @@ In accordance with the astrophysical parameters set by the user, radiation profi
   The minimum halo mass that needs to be considered is already constrained by the atomic cooling limit. Depending on the mass resolution of the simulation it might not even be reached.
 ].
 
-=== Painting with the binned approach
-The last step consists of applying the ionization and temperature distributions defined by the profiles onto a 3-d grid. This is done by iterating over the halos in the catalog and using their corresponding profile. For a given profile a 3-d kernel is generated and applied onto the grid via convolution using the `fftw` library. We refer to this procedure as "painting" since the addition of the contributions of each halo allows us to sequentially build up the final map. In general contributions from multiple halos can overlap without any restrictions. The ionization map is treated specially: In order to conserve the overall number of ionizing photons, we ensure that each cell is only ionized once. If multiple halos contribute to the ionization of a cell, the excess photons are redistributed to neighboring cells until they are either used up or reach the edge of the simulation box.
+=== Painting with the binned approach <painting>
+The last step consists of applying the ionization and temperature distributions defined by the profiles onto a 3-d grid. This is done by iterating over the halos in the catalog and using their corresponding profile. For a given profile a 3-d kernel is generated and mapped onto the grid via convolution using the fast fourier transform implemented by the `astropy` #cite(<astropy:2022>, form: "normal") package. We refer to this procedure as "painting" since the addition of the contributions of each halo allows us to sequentially build up the final map. In general, contributions from multiple halos can overlap without any restrictions. The ionization map is treated specially: In order to conserve the overall number of ionizing photons, we ensure that each cell is only ionized once.
+// If multiple halos contribute to the ionization of a cell, the excess photons are redistributed to neighboring cells until they are either used up or reach the edge of the simulation box.
+
+This final step ensures consistent painting: While the contributions to the temperature and coupling maps can be simply added, the ionization map requires a binary treatment. Cells cannot be ionized past a value of $x_"HII" = 1$. Discarding the excesses would violate photon conservation. A redistribution of the excess photons to neighboring cells is performed to ensure that the total number of ionizations matches the total number of emitted photons.
 
 The usage of precomputed profiles is crucial to the efficiency of the simulation but it introduces a discretization in halo mass since each halo is assigned the profile of the closest mass bin. The effect of this simplification has been shown to converge for a sufficient number of mass bins by @Schaeffer_2023.
 
 
+
 === Derivation of global quantities
-The global signal as well as the power spectrum are derived from the map data and compared to other models or observations. Being derived from a full 3-d #nbody simulation, the results are sensitive to the underlying cosmology and the detailed profile modeling. This means that the results depend on the underlying astrophysical model.
+The global signal as well as the power spectrum are derived from the map data to be compared to other models or observations. Being derived from a full 3-d #nbody simulation, the results are sensitive to the underlying cosmology and the detailed profile modeling. This means that the results also depend on the underlying astrophysical model.

references.bib

@@ -74,7 +74,7 @@
 archivePrefix = {arXiv},
        eprint = {astro-ph/0504097},
  primaryClass = {astro-ph},
-       adsurl = {https://ui.adsabs.harvard.edu/abs/2005Natur.435..629S},
+       url = {https://ui.adsabs.harvard.edu/abs/2005Natur.435..629S},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
@@ -265,3 +265,175 @@ archivePrefix = {arXiv},
    year={2021},
    month=jul, pages={3717–3733}
 }
+
+
+# c2ray
+@article{MELLEMA2006374,
+title = {C2-ray: A new method for photon-conserving transport of ionizing radiation},
+journal = {New Astronomy},
+volume = {11},
+number = {5},
+pages = {374-395},
+year = {2006},
+issn = {1384-1076},
+doi = {https://doi.org/10.1016/j.newast.2005.09.004},
+url = {https://www.sciencedirect.com/science/article/pii/S1384107605001405},
+author = {Garrelt Mellema and Ilian T. Iliev and Marcelo A. Alvarez and Paul R. Shapiro},
+keywords = {Cosmology: theory, Galaxies: formation, Galaxies: high-redshift, Intergalactic medium, Radiative transfer, Methods: numerical},
+abstract = {We present a new numerical method for calculating the transfer of ionizing radiation, called C2-ray (conservative, causal ray-tracing method). The transfer of ionizing radiation in diffuse gas presents a special challenge to most numerical methods which involve time- and spatial-differencing. Standard approaches to radiative transport require that grid cells must be small enough to be optically-thin while time steps are small enough that ionization fronts do not cross a cell in a single time step. This quickly becomes prohibitively expensive. We have developed an algorithm which overcomes these limitations and is, therefore, orders of magnitude more efficient. The method is explicitly photon-conserving, so the depletion of ionizing photons by bound-free opacity is guaranteed to equal the photoionizations these photons caused. As a result, grid cells can be large and very optically-thick without loss of accuracy. The method also uses an analytical relaxation solution for the ionization rate equations for each time step which can accommodate time steps which greatly exceed the characteristic ionization and ionization front crossing times. Together, these features make it possible to integrate the equation of transfer along a ray with many fewer cells and time steps than previous methods. For multi-dimensional calculations, the code utilizes short-characteristics ray tracing. The method scales as the product of the number of grid cells and the number of sources. C2-ray is well-suited for coupling radiative transfer to gas and N-body dynamics methods, on both fixed and adaptive grids, without imposing additional limitations on the time step and grid spacing. We present several tests of the code involving propagation of ionization fronts in one and three dimensions, in both homogeneous and inhomogeneous density fields. We compare to analytical solutions for the ionization front position and velocity, some of which we derive here for the first time. As an illustration, we apply C2-ray to simulate cosmic reionization in three-dimensional inhomogeneous cosmological density field. We also apply it to the problem of I-front trapping in a dense clump, using both a fixed and an adaptive grid.}
+}
+
+
+@article{Mebane_2020,
+   title={The effects of population III radiation backgrounds on the cosmological 21-cm signal},
+   volume={493},
+   ISSN={1365-2966},
+   url={http://dx.doi.org/10.1093/mnras/staa280},
+   DOI={10.1093/mnras/staa280},
+   number={1},
+   journal={Monthly Notices of the Royal Astronomical Society},
+   publisher={Oxford University Press (OUP)},
+   author={Mebane, Richard H and Mirocha, Jordan and Furlanetto, Steven R},
+   year={2020},
+   month=feb, pages={1217–1226}
+}
+
+
+@inproceedings{SKAlow,
+   series={AASKA14},
+   title={The Cosmic Dawn and Epoch of Reionisation with SKA},
+   url={http://dx.doi.org/10.22323/1.215.0001},
+   DOI={10.22323/1.215.0001},
+   booktitle={Proceedings of Advancing Astrophysics with the Square Kilometre Array — PoS(AASKA14)},
+   publisher={Sissa Medialab},
+   author={Koopmans, Leon and Pritchard, J and Mellema, G and Aguirre, J and Ahn, K and Barkana, R and van Bemmel, I and Bernardi, G and Bonaldi, A and Briggs, F and de Bruyn, A. G. and Chang, T. C. and Chapman, E and Chen, X and Courty, B and Dayal, P. and Ferrara, A. and Fialkov, A. and Fiore, F and Ichiki, K. and Illiev, I. T. and Inoue, S and Jelic, V and Jones, M and Lazio, J and Maio, U and Majumdar, S and Mack, K. J. and Mesinger, A. and Morales, M F. and Parsons, A. and Pen, U.L. and Santos, M and Schneider, R and Semelin, B and de Souza, R S and Subrahmanyan, R and Takeuchi, T and Vedantham, H and Wagg, J and Webster, R and Wyithe, S and Datta, Kanan Kumar and Trott, C.},
+   year={2015},
+   month=may, collection={AASKA14}
+}
+
+
+@article{HERA,
+   title={Hydrogen Epoch of Reionization Array (HERA)},
+   volume={129},
+   ISSN={1538-3873},
+   url={http://dx.doi.org/10.1088/1538-3873/129/974/045001},
+   DOI={10.1088/1538-3873/129/974/045001},
+   number={974},
+   journal={Publications of the Astronomical Society of the Pacific},
+   publisher={IOP Publishing},
+   author={DeBoer, David R. and Parsons, Aaron R. and Aguirre, James E. and Alexander, Paul and Ali, Zaki S. and Beardsley, Adam P. and Bernardi, Gianni and Bowman, Judd D. and Bradley, Richard F. and Carilli, Chris L. and Cheng, Carina and Acedo, Eloy de Lera and Dillon, Joshua S. and Ewall-Wice, Aaron and Fadana, Gcobisa and Fagnoni, Nicolas and Fritz, Randall and Furlanetto, Steve R. and Glendenning, Brian and Greig, Bradley and Grobbelaar, Jasper and Hazelton, Bryna J. and Hewitt, Jacqueline N. and Hickish, Jack and Jacobs, Daniel C. and Julius, Austin and Kariseb, MacCalvin and Kohn, Saul A. and Lekalake, Telalo and Liu, Adrian and Loots, Anita and MacMahon, David and Malan, Lourence and Malgas, Cresshim and Maree, Matthys and Martinot, Zachary and Mathison, Nathan and Matsetela, Eunice and Mesinger, Andrei and Morales, Miguel F. and Neben, Abraham R. and Patra, Nipanjana and Pieterse, Samantha and Pober, Jonathan C. and Razavi-Ghods, Nima and Ringuette, Jon and Robnett, James and Rosie, Kathryn and Sell, Raddwine and Smith, Craig and Syce, Angelo and Tegmark, Max and Thyagarajan, Nithyanandan and Williams, Peter K. G. and Zheng, Haoxuan},
+   year={2017},
+   month=mar, pages={045001} }
+
+
+
+
+@article{21cmfast,
+    author = {Mesinger, Andrei and Furlanetto, Steven and Cen, Renyue},
+    title = {21cmfast: a fast, seminumerical simulation of the high-redshift 21-cm signal},
+    journal = {Monthly Notices of the Royal Astronomical Society},
+    volume = {411},
+    number = {2},
+    pages = {955-972},
+    year = {2011},
+    month = {02},
+    abstract = {We introduce a powerful seminumeric modelling tool, 21cmfast, designed to efficiently simulate the cosmological 21-cm signal. Our code generates 3D realizations of evolved density, ionization, peculiar velocity and spin temperature fields, which it then combines to compute the 21-cm brightness temperature. Although the physical processes are treated with approximate methods, we compare our results to a state-of-the-art large-scale hydrodynamic simulation, and find good agreement on scales pertinent to the upcoming observations (≳1 Mpc). The power spectra from 21cmfast agree with those generated from the numerical simulation to within 10s of per cent, down to the Nyquist frequency. We show results from a 1-Gpc simulation which tracks the cosmic 21-cm signal down from z= 250, highlighting the various interesting epochs. Depending on the desired resolution, 21cmfast can compute a redshift realization on a single processor in just a few minutes. Our code is fast, efficient, customizable and publicly available, making it a useful tool for 21-cm parameter studies.},
+    issn = {0035-8711},
+    doi = {10.1111/j.1365-2966.2010.17731.x},
+    url = {https://doi.org/10.1111/j.1365-2966.2010.17731.x},
+    eprint = {https://academic.oup.com/mnras/article-pdf/411/2/955/4099991/mnras0411-0955.pdf},
+}
+
+
+# reference for the accretion value of the baseline
+@article{10.1093-mnras-stt1338,
+    author = {Dekel, A. and Zolotov, A. and Tweed, D. and Cacciato, M. and Ceverino, D. and Primack, J. R.},
+    title = {Toy models for galaxy formation versus simulations},
+    journal = {Monthly Notices of the Royal Astronomical Society},
+    volume = {435},
+    number = {2},
+    pages = {999-1019},
+    year = {2013},
+    month = {08},
+    abstract = {We describe simple useful toy models for key processes of galaxy formation in its most active phase, at z \> 1, and test the approximate expressions against the typical behaviour in a suite of high-resolution hydro-cosmological simulations of massive galaxies at z=4−1. We address in particular the evolution of (a) the total mass inflow rate from the cosmic web into galactic haloes based on the EPS approximation, (b) the penetration of baryonic streams into the inner galaxy, (c) the disc size, (d) the implied steady-state gas content and star formation rate (SFR) in the galaxy subject to mass conservation and a universal star formation law, (e) the inflow rate within the disc to a central bulge and black hole as derived using energy conservation and self-regulated Q ∼ 1 violent disc instability (VDI) and (f) the implied steady state in the disc and bulge. The toy models provide useful approximations for the behaviour of the simulated galaxies. We find that (a) the inflow rate is proportional to mass and to (1 + z)5/2, (b) the penetration to the inner halo is ∼50 per cent at z=4−2, (c) the disc radius is ∼5 per cent of the virial radius, (d) the galaxies reach a steady state with the SFR following the accretion rate into the galaxy, (e) there is an intense gas inflow through the disc, comparable to the SFR, following the predictions of VDI and (f) the galaxies approach a steady state with the bulge mass comparable to the disc mass, where the draining of gas by SFR, outflows and disc inflows is replenished by fresh accretion. Given the agreement with simulations, these toy models are useful for understanding the complex phenomena in simple terms and for back-of-the-envelope predictions.},
+    issn = {0035-8711},
+    doi = {10.1093/mnras/stt1338},
+    url = {https://doi.org/10.1093/mnras/stt1338},
+    eprint = {https://academic.oup.com/mnras/article-pdf/435/2/999/3468084/stt1338.pdf},
+}
+
+
+
+# evidence for better halo growth modelling requirements
+@article{McBride_2009,
+   title={Mass accretion rates and histories of dark matter haloes},
+   volume={398},
+   ISSN={1365-2966},
+   url={http://dx.doi.org/10.1111/j.1365-2966.2009.15329.x},
+   DOI={10.1111/j.1365-2966.2009.15329.x},
+   number={4},
+   journal={Monthly Notices of the Royal Astronomical Society},
+   publisher={Oxford University Press (OUP)},
+   author={McBride, James and Fakhouri, Onsi and Ma, Chung-Pei},
+   year={2009},
+   month=oct, pages={1858–1868}
+}
+
+
+
+# Size constraints for reionization simulations
+@article{Iliev_2014,
+   title={Simulating cosmic reionization: how large a volume is large enough?},
+   volume={439},
+   ISSN={0035-8711},
+   url={http://dx.doi.org/10.1093/mnras/stt2497},
+   DOI={10.1093/mnras/stt2497},
+   number={1},
+   journal={Monthly Notices of the Royal Astronomical Society},
+   publisher={Oxford University Press (OUP)},
+   author={Iliev, Ilian T. and Mellema, Garrelt and Ahn, Kyungjin and Shapiro, Paul R. and Mao, Yi and Pen, Ue-Li},
+   year={2014},
+   month=jan, pages={725–743}
+}
+
+
+@ARTICLE{astropy:2022,
+       author = {{Astropy Collaboration} and {Price-Whelan}, Adrian M. and {Lim}, Pey Lian and {Earl}, Nicholas and {Starkman}, Nathaniel and {Bradley}, Larry and {Shupe}, David L. and {Patil}, Aarya A. and {Corrales}, Lia and {Brasseur}, C.~E. and {N{"o}the}, Maximilian and {Donath}, Axel and {Tollerud}, Erik and {Morris}, Brett M. and {Ginsburg}, Adam and {Vaher}, Eero and {Weaver}, Benjamin A. and {Tocknell}, James and {Jamieson}, William and {van Kerkwijk}, Marten H. and {Robitaille}, Thomas P. and {Merry}, Bruce and {Bachetti}, Matteo and {G{"u}nther}, H. Moritz and {Aldcroft}, Thomas L. and {Alvarado-Montes}, Jaime A. and {Archibald}, Anne M. and {B{'o}di}, Attila and {Bapat}, Shreyas and {Barentsen}, Geert and {Baz{'a}n}, Juanjo and {Biswas}, Manish and {Boquien}, M{'e}d{'e}ric and {Burke}, D.~J. and {Cara}, Daria and {Cara}, Mihai and {Conroy}, Kyle E. and {Conseil}, Simon and {Craig}, Matthew W. and {Cross}, Robert M. and {Cruz}, Kelle L. and {D'Eugenio}, Francesco and {Dencheva}, Nadia and {Devillepoix}, Hadrien A.~R. and {Dietrich}, J{"o}rg P. and {Eigenbrot}, Arthur Davis and {Erben}, Thomas and {Ferreira}, Leonardo and {Foreman-Mackey}, Daniel and {Fox}, Ryan and {Freij}, Nabil and {Garg}, Suyog and {Geda}, Robel and {Glattly}, Lauren and {Gondhalekar}, Yash and {Gordon}, Karl D. and {Grant}, David and {Greenfield}, Perry and {Groener}, Austen M. and {Guest}, Steve and {Gurovich}, Sebastian and {Handberg}, Rasmus and {Hart}, Akeem and {Hatfield-Dodds}, Zac and {Homeier}, Derek and {Hosseinzadeh}, Griffin and {Jenness}, Tim and {Jones}, Craig K. and {Joseph}, Prajwel and {Kalmbach}, J. Bryce and {Karamehmetoglu}, Emir and {Ka{l}uszy{'n}ski}, Miko{l}aj and {Kelley}, Michael S.~P. and {Kern}, Nicholas and {Kerzendorf}, Wolfgang E. and {Koch}, Eric W. and {Kulumani}, Shankar and {Lee}, Antony and {Ly}, Chun and {Ma}, Zhiyuan and {MacBride}, Conor and {Maljaars}, Jakob M. and {Muna}, Demitri and {Murphy}, N.~A. and {Norman}, Henrik and {O'Steen}, Richard and {Oman}, Kyle A. and {Pacifici}, Camilla and {Pascual}, Sergio and {Pascual-Granado}, J. and {Patil}, Rohit R. and {Perren}, Gabriel I. and {Pickering}, Timothy E. and {Rastogi}, Tanuj and {Roulston}, Benjamin R. and {Ryan}, Daniel F. and {Rykoff}, Eli S. and {Sabater}, Jose and {Sakurikar}, Parikshit and {Salgado}, Jes{'u}s and {Sanghi}, Aniket and {Saunders}, Nicholas and {Savchenko}, Volodymyr and {Schwardt}, Ludwig and {Seifert-Eckert}, Michael and {Shih}, Albert Y. and {Jain}, Anany Shrey and {Shukla}, Gyanendra and {Sick}, Jonathan and {Simpson}, Chris and {Singanamalla}, Sudheesh and {Singer}, Leo P. and {Singhal}, Jaladh and {Sinha}, Manodeep and {Sip{H{o}}cz}, Brigitta M. and {Spitler}, Lee R. and {Stansby}, David and {Streicher}, Ole and {{{S}}umak}, Jani and {Swinbank}, John D. and {Taranu}, Dan S. and {Tewary}, Nikita and {Tremblay}, Grant R. and {Val-Borro}, Miguel de and {Van Kooten}, Samuel J. and {Vasovi{'c}}, Zlatan and {Verma}, Shresth and {de Miranda Cardoso}, Jos{'e} Vin{'i}cius and {Williams}, Peter K.~G. and {Wilson}, Tom J. and {Winkel}, Benjamin and {Wood-Vasey}, W.~M. and {Xue}, Rui and {Yoachim}, Peter and {Zhang}, Chen and {Zonca}, Andrea and {Astropy Project Contributors}},
+        title = "{The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package}",
+      journal = {The Astrophysical Journal},
+     keywords = {Astronomy software, Open source software, Astronomy data analysis, 1855, 1866, 1858, Astrophysics - Instrumentation and Methods for Astrophysics},
+         year = 2022,
+        month = aug,
+       volume = {935},
+       number = {2},
+          eid = {167},
+        pages = {167},
+          doi = {10.3847/1538-4357/ac7c74},
+archivePrefix = {arXiv},
+       eprint = {2206.14220},
+ primaryClass = {astro-ph.IM},
+       url = {https://ui.adsabs.harvard.edu/abs/2022ApJ...935..167A},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
+# evidence for inside out reionization
+@article{10.1111-j.1365-2966.2006.10502.x,
+    author = {Iliev, I. T. and Mellema, G. and Pen, U.-L. and Merz, H. and Shapiro, P. R. and Alvarez, M. A.},
+    title = {Simulating cosmic reionization at large scales – I. The geometry of reionization},
+    journal = {Monthly Notices of the Royal Astronomical Society},
+    volume = {369},
+    number = {4},
+    pages = {1625-1638},
+    year = {2006},
+    month = {07},
+    abstract = {We present the first large-scale radiative transfer simulations of cosmic reionization, in a simulation volume of (100 h−1 Mpc)3. This is more than a two orders of magnitude improvement over previous simulations. We achieve this by combining the results from extremely large, cosmological, N-body simulations with a new, fast and efficient code for 3D radiative transfer, c2-ray, which we have recently developed. These simulations allow us to do the first numerical studies of the large-scale structure of reionization which at the same time, and crucially, properly take account of the dwarf galaxy ionizing sources which are primarily responsible for reionization. In our realization, reionization starts around z∼ 21, and final overlap occurs by z∼ 11. The resulting electron-scattering optical depth is in good agreement with the first-year Wilkinson Microwave Anisotropy Probe (WMAP) polarization data. We show that reionization clearly proceeded in an inside-out fashion, with the high-density regions being ionized earlier, on average, than the voids. Ionization histories of smaller-size (5–10 comoving Mpc) subregions exabit a large scatter about the mean and do not describe the global reionization history well. This is true even when these subregions are at the mean density of the universe, which shows that small-box simulations of reionization have little predictive power for the evolution of the mean ionized fraction. The minimum reliable volume size for such predictions is ∼30 Mpc. We derive the power spectra of the neutral, ionized and total gas density fields and show that there is a significant boost of the density fluctuations in both the neutral and the ionized components relative to the total at arcmin and larger scales. We find two populations of H ii regions according to their size, numerous, mid-sized (∼10-Mpc) regions and a few, rare, very large regions tens of Mpc in size. Thus, local overlap on fairly large scales of tens of Mpc is reached by z∼ 13, when our volume is only about 50 per cent ionized, and well before the global overlap. We derive the statistical distributions of the ionized fraction and ionized gas density at various scales and for the first time show that both distributions are clearly non-Gaussian. All these quantities are critical for predicting and interpreting the observational signals from reionization from a variety of observations like 21-cm emission, Lyα emitter statistics, Gunn–Peterson optical depth and small-scale cosmic microwave background secondary anisotropies due to patchy reionization.},
+    issn = {0035-8711},
+    doi = {10.1111/j.1365-2966.2006.10502.x},
+    url = {https://doi.org/10.1111/j.1365-2966.2006.10502.x},
+    eprint = {https://academic.oup.com/mnras/article-pdf/369/4/1625/3799304/mnras0369-1625.pdf},
+}
+
+
+

results.typ

@@ -1,59 +1,62 @@
 #import "importer/main.typ": *
 #import "helpers.typ": *
 
-= Results <results>
-This section presents the results of the different simulation runs. We compare the effect of different accretion models on the global signal, map-level differences and statistical properties of the 21-cm brightness temperature field. We focus on three different implementations:
-- The fiducial model where the accretion rate is kept fixed independently of the halo and the redshift. This corresponds to the original implementation of #beorn where $alpha = 0.79$.
+= Impact of individual mass history modeling <results>
+This section presents the outputs of the different simulation runs. We compare the effect of different accretion models on the global signal, map-level differences, and statistical properties of the 21-cm brightness temperature field. We focus on three different implementations:
+- The fiducial model where the accretion rate is kept fixed, independently of the halo and the redshift. This corresponds to the original implementation of #beorn where $alpha = 0.79$.
 - A model where the accretion rate is computed individually for each halo based on its mass growth history and is considered during the painting of each halo.
 - A model where the accretion rate is computed individually for each halo but the considered value during the painting is set to the mean accretion rate of all halos at the respective redshift (effectively reducing the dynamic range of accretion rates).
 
 
-== Effect on the global signal
+== Impact on the global signal
 #let notebook = json("../workdir/11_visualization/simulation_signals.ipynb")
 #figure(
   image_cell(notebook, cell_id: "signal_comparison"),
   caption: [
-    Signal comparison between full runs with the different accretion models: Single value of $alpha$ for all halos according to the mean accretion rate (blue), individual accretion rates for each halo allowing a range from $alpha = 0$ to $alpha = 5$ (green), from $alpha = 0$ to $alpha = 2$ (yellow), and the previously model fixing $alpha = 0.79$ (purple, dashed).
+    Signal comparison between full runs with the different accretion models: Single value of $alpha$ for all halos according to the mean accretion rate (blue), individual accretion rates for each halo allowing a range from $alpha = 0$ to $alpha = 5$ (green), from $alpha = 0$ to $alpha = 3$ (yellow), and the previously model fixing $alpha = 0.79$ (purple, dashed).
     From _left_ to _right_:
     Evolution of the value of the coupling coefficient $x_alpha$.
     Evolution of the mean kinetic temperature $T_k$.
     History of the mean ionization fraction $x_"HII"$.
     Global evolution of the differential brightness temperature $d T_"b"$.
-    The bottom row shows the difference to the reference model - in this case we chose the model following the mean accretion rate. The comparison of the original model with fixed $alpha = 0.79$ is omitted for clarity.
+    The bottom row shows the relative difference to the reference model, chosen to be the model following the mean accretion rate. The comparison of the original model with fixed $alpha = 0.79$ is omitted for clarity.
   ],
 ) <fig:global_signal_combined>
 
-We first investigate the effect of the different accretion models on the global, i.e. averaged, quantities that consititute the 21-cm signal. @fig:global_signal_combined shows the evolution of the coupling coefficient $x_alpha$, the kinetic temperature $T_k$, the ionization fraction $x_"HII"$, and their combined effect on the differential brightness temperature $d T_"b"$. Moving away from the initial model where $alpha = 0.79$ for all halos, we see a clear delay in the evolution of all quantities. This is expected since the accretion rates are overall lower when computed individually for each halo. The more interesting comparison is between the simulation using the moving mean accretion rate and the one using the individual accretion rates. That is the difference which we illustrate in the bottom row of @fig:global_signal_combined. We see that heating is delayed by $Delta z approx 0.5$ whereas the coupling strength is initially lower but increases more rapidly at later times. This could be due to select high mass halos also experiencing high accretion and shifting the balance. This also explains the nearly identical ionization history since
-// TODO - HOW?
+We first investigate the effect of the different accretion models on the global, i.e. the averaged quantities that consititute the 21-cm signal. @fig:global_signal_combined shows the evolution of the coupling coefficient $x_alpha$, the kinetic temperature $T_k$, the ionization fraction $x_"HII"$, and their combined effect on the differential brightness temperature $d T_"b"$. Moving away from the initial model where $alpha = 0.79$ for all halos, we see a clear delay in the evolution of all quantities. This is expected since the accretion rates are overall lower when computed individually for each halo. The refined signals are noticeably less stable and show more fluctuations at early times. This is due to fluctuations in the tree fitting and is inherently linked to the difficulties of the #thesan simulation to resolve the earliest halos properly.
 
-Finally, these effects are summarized by differential brightness temperature: The absorption trough is shifted to later times because the cosmic dawn is delayed.
-This delayed heating results in a lower temperature. Even though the coupling is strong, the spin temperature remains closer to the CMB temperature, leading to a shallower absorption feature.
-The subsequent transition to emission is also delayed but drops to zero more rapidly, which is expected because the end of reionization occurs simultaneously for all models.
+The more interesting comparison is then between the simulation using the moving mean accretion rate and the ones using the individual accretion rates. That is the difference illustrated in the bottom row of @fig:global_signal_combined. Using the mean accretion model as a reference, we compare the two remaining model that consider individiual mass histories over $n=10$ snapshots. The first model allows a range of $0 <= alpha <= 5.0$ and the "reduced" model allows $0 <= alpha <= 3.0$.
 
-This comparison shows that even though the ionization history is largely unaffected by our refined treatment, the global signal is sensitive to the accretion model in ways that cannot be represented by only shifting the global accretion rate. An individual treatment of halos is the key to capture these effects.
+For the individual models we see that heating is delayed by $Delta z approx 0.5$ and the coupling strength is initially lower but increases more rapidly at later times. This could be due to the apparition of more and more high accretion halos which contribute more to the signal while the mean remains low because of the overall increase in halo number. Halos with high accretion rates also have a significant impact the ionization fraction due to the formation of large bubbles, which explains the closely matched ionization histories.
+
+Finally, a summary of these effects is seen in the differential brightness temperature $d T_"b"$: The absorption trough is shifted to later times because the cosmic dawn is delayed. The delayed heating from late star formation results in a lower temperature. Even though the coupling is strong, the spin temperature remains closer to the CMB temperature, leading to a shallower absorption feature.
+The subsequent transition to emission is now less delayed already. In high accretion halos the star formation is increased and heating above the CMB temperature is faster.
+Finally the emission is shorter and drops to zero more rapidly, which is expected because the end of reionization occurs simultaneously for all models.
+
+This comparison shows that even though the evolution of the ionization fraction is largely unaffected by our refined treatment, the global signal is nevertheless sensitive to the accretion model in ways that cannot be represented by only shifting the global accretion rate. An individual treatment of halos is the key to capture these effects.
 
 
 == Map-level investigation
 #let notebook = json("../workdir/11_visualization/simulation_maps.ipynb")
 
-Having established that the individual accretion model produces a distinct global signal, we now compare the map-level differences directly. For a fixed snapshot in time the original model and the model using the mean will create very similar maps since they use the same generalized trend. We therefore directly use the snapshot from the mean model as our reference so that the comparisons are not tainted by the timing differences to the original model.
-
-// TODO change map labels in figures
+Having established that the individual accretion model produces a distinct global signal, we now compare the map-level differences directly. When only considering a single fixed snapshot in time the original model and the model using the mean will create very similar maps since they use the same generalized trend. We therefore directly use a snapshot from the mean model as our reference so that the comparisons are not tainted by the timing differences to the original model.
 
 #figure(
   caption: [
-    Map slices of the core profiles applied onto the simulation grid for the different accretion models plotted at a fixed ionization fraction of $x_"HII" = 0.5$. From _top_ to _bottom_:
+    Map slices of the core profiles applied onto the simulation grid for the different accretion models plotted at a fixed ionization fraction of $x_"HII" = 0.5$. We compare the model that uses a single value of $alpha$ for all halos according to the mean accretion rate (left), individual accretion rates for each halo allowing a range from $alpha = 0$ to $alpha = 5$ (middle), and from $alpha = 0$ to $alpha = 3$ (right).
+
+    From _top_ to _bottom_:
     Map of the $x_alpha$ coupling coefficient and residual map when compared to the reference.
     Map of the kinetic temperature $T_k$ and residual map when compared to the reference.
     Map of the ionization fraction $x_"HII"$ and residual map when compared to the reference.
     In the residual maps blue regions correspond to values lower than the reference model while red regions are higher than the reference model.
   ]
 )[
-  #set image(height: 90%)
+  #set image(height: 87%)
   #image_cell(notebook, cell_id: "grids_and_diffs")
 ] <fig:grids_and_diffs>
 
-@fig:grids_and_diffs shows slices through the simulation box for the different accretion models. We explicitly fix the ionization fraction of $x_"HII" = 0.5$ which removes the effect of different timing of reionization. Thus we can focus on the spatial differences and to compare the morphology of the ionized regions
+@fig:grids_and_diffs shows slices through the simulation box for the different accretion models. We explicitly fix the ionization fraction of $x_"HII" = 0.5$ which removes the effect of different timing of reionization. Thus we can focus on the spatial differences and compare the morphology of the ionized regions
 #footnote[
   Since the models compared here all have a similar ionization history, the redshifts are identical in this case.
 ].
@@ -62,25 +65,18 @@ The maps resemble each other closely and we focus on the residual maps
 // rename?
 that highlight specific deviations produced when changing the accretion model. They show that fixing the mean accretion rate is not sufficient to fully represent the complex reionization behavior.
 
-The coupling coefficient map sees a decrease in all regions which is explained by an overall lower star formation rate compared to the reference case where $alpha = 0.56$. Only a select few halos with higher mass accretion rates produce a positive difference, which suggests that the bulk of the halos behaves similarly but that both positive and negative deviations occur.
+The coupling coefficient map sees an increase in all regions which is explained by the stronger emission of Lyman-$alpha$ photons at these late stages. At the considered redshift, the mean model uses $alpha approx 0.51$. Regions where this reversal has not occured are in theory possible but don't seem to appear with this particular halo population.
 
-This observation is reinforced by the kinetic temperature maps. Many regions are colder than in the fiducial case due to the lower heating by fewer stars. Nevertheless, some regions clearly stand out as being hotter than in the fiducial case. Again, the mean accretion rate model is not able to capture these differences.
-// The background temperature is higher in both alternative accretion models since fixing the ionization fraction means that we show earlier redshifts where the universe has not yet adiabatically cooled as much.
+As before, the kinetic temperature maps reflect the observation made for the signals. Most regions are colder than in the fiducial case due to the lower heating by fewer stars. The halos lag behind but some high accretion halos seem to catch up already, they practically vanish in the residual map.
+These differences are only visible because the mean model fails to capture this diversity of halo histories.
 
-Finally the ionization maps show the clearest differences due to the sharp bubble cutoff.
-// The ionized regions are more compact in the alternative accretion models and the strands of ionized gas connecting the larger bubbles are less pronounced.
-This hints towards a more individual ionization history where large structures of contiguous ionized gas are less common.
-// When we consider the difference to the fiducial model we see that using the mean accretion rate already captures this distinction well.
-There are  multiple bubbles where the detailed mass accretion history generates a clear contrast
-// find something better than "contrast"
-compared to the mean model. Capturing the diversity of halo histories is therefore important to generate maps with the realistic dynamic range.
-// not a fan of "dynamic" here.
+Finally, the ionization maps show the clearest differences due to the sharp bubble cutoff. There are multiple bubbles where the detailed mass accretion models generate a clear contrast, both positive and negative, compared to the mean model. The global picture remains largely unchanged, bubbles have formed and they are in the process of merging into larger contiguous structures. However, the individual bubbles show a nuanced morphology as a direct consequence of the individual halo histories. Capturing this diversity is therefore important to generate maps with the realistic range of existing structures.
 
+So far we have treated the model with a high $alpha$ range and the model with a lower range equally. Focusing on the differences between these two models throughout the previous maps, we see that they are mostly in agreement. This is expected since most halos are expected to have moderate accretion rates. The few very high accretion halos do however lead to small but visible differences that are most easily spotted in the ionization maps.
 
 #figure(
   caption: [
-    Map slices of the brightness temperature $d T_"b"$ for the different accretion models plotted at a fixed ionization fraction of $x_"HII" = 0.5$.
-    // of the 21 cm brightness temperature $d T_"b"$ at
+    Map slices of the brightness temperature $d T_"b"$ for the different accretion models plotted at a fixed ionization fraction of $x_"HII" = 0.5$. The layout is the same as in @fig:grids_and_diffs.
   ]
 )[
   #set image(width: 80%)
@@ -88,25 +84,26 @@ compared to the mean model. Capturing the diversity of halo histories is therefo
 ] <fig:dtb_maps>
 
 
-We give special attention
-// reformulate
-to the derived brightness temperature map. As a reminder, this is not a direct output of the simulation but the spatial distribution can be obtained from the local values of the simulated quantities via @eq:dTb. We present map slices and their comparison to the mean model in @fig:dtb_maps, as previously done for the individual fields. Our observations are compounded here since the contrasts of the fields are combined.
-// A little sentence describing the changes explicitly
+The derived brightness temperature maps are of particular interest. As a reminder, these are not a direct output of the simulation but the spatial distribution can be obtained from the local values of the simulated quantities via @eq:dTb. They correspond to the actual observations that can be made by 21-cm surveys.
+We present map slices and their comparison to the mean model in @fig:dtb_maps, as previously done for the individual fields.
 
+The obtained differences are a combination of the previously discussed effects. The brightness temperature in unionized regions of the IGM remains lower due to the overall lower heating. The bubbles themselves, which appear as dark regions without signal, display a range of morphologies. Some bubbles have grown larger due to the presence of high accretion halos while others are smaller. Beyond the immediate boundary of the bubbles, the temperature seems to be affected as well, as noticeable from the faint ring-like structures around some bubbles. The detailed maps show the richness of structures that can be obtained when considering individual halo histories. The subtle differences induced by this better modeling are expected to eventually be resolved and should definitely be taken into account when interpreting future observations.
 
-== Effect on statistic properties
-#let notebook = json("../workdir/11_visualization/simulation_signals.ipynb")
+We conclude this section by commenting on the differences between the two individual accretion models. The local fluctuations are hard to quantify from the maps alone. An inspection of the resulting power spectra, which exceeds the scope of this report, reveals a small increase in power at small scales when allowing for a larger range of accretion rates. Fluctuations like these are expected to not be visible in the signal measurements but increasing the dynamic range has no detrimental effect on the performance, which is why we recommend the wider range as a default choice.
 
-We also compare summary statistics of the $d T_b$ field. The time evolution of the power spectrum describes
-// what exactly?
+// == Comparison of summary statistics
+// #let notebook = json("../workdir/11_visualization/simulation_signals.ipynb")
 
+// #figure(
+//   caption: [
+//     Power spectra of the brightness temperature for the four models:
+//     Single value of $alpha$ for all halos according to the mean accretion rate (blue), individual accretion rates for each halo allowing a range from $alpha = 0$ to $alpha = 5$ (green), from $alpha = 0$ to $alpha = 3$ (yellow), and the previously model fixing $alpha = 0.79$ (purple, dashed).
+//     _Left_: Power spectra as a function of redshift for fixed scales of $k = 0.12 "Mpc"^(-1)$ and $k = 0.58"Mpc"^(-1)$ (dashed).
+//     ]
+// )[
+//   #set image(width: 70%)
+//   #image_cell(notebook, cell_id: "power_spectra_comparison")
+// ] <fig:power_spectra_comparison>
 
-// BIGG TODO
+// We briefly discuss how the power spectra of the brightness temperature are affected by the different accretion models.
 
-#figure(
-  caption: [
-    ]
-)[
-  #set image(width: 80%)
-  #image_cell(notebook, cell_id: "power_spectra_comparison"),
-] <fig:power_spectra_comparison>

template/template.typ

@@ -69,6 +69,7 @@
   // add space for heading
   show heading: it => v(0.3em) + it + v(0.3em)
   show heading.where(level: 1): it => it + v(0.3em)
+  show heading.where(level: 4): it => text(style: "italic", weight: "regular", it)
 
 
   set math.equation(numbering: "(1)", supplement: [Eq.])
@@ -84,12 +85,12 @@
 
 
 
-  let authors_block(authors, denomination: "Author") = {
-    if authors.len() == 0 {
+  let authors_block(people, denomination: "Author") = {
+    if people.len() == 0 {
       return
     }
     let prefix = denomination
-    if authors.len() > 2 {
+    if people.len() > 1 {
       prefix += "s"
     }
 
@@ -98,8 +99,8 @@
       text(prefix + ":  ", weight: 600),
       stack(
         dir: ttb,
-        spacing: 0.5em,
-        ..authors
+        spacing: 0.8em,
+        ..people
       )
     )
   }
@@ -391,8 +392,11 @@
   //
   // "First" page - abstract and TOC
   //
+  pad(x: 1.5em)[
+    #par(justify: true,  abstract)
+  ]
+
 
-  abstract
   v(2em)
   outline()
 

validation.typ

@@ -3,7 +3,7 @@
 
 
 = Validation <validation>
-We perform several validation tests to ensure the accuracy and reliability of our simulation results. This section details how we check that the new implementation remains robust and consistent with previous versions.
+We perform several validation tests to ensure the accuracy and reliability of our simulation results. This section details how we verify that the new implementation remains robust and consistent with previous versions.
 
 == Code validation
 #let notebook = json("../workdir/11_visualization/validation_simple_run.ipynb")
@@ -14,14 +14,12 @@ We perform several validation tests to ensure the accuracy and reliability of ou
     Global evolution of the differential brightness temperature $d T_"b"$.
     Evolution of the mean kinetic temperature $T_k$.
     History of the mean ionization fraction $x_"HII"$.
-    Dimensionless power spectrum of the $d T_"b"$ field as a function of redshift at $k = 0.12 "Mpc"^(-1)$.
+    21-cm power spectrum as a function of redshift at $k = 0.12 "cMpc"^(-1)$.
   ],
 
 ) <fig:validation_signal_comparison_old_v_new>
 
-We ensure consistency of the updated #beorn code with previous versions by running a series of simulations under identical conditions. We compare key outputs starting from the profiles of individual sources, to the ionization maps, and finally to the global reionization signals. This step-by-step comparison allows us to identify any discrepancy that may arise from the code changes. In @fig:validation_signal_comparison_old_v_new we present a comparison of the central observables of reionization between the previous and current versions of #beorn. The underlying simulation parameters are kept identical or updated to an equivalent setting where necessary. The results show an essentially unchanged evolution of the global brightness temperature, kinetic temperature, and ionization fraction. Note that this run uses an artifical halo catalog and an unphysically high star formation efficiency to ensure a rapid reionization within the simulation range. As such this setup is not representative of a realistic reionization scenario.
-// but ..
-// TODO - comment on the power spectra!
+We ensure the consistency of the updated #beorn code with previous versions by running a series of simulations under identical conditions. We compare key outputs starting from the profiles of individual sources, to the ionization maps, and finally to the global reionization signals. This step-by-step comparison allows us to identify any discrepancy that may arise from the code changes. In @fig:validation_signal_comparison_old_v_new we present a comparison of the central observables of reionization between the previous and current versions of #beorn. The underlying simulation parameters are kept identical or updated to an equivalent setting where necessary. The results show an essentially unchanged evolution of the global brightness temperature, kinetic temperature, and ionization fraction. The small changes in timing are attributed to the slightly shifted mass bins of the updated procedure. Note that this run uses an artifical halo catalog and an unphysically high star formation efficiency to ensure a rapid reionization within the simulation range. As such this setup is not representative of a realistic reionization scenario.
 
 Similarly we maintain backward compatibility with the input format used in previous #beorn runs (i.e. snapshots generated by #pkdgrav or `21cmfast`). This allows us to reproduce the earlier runs and match the results as described by @Schaeffer_2023.
 
@@ -32,18 +30,19 @@ Similarly we maintain backward compatibility with the input format used in previ
   image_cell(notebook, cell_id: "validation_signal_comparison"),
   caption: [
     Validation signal comparison between runs using a grid resolution of $128^3$ (blue) and $256^3$ (green) as well as a run using the #thesandark 1 simulation ($2100^3$ particles) with a grid resolution of $128^3$ (yellow).
-    // TODO - specfify each color
     From _left_ to _right_:
     Global evolution of the differential brightness temperature $d T_"b"$.
     Evolution of the mean kinetic temperature $T_k$.
     History of the mean ionization fraction $x_"HII"$.
-    Dimensionless power spectrum of the $d T_"b"$ field as a function of redshift at $k = 0.12 "Mpc"^(-1)$.
+    Dimensionless power spectrum of the $d T_"b"$ field as a function of redshift at $k = 0.12 "cMpc"^(-1)$.
   ],
 
 ) <fig:validation_signal_comparison>
 
-To ensure that our results are not sensitive to the numerical resolution of the simulation, we perform convergence tests. We compare the following variations of resolution: Firstly, we scrutinize the effect of the grid resolution by running simulations with $128^3$ and $256^3$ cells. Secondly, we investigate the impact of the mass resolution by comparing the results obtained from the #smallcaps[Thesan-Dark] 1 and 2 simulations, which have different particle masses, as mentioned in @procedure.
+To ensure that our results are not sensitive to the numerical resolution of the simulation, we perform convergence tests. We compare the following variations of resolution: Firstly, we scrutinize the impact of the grid resolution by running simulations with $128^3$ and $256^3$ cells. Secondly, we investigate how changing the mass resolution affects the results by comparing the outputs obtained from the #smallcaps[Thesan-Dark] 1 and 2 simulations, which have different particle masses, as mentioned in @procedure.
 
-@fig:validation_signal_comparison shows that there is no significant difference when using a finer grid resolution. This is expected since the ionized regions rapidly become larger than the cell size. The comparison between the different mass resolutions exhibits a deviation in the timing of reionization. This effect is expected and documented by @Kannan_2021: The lowest mass halos which are not resolved by #thesandark 2 form small bubbles quickly (as early as $z=10$) and contribute to the ionization budget at early times. We account for this by artificially increasing the lower mass cutoff during these validation runs but an imbalance between the two simulations remains. Globally the shapes of the signals are very similar and the power spectra match closely.
+@fig:validation_signal_comparison shows that there is no significant difference when using a finer grid resolution. This is expected since the ionized regions rapidly become larger than the cell size. The comparison between the different mass resolutions exhibits a deviation in the timing of reionization. This effect is expected and documented by @Kannan_2021: The lowest mass halos which are not resolved by #thesandark 2 form small bubbles quickly (as early as $z=10$) and contribute to the ionization budget at early times. We account for this by artificially increasing the lower mass cutoff during these validation runs but an imbalance between the two simulations remains. Additionally, the higher particle count per halo seems to enable a more consistent tree finding. The lack of good trees at the earliest snapshots is clearly visible as a dip in the $d T_b$ signal of the #thesandark 2 runs. When the accretion rate cannot be inferred from the tree, a fallback value is used instead. This leads to a systematic shift in the signal that abruptly becomes apparent when the fallback value is no longer needed. After observing this effect we tune the fallback value to be in line with the expected accretion rate as shown in @fig:alpha_evolution_vs_redshift. This largely removes the jump in the subsequent runs.
+
+Globally the shapes of the signals are very similar and, apart from the initially delayed cosmic dawn, the power spectra match closely after the discrepancy during the first few snapshots.
 
 An additional convergence test is performed by varying the binning of the accretion rate $alpha$ parameter. Keeping the overall range fixed between $alpha = 0$ and $alpha = 5.0$, we compare the signal generated when using $n = 5, 25, 50$ bins. A bulk of the halos where the accretion rate cannot be fitted is assigned a fallback value. This then falls into one of the bins, causing a global shift of the signal when changing the binning. However, the overall shape and features of the signal remain unchanged from $n = 25$ bins onwards.