with compiled pdf

.gitignore

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

abstract.typ

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

acknowledgements.typ

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

appendix.typ

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

conclusion.typ

@@ -3,12 +3,18 @@
 
 = Conclusion <conclusion>
 
-Upcoming refinements:
+With the highly anticipated detection of the reionization signal and upcoming observations of the conditions of the IGM during the cosmic dawn, the interpretation of these observations requires accurate predictions from theoretical models and simulations. #beorn
+by @Schaeffer_2023
+is a semi-numerical simulation framework that implements the halo model of reionization
+#cite(<schneider2023cosmologicalforecast21cmpower>, form: "normal") #cite(<Schneider_2021>, form: "normal").
+It uses flux profiles to express the emission of radiation by sources in terms of their host halo to efficiently simulate the reionization on large volumes and to generate predictions for the 21-cm signal. It excels in its computational efficiency and flexibility, allowing for fast and flexible execution.
 
-- 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. We also refactored the simulation procedure and achieved a faster execution time by a factor of two. This further enhances the usability of #beorn for large parameter studies.
 
-- The sensitivity of the results to the growth rate suggest that more refined halo finding and growth tracking algorithms should be investigated
-Where is rockstar?
+After validating the new procedure we have shown that the consistent modeling of halo growth produces simulation outputs which have distinct features compared to simpler models. We compared map outputs directly and also analyzed global quantities and their derived signal. The results are sensitive to the distribution of accretion rates, highlighting the importance of careful modeling of the halo growth.
 
 
+Works going beyond this proof-of-concept implementation should utilize more sophisticated history tracking that ensures the consistency of halo properties across mutlitple timesteps (e.g. the `rockstar` halo finder by @Behroozi_2012). We also highlighted the limitations incurred by the coarse mass resolution of the #thesan simulation, which is why subsequent research should be based on higher resolution simulations in order to benefit from accurate accretion rate matching down to the lowest halo masses.
+
+Furthermore, our analysis of halo growth shows that a simple modeling relying on the halo mass alone is insufficient. Many of the radiative properties in the halo model of reionization have similarly been expressed as a function of halo mass. Given the increasing evidence against simple mass-only models such as the stellar-to-halo-mass relation, the refinement of these models using stochasticity or additional halo properties is a promising avenue for future research. The impact of revised models, including our refinement to the halo growth will be subject to a future publication.
+

halo_mass_history.typ

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

helpers.typ

@@ -1,2 +1,12 @@
 // write the beorn name in a monospace font
 #let beorn = raw("BEoRN")
+
+#let lambdacdm = [$Lambda"CDM"$]
+
+#let pkdgrav = smallcaps[pkdgrav3]
+#let thesan = smallcaps[thesan]
+
+#let thesandark = smallcaps[thesan-dark]
+
+
+#let nbody = [N-body]

implementation.typ

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

introduction.typ

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

main.typ

@@ -13,51 +13,46 @@
   image(bytes(changed)),
   image("assets/uzh-logo.svg")
 )
-#set math.equation(numbering: "(1)", supplement: [Eq.])
+
+#import "importer/main.typ": *
+#let notebook = json("../workdir/11_visualization/very_high_res_snapshot.ipynb")
+
+#let front_image = image_cell(notebook, cell_id: "very_high_res_slice")
+// #let front_image = image("assets/bg-desat.jpg")
+
+
+
 #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",),
+  supervisors: ("Prof. Aurel Schneider", "Prof. Alexandre Refregier"),
   affiliation: "ETH Zürich, Universität Zürich",
   abstract: include("abstract.typ"),
   background-color: color.rgb(32, 64, 123),
   logos: logos,
-  // background-image: image("assets/background.png"),
-  background-image: image("assets/bg-desat.jpg"),
-  date: datetime.today().display("[day]. [month repr:long] [year]"),
+  background-image: front_image,
+  date: datetime.today().display("[day] [month repr:long] [year]"),
   font: "FreeSans",
   pre_content: muchpdf(read("assets/declaration-originality.pdf", encoding: none)),
 )
 
 
-
 //
 // Content
 //
 
-// Stars form early #cite(<10.1093>, form: "normal") but  @10.1093 state that they are bright.
-
-
-// #import "importer/main.typ": *
-// #let notebook = json("../workdir/00_baseline/reference_painting.ipynb")
-
-// #code_cell(notebook, cell_id: "handler")
-// #image_cell(notebook, cell_id: "profile_plot")
-// #image_cell(notebook, cell_id: "profile_plot_2")
-
-
 #include "introduction.typ"
 #include "procedure.typ"
 #include "halo_mass_history.typ"
 #include "implementation.typ"
 #include "validation.typ"
 #include "results.typ"
-// #include "outlook.typ"
 #include "conclusion.typ"
 
 #include "acknowledgements.typ"
 
 #bibliography("references.bib", style: "assets/the-astrophysical-journal.csl")
 
+#pagebreak()
 #include "appendix.typ"

outlook.typ

@@ -5,7 +5,7 @@
 
 
 - implement merger tree growth fitting based on a more sophisiticated growth model (e.g. based on the PS formalism)
--
+- experiment weighted binning of alphas (and masses for that matterr)
 
 - The sensitivity of the results to the growth rate suggest that more refined halo finding and growth tracking algorithms should be investigated
 Where is rockstar?

procedure.typ

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

references.bib

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

results.typ

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

template/template.typ

@@ -54,41 +54,53 @@
   set list(indent: 1em, marker: n => [#text(fill: background-color, "•")])
 
 
+  // 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
   set cite(
     form: "prose"
   )
+  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)
+  show heading.where(level: 4): it => text(style: "italic", weight: "regular", it)
+
+
+  set math.equation(numbering: "(1)", supplement: [Eq.])
+
+  set outline(depth: 2)
 
   //
   // Included content
   //
 
-  // figures
-  // set figure.caption(separator: [ --- ], position: top)
-
   // code blocks
   show raw.where(block: true) : it => h(0.5em) + box(fill: background-color.lighten(80%), outset: 0.5em, width: 100%, it) + h(0.5em)
 
 
 
-  let authors_block(authors, denomination: "Author") = {
-    if authors.len() == 0 {
+  let authors_block(people, denomination: "Author") = {
+    if people.len() == 0 {
       return
     }
     let prefix = denomination
-    if authors.len() > 2 {
+    if people.len() > 1 {
       prefix += "s"
     }
 
     stack(
       dir: ltr,
-      text(prefix + ": ", weight: 600),
+      text(prefix + ":  ", weight: 600),
       stack(
         dir: ttb,
-        spacing: 0.5em,
-        ..authors
+        spacing: 0.8em,
+        ..people
       )
     )
   }
@@ -102,12 +114,16 @@
     // the title page should not have any margin, this will be reset in the next page
     margin: 0pt,
   )
-
+[
+  #set image(
+      width: 105%,
+    )
   // set the image first so that it is the lowest layer
-  place(
+  #place(
     bottom,
     background-image
   )
+]
 
   //  define the base widht of a tile, as a tenth of the page width
   let tile_width = 1.51cm
@@ -302,8 +318,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)
-      align(center, text(subtitle, size: 2em, weight: 500))
+      v(2.5em, weak: true)
+      align(center, text(subtitle, size: 1.8em, weight: 500))
     }
     #pad(
       x: 6em,
@@ -376,8 +392,11 @@
   //
   // "First" page - abstract and TOC
   //
+  pad(x: 1.5em)[
+    #par(justify: true,  abstract)
+  ]
+
 
-  abstract
   v(2em)
   outline()
 

validation.typ

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