more results, first corrective fixes

abstract.typ

@@ -1,4 +1,4 @@
 #import "helpers.typ": *
 = Abstract
-We present an updated 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 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.

appendix.typ

@@ -1,7 +1,23 @@
+#import "helpers.typ": *
+#import "importer/main.typ": *
+
 #set heading(outlined: false, numbering: none)
 
 = Appendix
 
-== A - Generation of the cover image
 
-== B - Halo mass function of #smallcaps[Thesan-Dark 1] and #smallcaps[Thesan-Dark 2]
+== A - Halo mass function of #smallcaps[Thesan-Dark 1] and #smallcaps[Thesan-Dark 2]
+
+#let notebook = json("../workdir/11_visualization/halo_mass_function_thesan_1_2.ipynb")
+
+#figure(
+  image_cell(notebook, cell_id: "halo_mass_functions"),
+  caption: [
+  ]
+) <fig:halo_mass_functions>
+
+// TODO - comment
+
+
+== 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.

conclusion.typ

@@ -3,18 +3,16 @@
 
 = Conclusion <conclusion>
 
-Upcoming refinements:
+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
+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.
 
-- implement merger tree growth fitting based on a more sophisiticated growth model (e.g. based on the PS formalism)
+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.
--
 
-- The sensitivity of the results to the growth rate suggest that more refined halo finding and growth tracking algorithms should be investigated
+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.
-Where is rockstar?
 
-- compare directly to the results of THESAN which provides a star formation rate as well - so a consistent comparison can be made
+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??
-#lorem(200)
+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.
-
-#lorem(200)
-
-#lorem(100)

halo_mass_history.typ

@@ -7,7 +7,7 @@ This section shows the impact of the halo growth on the resulting radiation prof
 // Don't like refined simulation
 
 
-== Modelling mass accretion
+== Modeling mass accretion
 
 As described in @hmreio the fundamental assumption of #beorn is the halo model of reionization by @schneider2023cosmologicalforecast21cmpower.
 // no need to recite?
@@ -21,17 +21,19 @@ 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 modelled 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. For the first iteration halo growth was modeled through an 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 `@???` 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 `@???`
+// 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.
 
 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 both with 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.
+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.
 
 
 // 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
@@ -41,8 +43,8 @@ limitation is the assumption of a constant accretion rate $alpha$ for all halos,
 // // 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.
+// @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.
+// // e.g. papers like "2309...." suggest a revised halo mass growth.
 
 
 
@@ -67,7 +69,7 @@ In order to illustrate the necessity of a more precise treatment of the halo mas
 // 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 pursue the stochasticity of the accretion rate since the usage of n-body 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.
+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.
 
 
 
@@ -75,19 +77,19 @@ We note that the dominating factor when considering a distribution is the contri
 
 === Using THESAN
 
-In order to generate precise map-level predictions of the 21cm signal, #beorn combines the halo model of reionization with large scale N-body simulations which provide realistic snapshots of the dark matter distribution. They give a spatial context to the generated profiles.
+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 post-process the #pkdgrav
+As described in @procedure #beorn was initially used to postprocess the #pkdgrav
 // cite!
-simulation suite and 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.
+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.
 
 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 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.
 
-With a box length of $95.5 "cMpc"$ it provides a sufficient volume to avoid box size effects
+With a box length of $95.5 "cMpc"$ the simulation provides a sufficient volume to avoid box size effects
 // CITATION
 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.
@@ -95,41 +97,38 @@ while still allowing us to iterate quickly and test the refined model without ex
 
 
 // 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.
+// @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.
 
 
-// The simulation has two main limitations: First, the mass resolution of $3.12 * 10^7 "M_⊙"$ means that halos below a mass of $10^9 "M_⊙"$ are not resolved. This is particularly relevant as these low mass halos are expected to contribute significantly to the ionizing photon budget at high redshifts #cite(<Kannan_2021>, form: "normal"). To account for this, we use boosted models of star formation efficiency as described in section <sf_efficiency>. Second, the simulation only provides snapshots down to a redshift of $z=5.5$. As reionization is expected to be completed by this time, this does not impact our results.
+// Thesan halo catalog and the motivation to increase the cutoff.
 
-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.
-
-At the same time THESAN low mass halos seem overabundant which is why we use boosted models of star formation efficiency.
 
 
 
 === Main progenitor branch
 #let notebook = json("../workdir/11_visualization/show_trees.ipynb")
 
-Growth of structure in #lambdacdm is hierarchical: Small structures form first and merge to form larger structures. The growth of halos is reflected in merger trees. The central representation of halo mass evolution is given by merger trees.
+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.
-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.
-// Has this been shown to be true somewhere?
+// 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. 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.
 
 
 === Fitting procedure
-The restriction to the main progenitor means that we reduce the dimensionality of the mass history to a one-dimensional function of redshift compatible with the orginal assumption of an exponential growth model @eq:exponential_growth.
+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.
 
 #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 starting at $z = 10.3$ and looking back over $n=10$ snapshots. Select histories and their corresponding exponential fits are highlighted.
+    _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$.
   ]
 
@@ -137,10 +136,7 @@ The restriction to the main progenitor means that we reduce the dimensionality o
 
 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$.
 
-// TODO - determination of lookback
+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.
-How we determine lookback
-
-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. This means that 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.
 
 
 
@@ -150,18 +146,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: [
-    #lorem(50)
+    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$.
   ]
 ) <fig:alpha_evolution_vs_redshift>
-In order to obtain a sensible range of 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: far from all halos are well represented in the merger tree. Many of the histories are unphysical or incomplete and we describe their treatment 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.
 
-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 this also helps to absorb short term fluctuations introduced most likely by the halo finder. Another advantage of longer lookbacks is the reduced scatter in the resulting distribution which is due to the removal of incomplete, highly fluctuating trees. We note that these behaviors stabilize once we consider aroun 10 snapshots of lookback. Both the mean and standard deviation follow a stable trend and the mean settles at $alpha approx 0.6$.
+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.
 
-// comment on the high tail at high redshifts
+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.
-// explain how the dip at z=14 is probably more meaningful
-// comment on the contrast to $alpha = 0.79$ chosen before
 
-Should explain why different lookbacks produce slightly offset estimates
+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$.
-// Is this due to the fact that we discard cutoff trees (mostly due to _low_ masses)
+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.
 
-We attribute the slight offset of the means to the fact that the 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 lower 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.

helpers.typ

@@ -7,3 +7,6 @@
 #let thesan = smallcaps[thesan]
 
 #let thesandark = smallcaps[thesan-dark]
+
+
+#let nbody = [N-body]

implementation.typ

@@ -2,42 +2,41 @@
 
 = Implementation of changes <implementation>
 
-This section describes the adaptations that were necessary in order to utilize the individual treatment of halo mass accretion histories in #beorn. We distinguish between necessary changes that were required to reflect the underlying model and secondary changes that affect the quality of the simulation outputs indirectly.
+This section describes the adaptations that were necessary in order to utilize the individual treatment of halo mass accretion histories in #beorn. We distinguish between necessary changes that were required to implement the underlying model and secondary changes that affect the quality of the simulation outputs indirectly.
 
 == Profile generation taking into account halo mass history
-For each halo we require a flux profile that matches the halo properties which now include the accretion rate additionally to the mass and the redsift. The profiles are generated in a preprocessing step following the redshifts of the snapshots and the mass and accretion bins defined in the configuration. Since the dynamic range of accretion rates is large the resulting parameter space rapidlly expands. The computation of the profiles therefore utilizes vectorized operations to achieve reasonable runtimes.
+For each halo we require a flux profile that matches the halo properties which now include the accretion rate additionally to the mass and the redshift. The profiles are generated in a preprocessing step following the redshifts of the snapshots and the mass and accretion bins defined in the configuration.
+// 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
-Note that this introduces another "second degree" inconsistency: the flux profile attributes the halo a radiative behavior that is motivated by its history. This is repeated for each snapshot creating possible 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.
+// 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.
 
 
 == Parallel binned painting
-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. Through the addition of the accretion rate as a parameter we introduce a degeneracy that reduces the number of halos that can be treated simultaneously even though their mass is identical. To mitigate this effect we implement a parallelized version of the painting step that distributes the workload onto multiple processes
+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
 #footnote[
-  A rudimentary parallel implementation using `MPI` already exists. It leverages the fact that each snapshot can processed independently and distributes the snapshots onto multiple processes.
+  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.
 ].
-This implementation utilizes a shared memory approach and uses processes on a single node that use 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.
+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.
 
-Some of the painting procedure remains inherently sequential: the final ionization map enforces conservation of the photon count by distributing duplicate ionizations to neighboring cells and 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: 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.
 
 == Merger tree processing
-Fundamental changes include:
+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.
--
-Treatment of invalid trees with negative growth, etc.
-- checked that high alpha halos are rare and balid to be ignored since their history is erratic - probably misidentified progenitors
 
-#lorem(100)
+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 #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$.
 
 
 == Secondary changes
-#beorn was very opinionated in its assumptions and initial data. Since we intend it to create fast and reusable realisations we adapted the code to be more easily adjustable.
+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.
-- better io (proper hdf5 handling, cache)
-- better loading
-- refactoring for modularity
-- refined outputs for testing + validation
-- reduction of runtime
 
-Usage of `Pylians` for speedup @Pylians
+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.
-#lorem(100)
+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.

introduction.typ

@@ -3,40 +3,34 @@
 
 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.
 // Citation about an overview paper on ionization vs structure formation.
-Despite the milestones achieved in observational cosmology
+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.
-// Citation about CMB measurements, JWST, etc.
-, 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 the critical link between the late-time universe and the primordial conditions.
-// This period is crucial as it sets the stage for the subsequent evolution of the universe, including the formation of galaxies and large-scale structures.
-
-The epoch of reionization (EOR) spans the period from the end of the dark ages to the universe becoming fully ionized again. It simultaneously is affected by the fundamental mechanisms
-// reformulate
-and also affects the subsequent evolution of the universe.
-
-$=>$ reionization can serve as a constraint on cosmological models.
-// Paper by aurel on that
 
+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
+@Mao_2008
+@McQuinn_2006
+@schneider2023cosmologicalforecast21cmpower
+).
 Understanding and being able to model the EOR is therefore crucial for a comprehensive picture of cosmology.
 
-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. 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)$.
+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.
+// => argument that there is no "galaxy" in that sense below
 
-molecular cooling as a "workaround"
+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).
-but molecular hydrogen is destroyed by radiation from stars
-=> argument that there is no "galaxy" in that sense below
 
-The end of the dark ages is marked by the formation of the first generation of stars, called population III stars
 // Citation about Pop III stars and their role in the cosmic dawn.
 which...
 During the cosmic dawn ...
 
 
 
-The large amounts of neutral hydrogen in the intergalactic medium (IGM) 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
+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.
-
-
 
 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.
+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
 #cite(<Kannan_2021>, form: "normal")
@@ -45,20 +39,20 @@ State of the art simulations need to implement a range of physical processes, in
 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.
+These methods are computationally expensive, which limits their applicability for large-scale simulations.
 
 // 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 1-d profiles generated by this model to paint the 3-d space around sources which are obtained from a large scale N-body 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.
+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 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 modelled and its impact on the profiles.
+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 modelling assumed by #beorn and the steps required to produce a full 3-d lightcone simulation.
+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.
 // self-consistent treatment of mass accretion.
-@validation details the validation we perform on the refined procedure and in @results we compare the resulting signal to quantify the impact of mass accretion. @conclusion summarizes our findings and discusses potential future improvements.
+@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.
 
 
 
@@ -66,133 +60,6 @@ In @implementation we give an overview of the implementation of the modelling as
 Other points to mention
 - wouthuysen
 - cold reionization
-- IGM - introduce acronym
+- comoving distances - check consistency
 
 how #beorn compares to traditional approaches
-
-
-
-
-
-
-// The hyperfine transition of neutral hydrogen generates photons at
-// the wavelength of 21 cm, opening a new observational window into
-// the early Universe approximately one billion years after the Big
-// Bang. During this era, the radiation from the first stars and galaxies
-// pushes the spin temperature out of equilibrium before heating and
-// eventually ionising the neutral hydrogen of the intergalactic medium
-// (IGM). Next to the source properties, the 21-cm signal depends on
-// the clustering and temperature distribution of the neutral gas, the
-// primordial background radio emission, and the detailed interaction
-// processes between radiation and matter. It is therefore not surprising
-// that the 21-cm radiation from the cosmic dawn contains a wealth of
-// information about the properties of the first stars (Fialkov & Barkana
-// 2014; Mirocha et al. 2018; Ventura et al. 2023; Sartorio et al. 2023),
-// galaxies (Park et al. 2019; Reis et al. 2020; Hutter et al. 2021), and
-// black holes (Pritchard & Furlanetto 2007; Ross et al. 2019). It can
-// furthermore be used to constrain the cosmological model (Liu &
-// Parsons 2016; Schneider et al. 2023; Shmueli et al. 2023) and, in
-// particular, the dark sector, such as the nature of dark matter (Sitwell
-// et al. 2014; Chatterjee et al. 2019; Nebrin et al. 2019; Muñoz et al.
-// 2020; Jones et al. 2021; Giri & Schneider 2022; Hotinli et al. 2022;
-// Flitter & Kovetz 2022; Hibbard et al. 2022), interactions between
-// the dark and visible sector (Barkana et al. 2018; Fialkov et al. 2018;
-// Kovetz et al. 2018; Lopez-Honorez et al. 2019; Mosbech et al. 2023),
-// or potential exotic decay and annihilation processes (D’Amico et al.
-// 2018; Liu & Slatyer 2018; Mitridate & Podo 2018).
-// Reliable detection of the 21-cm signal at these redshifts has yet to
-// be achieved, but ongoing experiments, such as the Giant Metrewave
-// Radio Telescope (GMRT, Paciga et al. 2013), the Precision Array for
-// Probing the Epoch of Reionization (PAPER, Kolopanis et al. 2019),
-// the Murchison Widefield Array (MWA, Trott et al. 2020), the Low-
-// Frequency ARray (LOFAR, Mertens et al. 2020), and the Hydrogen
-// Epoch of Reionization Array (HERA, The HERA Collaboration et al.
-// 2023) have provided upper limits on the 21-cm power spectrum for
-// a broad range of redshifts. These bounds have been used to exclude
-// regions of the parameter space describing extreme properties of the
-// IGM during the epoch of reionisation (Ghara et al. 2020, 2021; Greig
-// et al. 2021a,b; The HERA Collaboration et al. 2022a).
-// The Square Kilometre Array (SKA), a next-generation radio in-
-// terferometer, is currently under construction in South Africa and
-// Western Australia. Its low-frequency component, SKA-low, has the
-// capability to not only measure the 21-cm power spectrum with high
-// signal-to-noise ratio but also provide sky images at redshifts around2
-// T. Schaeffer et al.
-// 𝑧 ≈ 5 − 25 (e.g. Mellema et al. 2015; Wyithe et al. 2015; Ghara
-// et al. 2017; Giri et al. 2018a; Bianco et al. 2021b). The potential
-// of SKA-low for studying the cosmic dawn and reionization era has
-// been extensively investigated in various studies, exploring properties
-// of the ionizing sources and the ionization structure of the universe
-// (e.g. Giri et al. 2018b; Zackrisson et al. 2020; Giri & Mellema 2021;
-// Gazagnes et al. 2021; Bianco et al. 2023). These studies highlight
-// the significant role that SKA-low will play in advancing our under-
-// standing of these critical cosmic epochs.
-// Next to the tremendous experimental effort, accurate and reliable
-// theoretical methods to model the 21-cm signal at the required accu-
-// racy level are currently being developed. Modelling the 21-cm signal
-// is challenging as it involves a broad dynamical range from minihaloes
-// to cosmological scales. It depends on the details of hydrodynamical
-// feedback processes for galaxies, the propagation of radiation through
-// large cosmological scales, and the detailed interaction processes of
-// photons with gas particles of the IGM (e.g., Iliev et al. 2006; Mellema
-// et al. 2006b; Trac & Cen 2007).
-// One option is to predict the 21-cm signal with the help of coupled
-// radiative-transfer hydrodynamic simulations, some well-known ex-
-// amples being the Cosmic Dawn (CoDA) (Ocvirk et al. 2016; Ocvirk
-// et al. 2020; Lewis et al. 2022), the 21SSD (Semelin et al. 2017),
-// and the THESAN simulations (Kannan et al. 2022; Garaldi et al.
-// 2022). Another option is to post-process N-body simulations with
-// ray-tracing algorithms, such as the Conservative, Causal Ray-tracing
-// code (C2 RAY; Mellema et al. 2006a) or the Cosmological Radiative
-// transfer Scheme for Hydrodynamics (CRASH; Maselli et al. 2003).
-// Full radiative-transfer numerical methods are fundamental to un-
-// derstanding the 21-cm signal and estimating the accuracy of more
-// approximate methods. However, they are very computationally ex-
-// pensive and can hardly be used to scan the vast cosmological and
-// astrophysical parameter space. To perform Bayesian inference anal-
-// ysis on a mock 21-cm data set, semi-numerical algorithms are often
-// used, better suited to generate thousands of realizations of the sig-
-// nal itself. They rely on the excursion set formalism (Furlanetto et al.
-// 2004), such as 21cmFAST (Mesinger et al. 2011) or SIMFAST21 (San-
-// tos et al. 2010).
-// In this paper, we present the new framework BEoRN which stands
-// for Bubbles during the Epoch of Reionisation Numerical simulator.
-// The code is based on a one-dimensional radiative transfer method
-// in which interactions between matter and radiation are treated in a
-// spherically symmetric way around sources. This approach is signifi-
-// cantly faster than full 3-d radiative transfer codes and arguably more
-// precise than semi-numerical algorithms which are not based on indi-
-// vidual sources. In this aspect, BEoRN is similar to other existing codes
-// such as BEARS (Thomas et al. 2009) or GRIZZLY (Ghara et al. 2018).
-// However, in contrast to other 1d radiative transfer codes, BEoRN self-
-// consistently accounts for the evolution of individual sources during
-// the emission of photons. This includes both the redshifting of pho-
-// tons due to the expansion of space and the increase of luminosity
-// caused by the growth of individual sources over time. Both effects
-// have a non-negligible influence on the radiation profile surrounding
-// sources.
-// The BEoRN framework allows for a flexible parametrisation to
-// model any source of radiation, such as e.g. Pop-III stars, galaxies,
-// or quasars. It produces a 3-dimensional (3D) light-cone realisation
-// of the 21-cm signal from the cosmic dawn to the end of reionisation
-// including redshift space distortion effects. The underlying gas density
-// field as well as the position of sources is directly obtained from
-// outputs of an 𝑁-body simulation. We have designed BEoRN to be
-// user-friendly and modular so that it can be applied in combination
-// with different gravity solvers or source models, for example.
-// MNRAS 000, 1–18 (2023)
-// The paper is structured as follows: Section 2 describes the BEoRN
-// code, while section 3 validates it by comparing its predictions with
-// the publicly available 21cmFAST code. In section 4, three benchmark
-// models are presented, calibrated to the latest observations, and the
-// evolution of the 21-cm signal during the cosmic dawn and epoch
-// of reionization is studied. The work concludes with a summary and
-// conclusion in section 5.
-// Note that throughout the paper, physical distance units are specified
-// with the prefix "𝑝", while co-moving distance units are specified
-// with the prefix "𝑐". The cosmological parameters used in this work
-// are consistent with Planck 2018 results (Planck Collaboration et al.
-// 2020), namely matter abundance Ωm = 0.31, baryon abundance
-// Ωb = 0.045, and dimensionless Hubble constant ℎ = 0.68. The
-// standard deviation of matter perturbations at 8ℎ −1 cMpc scale is
-// 𝜎8 = 0.81.

main.typ

@@ -23,8 +23,8 @@
 
 
 #show: tasteful-thesis.with(
-  // title: "BEoRN version 2",
   title: "Simulating the EOR with self-consistent growth of galaxies",
+  subtitle: "Master's Thesis",
   authors: ("Rémy Moll",),
   supervisors: ("Prof. Aurel Schneider",),
   affiliation: "ETH Zürich, Universität Zürich",
@@ -32,13 +32,12 @@
   background-color: color.rgb(32, 64, 123),
   logos: logos,
   background-image: front_image,
-  date: datetime.today().display("[day]. [month repr:long] [year]"),
+  date: datetime.today().display("[day] [month repr:long] [year]"),
   font: "FreeSans",
   pre_content: muchpdf(read("assets/declaration-originality.pdf", encoding: none)),
 )
 
 
-
 //
 // Content
 //

procedure.typ

@@ -4,46 +4,46 @@
 
 = Overview of the #beorn framework <procedure>
 
-This section describes the model describing the sources of radiation that drive reionization. We explain how #beorn implements this model to generate 3D maps of the IGM during the epoch of reionization. The code of #beorn as well as usage instructions are publicly available under #link("https://github.com/cosmic-reionization/BEoRN", "https://github.com/cosmic-reionization/BEoRN")#footnote[
+This section presents the model describing the sources of radiation that drive reionization. We explain how #beorn implements this model to generate 3-d maps of the IGM during the epoch of reionization. The code of #beorn as well as usage instructions are publicly available under #link("https://github.com/cosmic-reionization/BEoRN", "https://github.com/cosmic-reionization/BEoRN")#footnote[
   For an explicit overview of the changes referenced here, please refer to #link("https://github.com/moll-re/BEoRN")
 ].
 
 
 == The halo mass model of reionization <hmreio>
 
-The central action
+The central feature
 // don't like that word
-performed by #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 halos. 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 allows for 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 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 modelling
+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
-of the halo mass evolution is discussed in section @halo_mass_history, the model itself simply considers an arbitrary but known halo mass accretion history $M_"h" (z)$.
+of the halo mass evolution is discussed in @halo_mass_history. The model itself simply considers an arbitrary but known halo mass accretion history $M_"h" (z)$.
 
 The emission of radiation revolves around the star formation rate $dot(M)_star$ which is simply assumed to be proportional to the halo mass accretion rate via
 $
   dot(M)_star = f_star (M_"h") dot dot(M_"h")
 $ <eq:star_formation_rate>
 where the star formation efficiency $f_star$ introduces a mass dependence that enables the suppression of star formation in low mass halos and the implementation of a cooling limit.
-Motivated by abundance matching @schneider2023cosmologicalforecast21cmpower use the double power law
+Motivated by abundance matching, @schneider2023cosmologicalforecast21cmpower use the double power law
 $
   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.
+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.
-
 
+// TODO - make clear that this follows @Schneider
 
 === Expression of the profiles
-Deriving from the star formation rate the halo model predicts the production and distribution of photons in 3 distinct energy bands:
+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, ionizing UV photons, and X-ray photons. Each of these bands has a different effect on the IGM and is treated separately.
 // Not really sure that's true
 
-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
+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")
+#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 absorption of 21-cm photons before reionization. This is reflected in the absorption expected in the global signal before reionization.
 $
   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)
 $
@@ -77,34 +77,35 @@ $
 
 === 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 of the foreground with the CMB background.
+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.
 // 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.
 
 
 == 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 section @implementation.
+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.
 
-=== Halo catalog - n body simulations
+=== 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
+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.
 
 === Computation of radiation profiles
-In accordance with the astrophysical parameters set by the user, radiation profiles are computed to be applied onto halos in a subsequent step. We deliberately separate the computation of the profiles 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
+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
 #footnote[
   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 3D grid. This is done by iterating over the halos in the catalog and using their corresponding profile. For a given profile a 3D kernel is generated and applied onto the grid via convolution. 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.
+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.
 
 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 3D nbody simulation means that the results are sensitive to the underlying cosmology and the detailed profile modelling means that the results are sensitive to the underlying astrophysical model.
+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.

references.bib

@@ -199,3 +199,69 @@ archivePrefix = {arXiv},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+
+# ROCKSTAR
+@article{Behroozi_2012,
+   title={THE ROCKSTAR PHASE-SPACE TEMPORAL HALO FINDER AND THE VELOCITY OFFSETS OF CLUSTER CORES},
+   volume={762},
+   ISSN={1538-4357},
+   url={http://dx.doi.org/10.1088/0004-637X/762/2/109},
+   DOI={10.1088/0004-637x/762/2/109},
+   number={2},
+   journal={The Astrophysical Journal},
+   publisher={American Astronomical Society},
+   author={Behroozi, Peter S. and Wechsler, Risa H. and Wu, Hao-Yi},
+   year={2012},
+   month=dec, pages={109}
+}
+
+
+
+
+# Constraints on cosmology from reionization
+
+@article{Mao_2008,
+   title={How accurately can 21 cm tomography constrain cosmology?},
+   volume={78},
+   ISSN={1550-2368},
+   url={http://dx.doi.org/10.1103/PhysRevD.78.023529},
+   DOI={10.1103/physrevd.78.023529},
+   number={2},
+   journal={Physical Review D},
+   publisher={American Physical Society (APS)},
+   author={Mao, Yi and Tegmark, Max and McQuinn, Matthew and Zaldarriaga, Matias and Zahn, Oliver},
+   year={2008},
+   month=jul
+}
+
+
+@article{McQuinn_2006,
+   title={Cosmological Parameter Estimation Using 21 cm Radiation from the Epoch of Reionization},
+   volume={653},
+   ISSN={1538-4357},
+   url={http://dx.doi.org/10.1086/505167},
+   DOI={10.1086/505167},
+   number={2},
+   journal={The Astrophysical Journal},
+   publisher={American Astronomical Society},
+   author={McQuinn, Matthew and Zahn, Oliver and Zaldarriaga, Matias and Hernquist, Lars and Furlanetto, Steven R.},
+   year={2006},
+   month=dec, pages={815–834}
+}
+
+
+
+# importance of RSDs
+@article{Ross_2021,
+   title={Redshift-space distortions in simulations of the 21-cm signal from the cosmic dawn},
+   volume={506},
+   ISSN={1365-2966},
+   url={http://dx.doi.org/10.1093/mnras/stab1822},
+   DOI={10.1093/mnras/stab1822},
+   number={3},
+   journal={Monthly Notices of the Royal Astronomical Society},
+   publisher={Oxford University Press (OUP)},
+   author={Ross, Hannah E and Giri, Sambit K and Mellema, Garrelt and Dixon, Keri L and Ghara, Raghunath and Iliev, Ilian T},
+   year={2021},
+   month=jul, pages={3717–3733}
+}

results.typ

@@ -2,91 +2,111 @@
 #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:
+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$.
 - 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
 #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$ (purble, 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 = 2$ (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$.
-    Mean ionization fraction history $x_"HII"$.
+    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.
   ],
 ) <fig:global_signal_combined>
 
-We first compare 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 - since the overall star formation rate is reduced,
+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?
+
+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.
+
+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.
 
 
 == 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
+
 #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 of the $x_alpha$ coupling coefficient. Difference to the fiducial model.
+    Map of the $x_alpha$ coupling coefficient and residual map when compared to the reference.
-    Map of the kinetic temperature $T_k$. Difference to the fiducial model.
+    Map of the kinetic temperature $T_k$ and residual map when compared to the reference.
-    Map of the ionization fraction $x_"HII"$. Difference to the fiducial model.
+    Map of the ionization fraction $x_"HII"$ and residual map when compared to the reference.
-    In the difference plots blue regions correspond to values lower than the fiducial model while red regions are higher than the fiducial model.
+    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%)
-  #image_cell(notebook, cell_id: "grids_and_diffs"),
+  #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
+#footnote[
+  Since the models compared here all have a similar ionization history, the redshifts are identical in this case.
+].
+We omit the original model with $alpha = 0.79$ and directly compare the two alternative accretion models.
+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.
 
-@fig:grids_and_diffs shows slices through the simulation box for the different accretion models. Instead of fixing the redshift, we show slices at a fixed ionization fraction of $x_"HII" = 0.5$. This helps to remove the effect of different timing of reionization and focus on the spatial differences and to compare the morphology of the ionized regions. Since the slices we pick are identical the maps resemble each other closely so we focus on the difference plots that highlight specific deviations produced when changing the accretion model.
+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 maps sees a decrease in most regions which is explained by an overall lower star formation rate compared to the fiducial case where $alpha = 0.79$. Only a select few halos with high mass accretion rates produce a positive difference. Comparing directly between the individually modelled halos and the case where only the mean accretion rate is used, we see that the bulk of the halos are similar but that both positive and negative deviations occur. This shows that simply using the mean accretion rate is not sufficient to capture map level differences.
+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.
 
-This observation is repeated in the kinetic temperature maps. Many regions are colder than in the fiducial case due to the lower heating by fewer stars. Still 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.
+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.
-Finally the ionization maps show the clearest differences due to the modelling using a step function. 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 however multiple bubbles where the detailed mass accretion history generates a clear contrast compared to the mean model. Capturing the diversity of halo histories is therefore important to generate maps with the realistic dynamic range.
+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.
 
 
 #figure(
-  image_cell(notebook, cell_id: "results_lightcones"),
   caption: [
-    // Lightcone images showing brightness temperature slices as they evolve with time.
-    // Blue regions absorb the background radiation while orange regions correspond to IGM that has been heated above the background temperature. Black regions are ionized and do not affect the brightness temperature.
     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
   ]
-) <fig:results_lightcones>
+)[
+  #set image(width: 80%)
+  #image_cell(notebook, cell_id: "dtb_maps")
+] <fig:dtb_maps>
 
 
-// TODO - show a sense of scale
-
 We give special attention
 // reformulate
-to the derived brightness temperature map. 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:???`. @fig:results_lightcones shows slices and their comparison to the fiducial model, as previously done for the individual fields. Our observations are compounded here since the contrasts of the fields are combined.
+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
 
 
 == Effect on statistic properties
+#let notebook = json("../workdir/11_visualization/simulation_signals.ipynb")
 
 We also compare summary statistics of the $d T_b$ field. The time evolution of the power spectrum describes
 // what exactly?
 
-// #figure(
-//   image_cell(notebook, cell_id: "power_spectra_comparison"),
-//   caption: [
-//     ]
 
-// ) <fig:power_spectra_comparison>
+// BIGG TODO
 
-
+#figure(
-#lorem(50)
+  caption: [
-
+    ]
-Importance of RSD for the 21 cm signal
+)[
-https://arxiv.org/abs/2011.03558
+  #set image(width: 80%)
+  #image_cell(notebook, cell_id: "power_spectra_comparison"),
+] <fig:power_spectra_comparison>

template/template.typ

@@ -57,6 +57,7 @@
   // figure placement and caption style
   set figure(placement: top)
   show figure.caption: set text(size: 0.85em)
+  show figure.caption: set align(left)
 
 
   // citation style
@@ -66,11 +67,14 @@
   show cite: it => text(fill: background-color, it)
 
   // add space for heading
-  show heading.where(level:1): it => it + v(0.5em)
+  show heading: it => v(0.3em) + it + v(0.3em)
+  show heading.where(level:1): it => it + v(0.3em)
 
 
   set math.equation(numbering: "(1)", supplement: [Eq.])
 
+  set outline(depth: 2)
+
   //
   // Included content
   //
@@ -313,8 +317,8 @@
     #set text(font: font, fill: font-color)
     #align(center, text(title, size: 2.5em, weight: 600))
     #if subtitle != none {
-      v(1.5em, weak: true)
+      v(2.5em, weak: true)
-      align(center, text(subtitle, size: 2em, weight: 500))
+      align(center, text(subtitle, size: 1.8em, weight: 500))
     }
     #pad(
       x: 6em,

validation.typ

@@ -12,17 +12,19 @@ We perform several validation tests to ensure the accuracy and reliability of ou
   caption: [
     Validation signal comparison between the new implementation of #beorn (blue) and the old version (green, dashed). From _left_ to _right_:
     Global evolution of the differential brightness temperature $d T_"b"$.
-    Evolution of the mean kinetic temperature of the gas $T_k$.
+    Evolution of the mean kinetic temperature $T_k$.
-    Mean ionization fraction history $x_"HII"$.
+    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)$.
   ],
 
 ) <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 individuall 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 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.
+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 ...
+// but ..
+// TODO - comment on the power spectra!
+
+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.
 
-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.
 
 == Convergence tests
 #let notebook = json("../workdir/11_visualization/validation_convergence.ipynb")
@@ -33,15 +35,15 @@ Similarly we maintain backward compatibility with the input format used in previ
     // TODO - specfify each color
     From _left_ to _right_:
     Global evolution of the differential brightness temperature $d T_"b"$.
-    Evolution of the mean kinetic temperature of the gas $T_k$.
+    Evolution of the mean kinetic temperature $T_k$.
-    Mean ionization fraction history $x_"HII"$.
+    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)$.
   ],
 
 ) <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 compare effects of the grid resolution by running simulations with $128^3$ and $256^3$. 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 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.
 
-@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 sees 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. Globally the shapes of the signals are very similar and the power spectra match closely.
 
-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. Since a bulk of the halos where the accretion rate cannot be fitted is assigned a fallback value which then falls into one of the bin, there is a global shift to the signal when changing the binning. However, the overall shape and features of the signal remain unchanged from $n = 25$ bins onwards.
+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.