Growing the first galaxies’ merger trees

ABSTRACT

Modelling the growth histories of specific galaxies often involves generating the entire population of objects that arise in a given cosmology and selecting systems with appropriate properties. This approach is highly inefficient when targeting rare systems such as the extremely luminous high-redshift galaxy candidates detected by JWST. Here, we present a novel framework for generating merger trees with branches that are guaranteed to achieve a desired halo mass at a chosen redshift. This method augments extended Press Schechter theory solutions with constrained random processes known as Brownian bridges and is implemented in the open-source semi-analytic model galacticus. We generate ensembles of constrained merger trees to predict the growth histories of seven high-redshift JWST galaxy candidates, finding that these systems most likely merge ≈2 Gyr after the observation epoch and occupy haloes of mass ≳10¹⁴M_⊙ today. These calculations are thousands of times more efficient than existing methods, are analytically controlled, and provide physical insights into the evolution of haloes with rapid early growth. Our constrained merger tree implementation is publicly available at https://github.com/galacticusorg/galacticus.

galaxies: formation, galaxies: haloes, galaxies: high-redshift, cosmology: theory, dark matter

1 INTRODUCTION

The extended Press–Schechter (ePS) formalism (Press & Schechter 1974; Bond et al. 1991; Bower 1991; Lacey & Cole 1993) provides an analytic recipe for predicting the formation histories of dark matter haloes, both in standard ΛCDM cosmologies (e.g. Sheth, Mo & Tormen 2001; Giocoli et al. 2007) and beyond (e.g. Benson et al. 2013). The ePS formalism is based on the set of excursions δ(S) that a point in the linear overdensity field executes as it is smoothed on successively smaller scales, commonly paramterized by the mass variance S (see Fig. 1 for examples). Mergers in a halo’s growth history correspond to successive maxima of δ(S); the masses and redshifts of these mergers are determined by the values of the mass variance at which the linear overdensity reaches the threshold necessary for halo formation.

Figure 1.

Left: Constrained excursions (coloured lines) and a standard unconstrained excursion (grey line), that represent possible growth histories of a final halo of mass M₀ = 10¹⁴M_⊙ at redshift |$z$|₀ = 0. The constrained Brownian bridge excursions are forced to pass through an overdensity and mass variance corresponding to a halo of mass M₁ = 10¹²M_⊙ at |$z$|₁ = 8 (blue), M₁ = 10¹²M_⊙ at |$z$|₁ = 6 (green), and M₁ = 10¹¹ M_⊙ at |$z$|₁ = 8 (purple). After each constraint is achieved, growth histories are modeled with unconstrained excursions until |$z$|₀. M₀, |$z$|₀, M₁, and |$z$|₁ are free input parameters in our constrained solutions for any M₁ < M₀ and |$z$|₁ > |$z$|₀. Although excursions appear non-monotonic on the redshift axis, only their successive maxima correspond to steps in halo growth (shown in the right-hand panel). Right: Corresponding halo growth histories. The purple and green points from the left-hand panel are shown as lower limits on halo mass at their respective constrained masses and redshifts because these excursions first cross above the corresponding critical overdensities at larger mass scales. However, the constrained merger tree algorithm we present yields growth histories that precisely satisfy desired constraints, like the blue line in the right-hand panel.

Models that leverage the ePS formalism have been calibrated to numerical simulations and are remarkably useful for rapidly generating ensembles of halo growth histories and merger trees (for a review, see Zentner 2007 and references therein). However, existing ePS-based approaches are highly inefficient when modelling specific astrophysical systems. In particular, most ePS algorithms model excursions as uncorrelated random walks and thus cannot guarantee that constraints on a halo’s growth history beyond its present-day mass are satisfied. We will show that the ability to generate merger trees that satisfy specific constraints is desirable for modelling rare astrophysical systems such as the first detectable galaxies. In particular, we derive constrained excursion set solutions that allow us to directly generate the growth histories of rare high-redshift galaxies, rather than generating exponentially many realizations and selecting those with the desired properties. Our constrained merger tree realizations of the first galaxies are up to ∼10⁴ times more efficient than standard ePS calculations for sampling mean growth histories, and up to ∼10⁸ times faster for sampling outliers among this population.

Mathematically, excursions are usually modelled as uncorrelated random walks with normally distributed transition probabilities under the usual assumption of a Gaussian linear overdensity field.¹ However, as emphasized by Zentner (2007), the uncorrelated nature of ePS excursions is not a prediction of the theory; rather, it is a simplifying assumption that may not hold in practice. For example, any window function other than the sharp k-space filter correlates Fourier modes of the density field (Bond et al. 1991), and haloes’ large-scale environments also introduce correlations in their merger histories. Here, we go beyond the usual ePS assumption of uncorrelated random walks by implementing constrained excursions of the linear overdensity field in a merger tree algorithm. To achieve this, we use constrained random processes known as Brownian bridges. Unlike uncorrelated random walks, Brownian bridges are guaranteed to pass through desired points along their trajectories.

Physically, our Brownian bridge solutions allow us to fix a point in a halo’s growth history at a chosen mass scale below its final mass and construct excursions on larger mass scales that are consistent with this constraint. This implies that a branch of a final halo’s merger tree must exceed a desired mass by a specific redshift. In detail, we derive halo merger rates for excursions that begin at (δ₀, S₀) and are guaranteed to pass through a specific point (δ₁, S₁) for any S₁ > S₀ and δ₁ > δ₀. We implement these merger rates in a new merger tree algorithm, which guarantees that a halo with final mass M₀ at redshift |$z$|₀ has a progenitor of exactly mass M₁ at redshift |$z$|₁ for any M₁ < M₀ and |$z$|₁ > |$z$|₀. All four of these quantities (i.e. M₀, |$z$|₀, M₁, and |$z$|₁) are free input parameters in our approach. Thus, halo growth histories generated using our constrained merger tree algorithm may be arbitrarily rare relative to the unconstrained distribution for a final halo of mass M₀ at redshift |$z$|₀. Crucially, the distribution of these constrained trajectories is precisely equal to the subset of unconstrained excursions that satisfy the desired constraints; our predictions are therefore statistically unbiased, even for extremely rare systems.

The left-hand panel of Fig. 1 shows examples of constrained and unconstrained excursions that represent a final halo of mass M₀ = 10¹⁴M_⊙ at |$z$|₀ = 0; the right-hand panel shows the corresponding halo growth histories. More extreme constrained solutions – e.g. the blue model that exceeds M₁ = 10¹² M_⊙ at |$z$|₁ = 8, making it an ≈4σ outlier relative to the underlying halo population at that redshift – reach large overdensities more quickly along their trajectories. Because halo growth can be read off as the successive maxima of an excursion, followed from left to right on the left-hand panel of Fig. 1, such extreme constraints yield rapid early growth before |$z$|₁ followed by relatively quiescent merger histories. On the other hand, standard unconstrained ePS solutions, like the grey excursion in Fig. 1, often yield significant major mergers at late times for this final halo mass, as reflected in the corresponding growth history.

We leverage this ability to generate rare unbiased halo growth histories to grow the first large ensembles of merger trees representing the earliest detectable galaxies; this task would otherwise be computationally intractable. Due to observational selection effects, observed high-redshift galaxies are generally the most luminous objects that exist at early times, which implies that they likely occupy the most rapidly growing haloes at these epochs (e.g. Mason, Trenti & Treu 2023). For example, surprisingly luminous galaxies at redshifts |$z$| ≳ 10 have been identified in Hubble Space Telescope (HST; Oesch et al. 2016) and Hyper-Suprime-Cam Subaru Strategic Program (HSC-SSP; Harikane et al. 2022) data. Recently, similarly remarkable galaxy candidates have been detected in JWST data (e.g. Curtis-Lake et al. 2022; Labbe et al. 2022; Robertson et al. 2022; Naidu et al. 2022a) data, though the nature of these systems is still highly uncertain due to instrumental calibration (e.g. Boyer et al. 2022), dust, initial-mass function, and nebular-emission modelling (e.g. Zavala et al. 2023; Naidu et al. 2022b), as well as limited spectroscopic follow-up.

Estimates based on abundance matching and extreme value statistics suggest that the extremely luminous high-redshift galaxy candidates detected by JWST occupy haloes with masses ≳10¹¹ at |$z$| ≳ 8, implying that they are at least ≈3σ outliers relative to a standard ΛCDM halo population (e.g. Lovell et al. 2022); as a result, the very existence of these early luminous objects may pose challenges to ΛCDM itself (Boylan-Kolchin 2022). Our results allow us to statistically predict the merger rates and descendant mass distributions of these early galaxies, assuming they are hosted by rare haloes in a standard ΛCDM cosmology. Note that, we do not attempt to forward model the baryonic components of rare early haloes in this paper; however, we intend to explore such modelling in future work.

Beyond the high-redshift applications considered in this paper, we anticipate that constrained excursion set solutions will generally be useful for modelling other kinds of astrophysical systems. For example, near-field cosmological observations provide detailed constraints on the merger histories of systems like the Milky Way (e.g. see Helmi 2020 for a review). Such systems are also ‘rare’, in the sense that halo growth histories generated using standard ePS theory are very unlikely to match the specific set of relevant constraints. Thus, even if a system’s growth history is not extreme relative to the mean of the underlying distribution, our algorithm yields merger trees that are guaranteed to match its specific properties and sample the remaining scatter in possible growth histories. We set the stage for such applications in this work by demonstrating that certain choices of our Brownian bridge constraint parameters yield merger events that resemble the infall of the Large Magellanic Cloud into the Milky Way. We also intend to explore other classes of constrained excursion set solutions, along with relevant astrophysical and cosmological applications, in future studies.

This paper is organized as follows. In Section 2, we review essential ingredients of excursion set theory and derive the first-crossing distribution for Brownian bridge excursions; in Section 3, we describe our method for constructing constrained merger trees based on these first-crossing distributions; in Section 4, we present halo growth histories, formation time distributions, and merger rates for constrained merger trees; in Section 5, we use our technique to predict the merger rates and descendants of high-redshift JWST galaxy candidates; in Section 6, we model the Milky Way–Large Magellanic Cloud merger; we address potential caveats in Section 7, and we conclude in Section 8. We also present several Appendices: in Appendix  A, we present the full integral equations that determine our unconstrained and constrained first-crossing rate distributions; in Appendix  B, we show that our constrained solutions converge to the unconstrained distribution of halo growth histories in appropriate limits, and we demonstrate the numerical convergence of our constrained results; in Appendix  C, we derive the distribution of descendant halo masses that enters our high-redshift galaxy calculations, given the |$z$| = 0 halo mass function.

Throughout, we adopt the best-fit Planck cosmology with H₀ = 67.66 km s⁻¹ Mpc⁻¹, Ω_m = 0.3111, Ω_b = 0.049, Ω_Λ = 0.6889, n_s = 0.9665, and σ₈ = 0.8102 (Planck VI Collaboration 2020). Note that, because we construct dark matter-only halo growth histories in this paper, we set Ω_b = 0 when computing halo masses. Throughout, log  denotes the base-10 logarithm.

2 EXCURSION SET THEORY

In this section, we describe our excursion set theory procedure for modelling constrained dark matter halo growth histories. In Section 2.1, we briefly review standard ePS theory ingredients, including an expression for the first-crossing distribution of excursions given the collapse barrier. We use this expression to derive the first-crossing distribution for standard unconstrained excursions in Section 2.2 and our constrained Brownian bridge excursions in Section 2.3. We present examples of the resulting first-crossing rate distributions in Section 2.4.

2.1 Standard ingredients

The first ingredient is the linear power spectrum of matter density fluctuations P(k), which we use to define the mass variance |$\sigma (M)\equiv \sqrt{S(M)}$| (i.e. the fractional root-variance of the linear overdensity field extrapolated to |$z$| = 0), according to

$$ \begin{eqnarray*} S(M)\equiv \sigma ^2(M) = \frac{1}{2\pi ^2}\int _{0}^{\infty }4\pi k^2 P(k) W^2(k|M) \mathrm{d}k, \end{eqnarray*} $$

(1)

where k is the comoving wavenumber and W(k|M) is a window function that smooths the linear overdensity field δ on a mass scale M. We assume a sharp k-space window function

$$ \begin{eqnarray*} W_k(k|M) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 &k\le k_S(M)\\ 0 &k\gt k_S(M), \end{array}\right. \end{eqnarray*} $$

(2)

where the subscript S indicates that the halo mass M is evaluated relative to the root-variance defined by equation (1). In standard ePS theory and in all of our calculations, steps along the trajectory δ(S) correspond to draws from a Gaussian random field in Fourier space and are thus uncorrelated for a sharp k-space window function. However, to compute σ(M), we follow Benson et al. (2013) by using a real-space top-hat window function,

$$ \begin{eqnarray*} W_R(k|M) = \frac{3[\sin (kR)-kR\cos (kR)]}{(kR)^3}, \end{eqnarray*} $$

(3)

where the filtering scale R is related to halo mass M via

$$ \begin{eqnarray*} R = \left(\frac{3M}{4\pi \bar{\rho }}\right)^{1/3}, \end{eqnarray*} $$

(4)

and |$\bar{\rho }$| is the mean matter density at |$z$| = 0 (Lacey & Cole 1993). We use a real-space top-hat window function because it allows for this unambiguous mass definition and corresponds to a physically intuitive scenario in which localized overdensities collapse into haloes.

The next ingredient is the critical threshold overdensity δ_c(⁠|$z$|⁠), which represents a condition for the gravitational collapse of a density perturbation in linear theory. This quantity acts as an absorbing barrier B(S), because trajectories that first up-cross the barrier at variance S collapse into a halo of mass M[S].² We follow Bond et al. (1991) by assuming that

$$ \begin{eqnarray*} B(S,t) = \delta _c(M[S],t)=\frac{\delta _{c,0}}{D(t)}, \end{eqnarray*} $$

(5)

where δ_{c, 0} ≈ 1.686 is the linear-theory collapse threshold and D(t) is the linear growth factor (Carroll, Press & Turner 1992). In particular, following Benson et al. (2013), we use the barrier in equation (5) rather than the ‘remapped’ barrier introduced in Sheth et al. (2001) that is calibrated to match halo mass functions from cosmological simulations because we only aim to compute merger rates rather than halo mass functions.

To predict the progenitors of a halo of mass M₀ at redshift |$z$|₀, we need to compute the conditional distribution of halo masses that would lead to this final mass (i.e. the conditional mass function; Lacey & Cole 1993). In the excursion set approach, this corresponds to solving for the first up-crossing of δ(S) given excursions that begin from the point (δ₀, S₀) corresponding to the final halo in question. As described in Benson et al. (2013), computing merger rates as a function of the relative variance,

$$ \begin{eqnarray*} \tilde{S}\equiv S-S_0, \end{eqnarray*} $$

(6)

amounts to solving the first-crossing problem with a modified barrier

$$ \begin{eqnarray*} \tilde{B}(\tilde{S},t^{\prime },t_0)=B(\tilde{S}+S_0,t^{\prime })-B(S_0,t_0), \end{eqnarray*} $$

(7)

where t₀ is defined implicitly (given δ₀) via equation (5) and t′ < t is an earlier time in this halo’s history. We then estimate the first-crossing rate using the finite-difference method from Benson et al. (2013). In particular, we set t′ = (1 − ϵ)t₀, where ϵ ≪ 1, which implies that the first-crossing rate of the effective barrier |$\tilde{B}(\tilde{S},t^{\prime },t_0)$| is then df/dt = f(S)/ϵt₀, where f(S)dS is the fraction of excursions that first up-cross the barrier between S and S + dS.

Thus, the first-crossing rate distribution df(S)/dt gives the rate at which excursions first achieve the critical overdensity, and physically corresponds to the branching rate of dark matter haloes into lower-mass progenitors. The problem of constructing halo merger histories therefore boils down to finding the first-crossing distribution f(S). Zhang & Hui (2006) formulated a generic integral equation for f(S),

$$ \begin{eqnarray*} 1 = \int _0^S f(\hat{{S}})\mathrm{d}\hat{{S}} + \int _{-\infty }^{\tilde{B}(S)} P(\delta ,S)\mathrm{d}\delta , \end{eqnarray*} $$

(8)

where P(δ, S)dδ is the probability that a trajectory lies between δ and δ + dδ at variance S, having never crossed the barrier at any smaller S. Equation (8) follows from mass conservation, because all trajectories must either (i) cross B(S) at |$0\,\lt\, \hat{\boldsymbol{ S}}\,\lt\, S$|⁠, or (ii) lie below the barrier at S, having never crossed it at smaller S (Benson et al. 2013). P(δ, S) can be further decomposed according to

$$ \begin{eqnarray*} P(\delta ,S) = P_0(\delta ,S)-\int _{0}^S f(\hat{{ S}})P_0[\delta ,S|\tilde{B}(\hat{{ S}}),\hat{{ S}}]\mathrm{d}\hat{{ S}}, \end{eqnarray*} $$

(9)

where the first term is the distribution of trajectories in the absence of a barrier, and the second term is the fraction of trajectories that first up-crossed the barrier at |$\hat{\boldsymbol{ S}}\,\lt\, S$| and then reach (δ, S). Combining equations (8) and (9) yields

$$ \begin{eqnarray*} 1 &= & \int _0^S f(\hat{{ S}})\mathrm{d}\hat{{ S}} \nonumber \\ &&+\, \int _{-\infty }^{\tilde{B}(S)} \left[P_0(\delta ,S)-\int _{0}^S f(\hat{{ S}})P_0[\delta ,S|\tilde{B}(\hat{{ S}}),\hat{{ S}}]\mathrm{d}\hat{{ S}}\right]\mathrm{d}\delta . \end{eqnarray*} $$

(10)

We follow the strategy for numerically solving equation (10) by discretizing the integrals over grids of |$\hat{\boldsymbol{ S}}$|⁠, as described in Appendix A of Benson et al. (2013). Thus, in this approach, we must specify P₀(δ, S) and |$P_0[\delta ,S|\tilde{B}(\hat{\boldsymbol{ S}}),\hat{\boldsymbol{ S}}]$| to define an excursion set solution.

2.2 Unconstrained excursions

In the standard excursion set problem, the distribution of trajectories P₀(δ, S) in the absence of a barrier (i.e. the first term in equation 9) is a normal distribution with zero mean and variance S,

$$ \begin{eqnarray*} P_0(\delta ,S) = \mathcal {N}(\delta ,S) = \frac{1}{\sqrt{2\pi S}}\exp \left(-\frac{\delta ^2}{2S}\right). \end{eqnarray*} $$

(11)

Next, we consider trajectories that start at |$(\tilde{B}(\hat{\boldsymbol{ S}}),\hat{\boldsymbol{ S}})$| and end at (δ, S) to derive the distribution inside the integrand of the second term in equation (9). For uncorrelated random walks, there is no net drift in the trajectories, so the effective difference in δ between these two points is simply

$$ \begin{eqnarray*} \Delta \delta = \delta -\tilde{B}(\hat{{ S}}). \end{eqnarray*} $$

(12)

The covariance between these two points is |$\mathrm{Cov}(\hat{\boldsymbol{ S}},S) = \min (\hat{\boldsymbol{ S}},S) = \hat{\boldsymbol{ S}}$| (since |$\hat{\boldsymbol{ S}}\,\lt\, S$| by construction). Thus, the residual variance is

$$ \begin{eqnarray*} \Delta S = S-\mathrm{Cov}(\hat{{ S}},S) = S-\hat{{ S}}. \end{eqnarray*} $$

(13)

The desired conditional distribution function is therefore

$$ \begin{eqnarray*} P_0[\delta ,S|\tilde{B}(\hat{{ S}}),\hat{{ S}}] &=& \mathcal {N}(\Delta \delta ,\Delta S) \nonumber\\ &=& \frac{1}{\sqrt{2\pi (S-\hat{{ S}})}}\exp \left(-\frac{(\delta -\tilde{B}(\hat{{ S}}))^2}{2(S-\hat{{ S}})}\right). \end{eqnarray*} $$

(14)

Plugging equations (11) and (14) into equation (10) and evaluating the integral in the second term over dδ analytically yields the final form of our integral equation for f(S). This equation is presented explicitly in Appendix A1. We solve for f(S) following the numerical scheme presented in Benson et al. (2013), using the midpoint integration method from Du, Behrens & Niemeyer (2017). We then use the resulting first-crossing distributions to generate unconstrained halo growth histories as described in Section 3. Two free quantities specify the distribution of unconstrained halo growth histories:

• M₀, the final halo mass, which enters equation (10) through the definition |$\tilde{S}=S-S_0$| and through the effective barrier |$\tilde{B}(S)$| implicitly, via δ_c(M[S], t);

• |$z$|₀, the redshift at which the final halo exists, which only enters equation (10) via the dependence of the effective barrier on δ_c(M[S], t). The final redshift |$z$|₀ need not equal zero.

2.3 Constrained Brownian bridge excursions

We construct constrained merger trees with branches that begin at (δ₀, S₀) and are forced to pass through the point (δ₁, S₁) using the Brownian bridge first-crossing rate distribution for S₀ ≤ S ≤ S₁, and then switch to an unconstrained solution for S > S₁. We justify this procedure in Section 3. Note that, Brownian bridge excursions may up-cross δ₁ at S < S₁. Thus, such excursions formally place a lower limit on the mass of the progenitor on a single branch of the merger tree at that redshift; we refer to this as the ‘constrained’ branch. Our implementation of Brownian bridge excursions in galacticus exploits symmetries of the excursion set problem to guarantee that this constrained branch reaches a mass of precisely M₁ at redshift |$z$|₁. We discuss these points in detail below, and focus on the solution for the constrained first-crossing rate here.

For the purposes of computing the distributions that enter our solution for the first-crossing distribution in equation (10), there are two main differences between the Brownian bridge and unconstrained cases. First, unlike uncorrelated random walks, Brownian bridge excursions have a net drift characterized by

$$ \begin{eqnarray*} \mu (S) = \langle \delta (S) \rangle = \delta _0 + \frac{S-S_0}{S_1-S_0}(\delta _1-\delta _0). \end{eqnarray*} $$

(15)

Second, the covariance for two points |$(\tilde{B}(\hat{S}),\hat{S})$| and (δ, S) along a Brownian bridge is

$$ \begin{eqnarray*} \mathrm{Cov}(\tilde{B}(\hat{{ S}}),\delta) = \frac{(S_1-S)(\hat{{ S}}-S_0)}{S_1-S_0}, \end{eqnarray*} $$

(16)

for |$\hat{{ S}}\lt S$|⁠. Note that, as S → S₁, μ(S) → δ₁ and |$\mathrm{Cov}(\tilde{B}(\hat{{ S}}),\delta)\rightarrow 0$|⁠, because we can recursively apply the Brownian bridge drift and covariance as we update δ and S at each step along an excursion. In particular, any trajectory that starts at (S₀, δ₀) and passes through |$(\hat{{ S}},\hat{\delta })$| on its way to (S₁, δ₁), where |$S_0\lt \hat{{ S}}\lt S_1$|⁠, is itself a Brownian bridge that starts at |$(\hat{{ S}},\hat{\boldsymbol{ \delta} })$| and ends at (S₁, δ₁). Thus, steps of mean size δ₁ are forced to occur as S₁ is approached, ensuring that the trajectories up-cross δ₁ before S₁ is reached, as desired.

Next, to compute the effective difference in δ between the points |$(\tilde{B}(\hat{{ S}}),\hat{{ S}})$| and (δ, S) for Brownian bridge excursions, we subtract the appropriate drift term from each point,

$$ \begin{eqnarray*} \Delta \delta ^{\prime } &=& [\delta -\mu (S)] - [\tilde{B}(\hat{{ S}})-\mu (\hat{{ S}})] \nonumber \\ &=& \delta -\tilde{B}(\hat{{ S}})-\frac{S-\hat{{ S}}}{S_1-S_0}(\delta _1-\delta _0). \end{eqnarray*} $$

(17)

Similarly, the residual variance is given by

$$ \begin{eqnarray*} \Delta S^{\prime } &=& \mathrm{Var}(S)-\mathrm{Cov}(\tilde{B}(\hat{{ S}}),\delta) \nonumber \\ &=& \frac{(S_1-S)(S-S_0)}{S_1-S_0}-\frac{(S_1-S)(\hat{{ S}}-S_0)}{S_1-S_0} \nonumber \\ &=& \frac{(S_1-S)(S-\hat{{ S}})}{S_1-S_0}. \end{eqnarray*} $$

(18)

Note that Δδ′ → Δδ and ΔS′ → ΔS in the limit S₁ → ∞. Physically, this limit represents constrained trajectories that must exceed a very small halo mass at |$z$|₁, in which case the constrained first-crossing rate distribution equals that for the unconstrained case, as expected.

With equations (16) and (18) in hand, we can now re-write the distribution functions that appear in equation (9) for constrained Brownian bridge excursions. Specifically, the distribution of trajectories in the absence of a barrier becomes

$$ \begin{eqnarray*} P^{\prime }_0(\delta ,S) &=& \mathcal {N}(\Delta \delta ^{\prime }|_{\hat{S}\rightarrow S_0},\Delta S^{\prime }|_{\hat{S}\rightarrow S_0}) \nonumber \\ & =& \mathcal {N}\left(\delta {-}\tilde{B}(S_0)-\frac{S-S_0}{S_1-S_0}(\delta _1{-}\delta _0),\frac{(S_1-S)(S-S_0)}{S_1-S_0}\right).\nonumber\\ \end{eqnarray*} $$

(19)

The conditional distribution function appearing in equation (9) follows analogously, without taking the |$\tilde{S}\rightarrow S_0$| limit in the expressions for |$\Delta \hat{\boldsymbol{ \delta} }$| and |$\Delta \hat{{ S}}$|

$$ \begin{eqnarray*} &&P^{\prime }_0[\delta ,S|\tilde{B}(\hat{{ S}}),\hat{{ S}}]\nonumber\\ &&\quad = \mathcal {N}(\Delta \delta ^{\prime },\Delta S^{\prime }) \nonumber \\ &&\qquad = \mathcal {N}\left(\delta -\tilde{B}(\hat{S})-\frac{S-\hat{{ S}}}{S_1-S_0}(\delta _1-\delta _0),\frac{(S_1-S)(S-\hat{{ S}})}{S_1-S_0}\right). \nonumber\\ \end{eqnarray*} $$

(20)

Plugging equations (19) and (20) into equation (10) and evaluating the integral in the second term over dδ analytically yields the final form of our integral equation for f(S) in the constrained case. This equation is presented explicitly in Appendix A2, and we solve for f(S) in the Brownian bridge case following the same numerical scheme as for the unconstrained problem. Four free quantities specify the distribution of constrained halo growth histories:

• M₀, the final halo mass, which (as in the unconstrained case) enters equation (10) through the definition |$\tilde{S}=S-S_0$| and through the effective barrier |$\tilde{B}(S)$| implicitly, via δ_c(M[S], t), and (only in the constrained case) enters equations (19) and (20) via S₀;

• |$z$|₀, the redshift at which the final halo exists, which (as in the unconstrained case) enters equation (10) via the dependence of the effective barrier on δ_c(M[S], t), and (only in the constrained case) enters equations (19) and (20) via δ₀;

• M₁, the mass that the constrained trajectory exceeds at redshift |$z$|₁, which enters equations (19) and (20) via S₁;

• |$z$|₁, the redshift at which the constrained trajectory exceeds mass M₁, which enters equations (19) and (20) via δ₁.

2.4 Examples of first-crossing rate distributions

Fig.2 shows examples of the first-crossing rate distribution, obtained using our numerical solution method described in Appendix  A, for a Brownian bridge excursion with M₁ = 10¹²M_⊙ and |$z$|₁ = 8 that reaches M₀ = 10¹⁴M_⊙ at |$z$|₀ = 0 (this corresponds to the blue constraint shown in Fig. 1). We plot the first-crossing rate distribution as a function of ΔS = S[M] − S[M_p], where M_p is the mass of the parent halo, i.e. the halo branching into smaller objects at a given redshift as we work backwards in time.

For early times and large parent halo masses, the constrained first-crossing rate is very similar to the unconstrained solution, except that it cuts off at large ΔS. This is expected, because branching events to much smaller masses at early times would violate the Brownian bridge constraint. On the other hand, as M_p approaches M₁ and |$z$|_p approaches |$z$|₁, the constrained first-crossing rates begin to peak at a ΔS corresponding to the branching event necessary to satisfy the constraint (formally, in the limit M_p → M₁ and |$z$|_p → |$z$|₁, the constrained first-crossing rate approaches a δ function at the appropriate ΔS). These features lend confidence to our solutions for the constrained first-crossing rate.

3 CONSTRUCTING MERGER TREES

We now use the unconstrained and constrained first-crossing rate distributions derived in the previous section to construct merger trees. We describe the standard algorithm for constructing unconstrained merger trees in Section 3.1, our modifications to this algorithm for constructing constrained merger trees in Section 3.2, and our implementation of the modified algorithm in galacticus in Section 3.3. We present examples of constrained and unconstrained merger trees in Section 3.4.

3.1 Standard algorithm

To construct merger trees, we generally follow the procedure described in Benson et al. (2013), which is based on the Parkinson, Cole & Helly (2008) algorithm. We briefly summarize this procedure below in order to provide context for our constrained merger tree algorithm.

3.1.1 Branching rates

In the standard Parkinson et al. (2008) algorithm, the probability per unit time of a binary merger occurring for a halo of mass M is

$$ \begin{eqnarray*} \frac{\mathrm{d}p}{\mathrm{d}\omega } = \int _{M_{\mathrm{res}}}^{M/2} \frac{M}{M^{\prime }}\frac{\mathrm{d}f[M^{\prime }]}{\mathrm{d}t}\frac{\mathrm{d}S[M^{\prime }]}{\mathrm{d}M^{\prime }}\left|\frac{\mathrm{d}t}{\mathrm{d}\omega }\right|G[\omega ,\sigma (M),\sigma (M^{\prime })]\mathrm{d}M^{\prime }, \nonumber\\ \end{eqnarray*} $$

(21)

where ω ≡ δ_{c, 0}/D(t), M_res is the minimum mass resolution of the merger tree, df[M′]/dt is the unconstrained first-crossing rate distribution obtained following the procedure in Appendix A1, and |$G[\omega ,\sigma (M),\sigma (M^{\prime })]\equiv G[\omega ,\sigma _M,\sigma _{M^{\prime }}]$| is an empirical modifier function with parameters calibrated to reproduce the progenitor mass function from cosmological simulations. Throughout, we use a generalized form of the Parkinson et al. (2008) modifier function,

$$ \begin{eqnarray*} G[\omega ,\sigma _1,\sigma _2]=G_0(\sigma _2/\sigma _1)^{\gamma _1}(\omega /\sigma _1)^{\gamma _2}\left(1-{\sigma _2^2 \over \sigma _1^2}\right)^{\gamma _3}, \end{eqnarray*} $$

(22)

with G₀ = 1.14, γ₁ = −0.33, γ₂ = 0.06, and γ₃ = 0.65 based on a preliminary recalibration to recent simulations (A. Benson et al., in preparation). Because we focus on high-mass haloes in this study, the details of the modifier function do not affect our results.

Accretion of haloes below the resolution limit contributes mass to all haloes in the merger tree at a rate

$$ \begin{eqnarray*} \frac{\mathrm{d}R}{\mathrm{d}\omega } = \int _{M_{\mathrm{min}}}^{M_{\mathrm{res}}} \frac{\mathrm{d}f[M^{\prime }]}{\mathrm{d}t}\frac{\mathrm{d}S[M^{\prime }]}{\mathrm{d}M^{\prime }}\left|\frac{\mathrm{d}t}{\mathrm{d}\omega }\right|G[\omega ,\sigma (M),\sigma (M^{\prime })]\mathrm{d}M^{\prime }, \nonumber\\ \end{eqnarray*} $$

(23)

where M_min is a minimum mass element, chosen to be sufficiently small such that the resulting progenitor mass distributions converge at the per cent level (also see Benson et al. 2013).

3.1.2 Merger tree construction

To construct unconstrained merger trees, we select a final halo mass M₀ and redshift |$z$|₀ and we probabilistically generate a series of branching events in which the halo splits into progenitors of lower masses. In particular, when a branching event is determined to occur by comparing a randomly drawn number to the unconstrained branching probability (equation 21), the mass of one of the progenitors is drawn from d²f/dωdM′ over the interval [M_res, M/2], and the remaining progenitor’s mass is determined by mass conservation.

We define the main branch as the branch of the merger tree obtained by choosing the more massive progenitor of each merger event (working backwards in time) until the mass resolution limit is reached. In practice, we first construct the main branch and then revisit all secondary progenitors that were generated during construction of the main branch, repeating this procedure to build out the remainder of the unconstrained merger tree.

3.2 Modifications to the merger tree construction algorithm

We modify the Parkinson et al. (2008) algorithm to build constrained merger trees as follows.

3.2.1 Branching rates

As described by Cole et al. (2000), the integral in equation (21) spans [M_res, M/2] rather than the full mass range [M_res, M] of potential progenitor haloes. We note that this is an empirical choice, based on which ePS merger trees have been calibrated to cosmological simulations, rather than a requirement of the theory. None the less, it is usually a reasonable approximation because halo merger rates derived from the unconstrained branching rate is almost symmetric about M/2 (e.g. see Benson, Kamionkowski & Hassani 2005; Jiang & van den Bosch 2014). However, our constrained excursions yield branching rates that are generally asymmetric about M/2 and can thus break down using the usual [M_res, M/2] interval. For example, as a Brownian bridge excursion approaches the constrained point, the rate approaches a δ-function that guarantees a branching event resulting in the desired progenitor mass M₁. Symmetrizing the branching rate by sampling from [M_res, M/2] is therefore not a valid approximation given our constrained first-crossing rate distributions.

To address this, we compute branching rates using the full mass range of [M_res, M − M_res] throughout this work,

$$ \begin{eqnarray*} \frac{\mathrm{d}p}{\mathrm{d}\omega } = \frac{1}{2}\int _{M_{\mathrm{res}}}^{M-M_{\mathrm{res}}} \frac{M}{M^{\prime }}\frac{\mathrm{d}f[M^{\prime }]}{\mathrm{d}t}\frac{\mathrm{d}S[M^{\prime }]}{\mathrm{d}M^{\prime }}\left|\frac{\mathrm{d}t}{\mathrm{d}\omega }\right| \tilde{G}[\omega ,\sigma _M,\sigma _{M^{\prime }}]\mathrm{d}M^{\prime }, \nonumber\\ \end{eqnarray*} $$

(24)

where the factor of 1/2 is needed to normalize the branching rate because we draw progenitor masses from the full interval.³ In equation (24), |$\tilde{G}[\omega ,\sigma _M,\sigma _{M^{\prime }}]$| is our Parkinson et al. (2008)-like modifier function symmetrized about M/2 using the M < M/2 half of equation (22), that is

$$ \begin{eqnarray*} \tilde{G}[\omega ,\sigma _M,\sigma _{M^{\prime }}] = \left\lbrace \begin{array}{ll}G[\omega ,\sigma _M,\sigma _{M^{\prime }}] & \hbox{if } M^{\prime } \le M/2, \\ G[\omega ,\sigma _M,\sigma _{M-M^{\prime }}] & \hbox{if } M^{\prime } \gt M/2. \end{array} \right. \end{eqnarray*} $$

(25)

We make this choice because the Parkinson et al. (2008) modifier function (and our version of it in equation 22) is calibrated over the interval [0, M/2] under the assumption of binary mergers.⁴ We accordingly modify the sub-resolution accretion rate expression to

$$ \begin{eqnarray*} \frac{\mathrm{d}R}{\mathrm{d}\omega } &=& \frac{1}{2} \Bigg (\int _{M_{\mathrm{min}}}^{M_{\mathrm{res}}} \frac{\mathrm{d}f[M^{\prime }]}{\mathrm{d}t}\frac{\mathrm{d}S[M^{\prime }]}{\mathrm{d}M^{\prime }}\left|\frac{\mathrm{d}t}{\mathrm{d}\omega }\right|\tilde{G}[\omega ,\sigma _M,\sigma _{M^{\prime }}]\mathrm{d}M^{\prime }\nonumber \\ && +\, \int _{M-M_{\mathrm{res}}}^{M-M_{\mathrm{min}}} \frac{\mathrm{d}f[M^{\prime }]}{\mathrm{d}t}\frac{\mathrm{d}S[M^{\prime }]}{\mathrm{d}M^{\prime }}\left|\frac{\mathrm{d}t}{\mathrm{d}\omega }\right|\tilde{G}[\omega ,\sigma _M,\sigma _{M^{\prime }}]\mathrm{d}M^{\prime } \Bigg).\nonumber\\ \end{eqnarray*} $$

(26)

Thus, when a branching event is determined to occur by comparing a randomly drawn number to the branching probability (equation 24), the mass of one of the progenitors is drawn from d²f/dωdM′ over the interval [M_res, M − M_res], and the remaining progenitor’s mass is determined by mass conservation; in the following subsection, we describe how the constrained first-crossing rate distribution is used in these calculations. For consistency, note that, we also calculate branching rates and sub-resolution accretion using equations (24) and (26) (rather than equations 21 and 23) and draw progenitor masses from the full mass interval when constructing unconstrained merger trees (and the unconstrained portions of constrained trees). Given our modifications to the progenitor mass interval, the unconstrained merger trees we construct are not guaranteed to yield halo growth histories that match previous ePS predictions calibrated to cosmological simulations. As discussed below, we defer such recalibration to future work; here, we aim to demonstrate that our constrained merger trees are well converged and satisfy the desired constraints.

3.2.2 Merger tree construction

To construct merger trees that satisfy a Brownian bridge constraint, we proceed as follows. We first select the final halo mass, M₀, and redshift, |$z$|₀, as the starting point for our tree construction. We also select the Brownian bridge constrained mass, M₁, and redshift, |$z$|₁. We then build a specific branch of the tree, referred to as the ‘constrained’ branch, by drawing branching events from our constrained first-crossing rate distribution (obtained by plugging equations 19 and 20 into equation 10). We use the branching and sub-resolution accretion rates in equations (24) and (26), respectively, to build the constrained branch back in time to redshift |$z$|₁.

We then force the constrained branch to split such that one of its progenitors has a mass of precisely M₁ at redshift |$z$|₁. This is allowed due to the properties of our constrained Brownian bridge excursions. In particular, these excursions guarantee that δ₁ is first up-crossed at S ≤ S₁ (for example, see the purple excursion in Fig. 1 for a case that clearly first up-crosses δ₁ at S < S₁). If δ₁ is first up-crossed at S = S₁, then a branching event to a progenitor of mass M₁ at redshift |$z$|₁ is directly predicted to occur. Meanwhile, if the constrained excursion reaches δ₁ at S < S₁ and then steps towards some δ > δ₁, this step can be replaced by an equally likely, reflected step towards smaller δ of equal magnitude, following the classic approach to the ‘cloud-in-cloud’ problem in Bond et al. (1991), based on Chandrasekhar (1943).⁵ Such reflections can be applied for any step that would exceed δ₁ for S < S₁, until S₁ is reached, guaranteeing a branching event to a progenitor of mass precisely M₁ at redshift |$z$|₁. Thus, we effectively create constrained excursions that first up-cross δ₁ at S₁.⁶ Fig. 3 illustrates an example of this reflection principle for possible excursions corresponding to a M₁ = 10¹²M_⊙, |$z$|₁ = 8 constraint.

Figure 2.

First-crossing rate distributions for constrained Brownian bridge excursions with M₁ = 10¹² M_⊙, |$z$|₁ = 8, M₀ = 10¹⁴ M_⊙, and |$z$|₀ = 0, versus standard, unconstrained ePS trajectories with the same final halo mass and redshift (grey). Solutions are shown for parent haloes with M_p = 5 × 10¹³ M_⊙ and |$z$|_p = 1 (solid), M_p = 5 × 10¹²M_⊙ and |$z$|_p = 7 (dashed), and M_p = 5 × 10¹³ M_⊙ and |$z$|_p = 7 (dot–dashed); the dot–dashed and dashed unconstrained distributions are identical. Vertical lines show progenitor masses of 10¹²M_⊙, corresponding to the constrained mass, for M_p = 5 × 10¹³ M_⊙ (dashed) and M_p = 5 × 10¹² M_⊙ (dot–dashed).

Figure 3.

Possible Brownian bridge excursions for a desired constraint of M₁ = 10¹² M_⊙ at |$z$|₁ = 8 (blue point), for a final halo of mass M₀ = 10¹⁴M_⊙ at redshift |$z$|₀ = 0. The dashed pink portion of the excursion first exceeds δ₁ at the pink point, corresponding to a growth history that overshoots the desired constraint. We reflect such excursions about δ₁ to obtain equally probable trajectories, resulting in the solid blue portion of the excursion for the example shown here. This symmetry allows us to insert a merger to a halo of mass exactly M₁ at redshift |$z$|₁ during our constrained merger tree construction. Note that, our constrained merger tree algorithm always uses the reflection principle (corresponding to the blue solid line for the excursion shown here). The bifurcation points are shown for illustrative purposes and do not correspond to physical interactions.

Two conditions must be met for these reflected steps to be equally likely compared to steps that first up-cross δ₁ for S < S₁. First, the barrier must be constant as a function of ΔS, which holds for all scenarios we consider. Second, the Brownian bridge drift term (equation 17) must be zero, which holds because a constrained excursion that has reached (S, δ₁) at S < S₁ is also a Brownian bridge that ends at (S₁, δ₁), with a new starting point of (S, δ₁). Thus, the constrained branch is allowed to split to a progenitor of mass precisely M₁ at redshift |$z$|₁ (working backwards in time) while maintaining the correct progenitor mass distribution for S < S₁. We emphasize that constrained merger trees cannot self-consistently be constructed simply by forcing such branching events to occur using the unconstrained first-crossing rate distribution. In particular, our solution for the constrained first-crossing rate distribution guarantees that the resulting progenitor mass distribution is equivalent to the distribution of unconstrained trajectories down-sampled to satisfy the relevant constraint.

In practice, when a branching event to a progenitor with M > M₁ is drawn along the constrained branch, we manually insert a branching event to a progenitor with a mass of precisely M₁ at a redshift of |$z$|₁. The remainder of the constrained branch (i.e. the |$z$| > |$z$|₁ portion) is constructed using the unconstrained first-crossing rate distribution. We then return to all secondary haloes and build out the merger tree using the unconstrained first-crossing rate distribution. Thus, except for the constrained branch at redshifts |$z$| < |$z$|₁, all branches of the merger tree are built using the unconstrained first-crossing rate distribution. As in the unconstrained case, all branches (including the constrained branch) are built backwards in time until they reach the mass resolution limit.

Our procedure to stitch together constrained and unconstrained first-crossing rate distributions along the constrained branch, at |$z$|₁, warrants discussion. Within the scope of the standard uncorrelated random walk excursion set model, this choice implies that trajectories for S > S₁ are fairly sampled from all possible trajectories that lead to (δ₁, S₁). Note that, for S > S₁, the statistics of excursions obtained by stitching unconstrained solutions onto our constrained trajectories differ from the distribution of unconstrained excursions only because their ‘origin’ is shifted to (δ₁, S₁); for example, see the trajectories in Fig. 1. Thus, unconstrained trajectories that happen to pass through (δ₁, S₁) are drawn from the same distribution as our constrained trajectories for S > S₁. Furthermore, for S ≫ S₁, our constrained and unconstrained excursions tend towards a mean difference of δ₁ that becomes negligible compared to the variance S − S₁, implying that our constrained excursions converge to the unconstrained solutions in this limit. Similarly, for large S₁ and sufficiently large δ₁ (meaning that the constrained branch must exceed a small mass at early times), the statistics of our constrained and unconstrained excursions precisely match. We demonstrate that our constrained halo growth histories converge to the unconstrained results in the relevant limits in Appendix  B.

3.3 Implementation in |$ {\small Galacticus}$|

We implement the algorithm described above in the open-source semi-analytic model galacticus (Benson 2012), and we make our code publicly available at https://github.com/galacticusorg/galacticus.⁷ Here, we briefly describe relevant implementation details, extensions to standard galacticus excursion set solvers, and numerical choices, and we refer the reader to the extensive online documentation for further details. Note that all results presented here are from dark matter – only merger trees; galacticus semi-analytic galaxy formation models have not been run on top of these merger trees, although the dark matter merger tree structure is unchanged when baryons are added.

In galacticus, we add functionality that marks branches of constrained merger trees according to whether they were drawn from the constrained branching rate distribution or not. ePS-based merger tree construction algorithms often analyse the ‘main’ branch of each merger tree – namely, the branch defined by selecting the more massive halo progenitor at every branching event, working backwards in time from the final halo at |$z$|₀. In practice, for constrained solutions with parameters (M₁, |$z$|₁), we define the ‘constrained’ branch of each merger tree as the branch that directly follows the constrained solution for S ≤ S₁, joined with the main branch rooted at (M₁, |$z$|₁) (the latter portion of each tree is generated using the unconstrained solution for S > S₁).

During merger tree construction, we choose sufficiently small time-steps such that neither the maximum probability for a binary merger nor the maximum fractional change in mass due to sub-resolution accretion exceed 10  per cent at any time. We also employ the time-step optimization algorithm described in Appendix A of Benson, Ludlow & Cole (2019), which speeds up our merger construction at the 10  per cent level without loss of accuracy. Throughout, we set M_min = 10⁶ M_⊙, which is well below any halo mass that directly enters the merger trees we present.

3.4 Examples of constrained merger trees and halo growth histories

Fig.4 shows examples of constrained (left) and unconstrained (right) merger trees for a final halo of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0. The constrained merger tree is generated using M₁ = 10¹²M_⊙ and |$z$|₁ = 8. This extreme constraint forces the constrained branch to grow to 10¹² M_⊙ in only ≈0.62 Gyr, corresponding to the blue point shown in Fig. 1 (however, the merger tree in Fig. 4 does not correspond to the exact excursion in Fig. 1). The constrained branch of this tree is marked in blue, while the main branch of the unconstrained tree is marked in grey. Because we only show progenitors above a mass threshold of 10¹¹M_⊙ at a discrete set of redshifts in Fig. 4, the constrained branch coincides with the main branch of the constrained tree in the visualization. In general, the constrained branch is not guaranteed to merge directly onto the main branch; we discuss this point below.

Figure 4.

Examples of constrained (left) and unconstrained (right) merger trees for a final halo of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0, which corresponds to the bottom node of each tree. The constrained achieves a mass of precisely M₁ = 10¹² M_⊙ by |$z$|₁ = 8; the progenitor at this time is marked by a star. The size of each halo is scaled by the logarithm of its mass, and the edge lengths are proportional to the logarithm of the time elapsed between successive redshifts. The constrained branch of the constrained tree is marked in blue, and the main branch of the unconstrained tree is marked in grey. For illustration, we only show haloes down to a mass of 10¹² M_⊙; progenitor haloes above this mass in the constrained (unconstrained) tree first appear at |$z$| ≈ 12 (⁠|$z$| ≈ 8).

The growth history of the constrained merger tree in Fig. 4 is clearly very different than a typical unconstrained growth history for this final halo. In particular, the constrained merger tree experiences very little late-time growth, and a significant portion of its final mass is in place by the constrained redshift of |$z$|₁ = 8. As a consequence, the progenitors of the constrained final halo above the mass resolution limit of 10¹¹ M_⊙ that we use to visualize the merger tree extend significantly farther back in time, to |$z$| ≈ 12, relative to the unconstrained progenitors above this mass that only extend to |$z$| ≈ 8. Visual inspection of many constrained and unconstrained merger trees indicates that this behaviour is typical for the extreme constraint used in Fig. 4. We explore all of the trends discussed above quantitatively in Section 4.

Fig.5 shows several growth history realizations for a final halo of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0 using the sets of constraint parameters shown in Fig. 1 and a mass resolution of 10¹⁰M_⊙. These constrained growth histories exhibit the expected behaviour, in that constraints with larger M₁ and/or larger |$z$|₁ yield growth histories that are more extreme relative to the unconstrained distribution. Decreasing either M₁ or |$z$|₁ tends to push growth histories towards the unconstrained distribution by decreasing the amount of early growth, as quantified in Section 4. Note that, the constrained branches pass precisely through their respective (M₁, |$z$|₁) constraints, as desired.

We reiterate that the main branches shown in Fig. 5 do not necessarily exceed a mass of M₁ at redshift |$z$|₁, because (unlike for the constrained merger tree shown in Fig. 4) each constrained branch does not necessarily merge directly onto its tree’s main branches. The choice of a ‘preferred’ branch for analysis is less obvious – and more application-dependent – for constrained merger trees compared to standard, unconstrained trees. In particular, the constraints we consider often increase the rate of early major mergers, for which identifying an unambiguous ‘main’ progenitor can be difficult. For example, we could define a branch based on the most massive progenitor of the parent halo at every redshift, regardless of whether this selects a connected sequence of merger events. The quantitative differences between this approach and our calculations based on standard main branches are generally small; the differences are only evident in the rare cases when main branches build up their mass at late times and experience a low-redshift major merger with the constrained branch (e.g. see the latest-forming main branches in the left-hand panel of Fig. 5).

Figure 5.

Growth histories for a final halo of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0. Five constrained solutions each are shown for (M₁, |$z$|₁) = (10¹² M_⊙, 8) (blue), (10¹² M_⊙, 6) (green), and (10¹¹M_⊙, 8) (purple); the dark (light) grey band shows 68 per cent (95 per cent) quantiles of the distribution of unconstrained main branch growth histories. The left-hand panel shows the main branches of the constrained merger trees, and the right-hand panel shows the corresponding constrained branches (see Section 3.3 for definitions). The constrained branches pass precisely through the constrained points; in our implementation, these branches are guaranteed to achieve a mass of precisely M₁ at redshift |$z$|₁. However, note that the constrained branch progenitor is not necessarily on the main branch at |$z$|₁; this implies that main branch progenitors can be less massive than M₁ at redshift |$z$|₁.

Thus, we use main branches when studying the statistical properties of our constrained solutions below. In our high-redshift galaxy application (Section 5), we explicitly track both the constrained and main branches to determine the fates of early galaxies represented by constrained branches in our merger trees.

4 RESULTS

We now describe our main results, focusing on our predictions for constrained halo growth histories in Section 4.1, formation times in Section 4.2, and merger rates in Section 4.3. Throughout this section, we only analyse the main branches of merger trees, allowing for a direct comparison between the constrained and unconstrained cases.

4.1 Halo growth histories

Fig.6 illustrates the dependence of halo growth histories on the four free parameters chosen to construct constrained solutions: the mass M₁ that the constrained branch must reach (top-left-hand panel), the constrained redshift |$z$|₁ (top-right-hand panel), the mass of the final halo M₀ (bottom-left-hand panel), and the redshift of the final halo |$z$|₀ (bottom-right-hand panel). All panels show 68 per cent quantiles of constrained halo growth history distributions, calculated using 1000 constrained realizations, along with the corresponding quantiles for a fixed set of 1000 unconstrained realizations that use M₀ = 10¹⁴ M_⊙ and |$z$|₀ = 0, for reference. All realizations use a mass resolution of 10¹⁰M_⊙ (see Appendix  B for convergence tests).

Beginning with the upper-left-hand panel of Fig. 6, we vary M₁ for a final halo of mass M₀ = 10¹⁴M_⊙ and redshift |$z$|₀ = 0, at fixed |$z$|₁ = 8. We find that increasing M₁ shifts halo growth histories systematically earlier, amplifying haloes’ phases of rapid early growth and leading to more quiescent behaviour at late times. These qualitative trends are expected based on the requirement to achieve masses of M₁ and M₀ at specific redshifts, but the detailed shape of the constrained growth histories is a non-trivial prediction of our excursion set solutions. For example, sudden major mergers onto the main branch near |$z$|₁ that do not shift the entire growth history earlier would also satisfy our constraints. However, such histories are very unlikely compared to ones that feature more overall growth than unconstrained excursions at earlier times, which can be seen from the sample of individual growth histories shown in Fig. 5. As a result, for a final halo with M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0 and for fixed |$z$|₁ = 8, the first haloes along the main branch above 10¹¹ M_⊙ (i.e. 10⁻³ times the final halo mass) appear at |$z$| ≈ 10, 9, and 8 for M₁ = 10¹², 5 × 10¹¹, and 10¹¹ M_⊙, respectively, compared to |$z$| ≈ 6.5 in the unconstrained case. We also see that, at fixed |$z$|₁, decreasing M₁ causes the constrained growth histories to approach the unconstrained distribution, as expected. We explicitly demonstrate convergence to the unconstrained solutions in this limit in Appendix  B.

Figure 6.

Growth histories for variations in constraint parameters (top panels) or final halo parameters (bottom panels). Top-left: Growth histories for final haloes of M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0 whose constrained branches reach a mass of 10¹²M_⊙ (blue), 5 × 10¹¹M_⊙ (pink), and 10¹¹ M_⊙ (purple) at |$z$|₁ = 8, listed from earliest to latest-forming. Top-right: Growth histories for final haloes of M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0 whose constrained branches reach a mass of 10¹² M_⊙ at |$z$|₁ = 8 (blue), 6 (green), and 4 (red). Bottom-left: Growth histories for final haloes at |$z$|₀ = 0 whose constrained branches reach a mass of M₁ = 10¹² M_⊙ at |$z$|₁ = 8, for final halo masses of M₀ = 10¹⁴M_⊙ (blue), 5 × 10¹³M_⊙ (brown), and 10¹³M_⊙ (red). Bottom-right: Growth histories for final haloes of M₀ = 10¹⁴M_⊙ whose constrained branches reach a mass of 10¹² M_⊙ at |$z$|₁ = 8, for final halo redshifts of |$z$|₀ = 4 (dark green), 2 (olive), and 0 (blue). In all panels, the coloured bands shows 68 per cent quantiles of the growth history distribution for 1000 constrained realizations, and the dotted grey lines indicate the corresponding quantiles for 1000 unconstrained realizations with a final halo mass of 10¹⁴ M_⊙ at |$z$|₀ = 0, for reference.

Next, in the upper-right-hand panel of Fig. 6, we vary |$z$|₁ for the same final halo properties, at fixed M₁ = 10¹² M_⊙. Increasing |$z$|₁ tends to shift growth histories systematically earlier in a coherent manner similar to the effects of increasing M₁. Quantitatively, for a final halo with M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0 and for fixed M₁ = 10¹² M_⊙, the first haloes along the main branch above 10¹¹ M_⊙ appear at |$z$| ≈ 10, 9, and 8 for |$z$|₁ = 8, 6, and 4, respectively, again compared to |$z$| ≈ 6.5 in the unconstrained case. At fixed M₁, decreasing |$z$|₁ causes the constrained solutions to approach the unconstrained distribution; again, we demonstrate convergence in this limit in Appendix  B.

For variations in both M₁ and |$z$|₁, we find that a quantitative measure of extremity based on the gradient of the Brownian bridge drift term in equation (16) predicts the magnitude of the resulting shifts in growth histories towards earlier or later times fairly well. In particular, the quantity |$\mathcal {E}\equiv (\delta _1-\delta _0)/(S_1-S_0)$| correlates with the mean temporal offset of the resulting constrained growth histories relative to the unconstrained distribution.

Finally, the bottom-left-hand and bottom-right-hand panels of Fig. 6 fix constraint parameters of M₁ = 10¹² M_⊙ and |$z$|₁ = 8, and vary M₀ and |$z$|₀, respectively. Variations in M₀ do not significantly affect the early growth, and instead simply cause the histories to flatten out (i.e. stop growing) at different redshifts. This follows because changing M₀ shifts the origin of our constrained excursions but does not affect the statistical properties of the random walks for variances larger than S₁ (corresponding to the unconstrained portion of our constrained excursions, as shown in e.g. Fig. 1). Meanwhile, increasing |$z$|₀ simply cuts off the growth histories at |$z$| > 0. As expected, we find that constrained growth histories at fixed M₁ and |$z$|₁ are greater outliers relative to the corresponding unconstrained distributions as M₀ is decreased (and approaches M₁), and as |$z$|₀ is increased (and approaches |$z$|₁).

We note that, based on the halo mass function predicted by our unconstrained ePS calculations, haloes that satisfy our high-redshift constraints are expected to descend to haloes with a specific mass distribution at |$z$| = |$z$|₀; these descendant mass distributions are derived in Appendix  C. Choosing values of M₀ well outside the expected range of descendant masses alters the resulting growth histories relative to the unconstrained distribution, effectively by making the constraint more extreme. However, even when comparing constrained and unconstrained solutions that descend to a likely halo mass for a given high-redshift constraint, the constrained growth histories statistically differ from their unconstrained counterparts. For example, this is the case for our fiducial choice of M₀ = 10¹⁴M_⊙ and |$z$|₀ = 0 given M₁ = 10¹²M_⊙, and |$z$|₁ = 8, and is reflected by a change in the overall morphology of constrained versus unconstrained merger trees, where the constrained cases form earlier (e.g. see Fig. 4).

4.2 Halo formation times

Fig.7 shows distributions of half-mass formation redshift, defined as the redshift, |$z$|_1/2, at which M(⁠|$z$| = |$z$|_1/2)/M(⁠|$z$| = 0) = 0.5, for final haloes of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0. In particular, we compare |$z$|_1/2 distributions from the sample of 1000 unconstrained merger trees described above to constrained predictions with the same M₀ and |$z$|₀, with 1000 realizations each for several values of M₁ (left-hand panel) and |$z$|₁ (right-hand panel). The unconstrained |$z$|_1/2 distributions have a mean and standard deviation of 0.61 ± 0.27 and are in reasonably good agreement with predictions from cosmological simulations (e.g. Wechsler et al. 2002; Nadler et al. 2022).

Figure 7.

Cumulative distributions of formation time, defined as the redshift at which a halo reaches half of its |$z$| = 0 mass, for final haloes of mass M₀ = 10¹⁴ M_⊙ at redshift |$z$|₀ = 0. Constrained realizations for various Brownian bridge parameters are shown by coloured lines, and the cumulative distribution for standard unconstrained excursions is shown by the grey line; each distribution is computed using 1000 constrained or unconstrained realizations. Left: Half-mass formation redshift distributions for constrained merger trees whose constrained branches reach a mass of 10¹² M_⊙ (blue), 5 × 10¹¹ M_⊙ (pink), and 10¹¹ M_⊙ (purple) at |$z$|₁ = 8, listed from earliest to latest-forming. Right: Same as the left-hand panel, for constrained merger trees whose constrained branches reach a mass of 10¹² M_⊙ at |$z$|₁ = 8 (blue), 6 (green), and 4 (red).

The constrained |$z$|_1/2 distributions behave as expected based on the discussion in the previous subsection. In particular, at fixed |$z$|₁ = 8, these distributions have a mean and standard deviation of 1.9 ± 0.9, 1.6 ± 0.8, and 1.1 ± 0.6 for M₁ = 10¹², 5 × 10¹¹, and 10¹¹ M_⊙, respectively. These shifts towards earlier formation times result from more rapid early accretion, which is necessary to achieve the constrained mass. Quantitatively, these constrained realizations are thus ≈4σ, 3σ, and 2σ outliers relative to the unconstrained formation time distribution, respectively.

Meanwhile, at fixed M₁ = 10¹² M_⊙, the constrained formation time distributions have a mean and standard deviation of 1.9 ± 0.9, 1.5 ± 0.8, and 1.1 ± 0.5 for |$z$|₁ = 8, 6, and 4, again making them ≈4σ, ≈3σ, and 2σ outliers in the mean relative to the unconstrained formation time distribution, respectively. Moreover, our constrained solutions produce individual merger trees with formation times that are extreme outliers relative to the unconstrained distribution; for example, the realizations in our |$z$|₁ = 8, M₁ = 10¹² M_⊙ case with |$z$|_1/2 = 4 that occur about one in a hundred times arise by chance in the corresponding unconstrained trees roughly one in ten billion times, representing a ∼10⁸ times speed up.

4.3 Halo merger rates

Having shown that halo growth histories are shifted systematically earlier in the presence of constraints that guarantee the existence of a massive progenitor at early times, we now investigate the distributions of mergers that build up the final haloes in constrained scenarios. In particular, we will demonstrate that high-redshift constraints lead to systematically earlier merging activity and an enhanced rate of minor mergers after the constraint is achieved.

Fig.8 shows the distribution of merger events along the main branch as a function of redshift, calculated as ΔM/M, where M is the mass of the main branch progenitor at a given redshift and ΔM is the amount by which M increases during a given merger. For clarity, we only show mergers that occur while the main branch has achieved at least 10 per cent of its final mass, corresponding to a mass of 10¹³ M_⊙ for the M₀ = 10¹⁴M_⊙, |$z$|₀ = 0 haloes used to generate the results shown in Fig. 8. Because the resolution of these merger trees is M_res = 10¹⁰M_⊙, we can only measure merger events with ΔM/M > 10⁻⁴ for the solutions discussed here.

Figure 8.

The mass and redshift distribution of mergers along the main branch for 1000 realizations of constrained merger trees with M₁ = 10¹² M_⊙ (blue) and 10¹¹ M_⊙ (purple), both with |$z$|₁ = 8, for a final halo of mass M₀ = 10¹⁴ M_⊙ at redshift |$z$|₀ = 0. Lines show four iso-proportion contours of the joint probability density, and the corresponding unconstrained distribution is shown as a grey histogram; the top and right-hand panels show the corresponding marginal distributions, where the spread encompasses the four iso-proportion contours. Only mergers that take place when the main branch has achieved at least 10 per cent of its final mass (i.e. M > 10¹³M_⊙) are considered; there are typically a few thousand such mergers above our resolution limit of ΔM/M > 10⁻⁴.

Fig.8 shows all individual mergers from 1000 realizations of constrained merger trees for M₁ = 10¹² and 10¹¹ M_⊙, both using |$z$|₁ = 8, and we compare these to the underlying distribution for the same number of realizations of unconstrained merger trees. We find that merger times are shifted systematically earlier as M₁ increases; even for M₁ = 10¹¹ M_⊙, for which the constrained and unconstrained halo growth history 68 per cent quantiles partially overlap (see the top-left-hand panel of Fig. 6), this shift is significant. For example, the mean and standard deviation of the redshift at which any merger with ΔM/M > 10⁻² occurs is |$z$| = 1.1 ± 0.6, |$z$| = 1.7 ± 0.9, and |$z$| = 2.8 ± 1.7 for the unconstrained, M₁ = 10¹¹ M_⊙, and M₁ = 10¹² M_⊙ case, respectively. These shifts are driven by the requirement that the constrained branch progenitor achieves an unusually large mass at early times and then (eventually) merges onto the main branch. To guarantee conservation of mass, there are then fewer late-time mergers in the constrained trees, and the late-time mergers that do occur are biased towards lower mass ratios relative to the unconstrained distribution.

We find that there is a slightly enhanced tail of small mergers in the constrained cases that compensates for major merger events with the constrained branches at early times. None the less, the overall mass distribution of mergers along the main branch does not significantly differ among the cases shown in Fig. 8, despite their different distributions of merger redshifts. The overall similarity in merger mass distributions is reasonable because most merger events are minor, and the distribution of first-upcrossings of the ePS barrier after shifting the origin of our excursions appropriately, for each constrained case, is not significantly influenced by the constraints.

5 THE DESCENDANTS OF THE FIRST GALAXIES

We now illustrate their utility in a specific astrophysical context. In particular, we consider the extremely luminous high-redshift galaxy candidates detected by JWST. Because of observational selection effects, these systems represent the most luminous and massive objects in existence at early times, and are therefore outliers relative to the underlying high-redshift galaxy and halo populations. In Section 5.1, we estimate the halo masses of these high-redshift systems; we interpret these masses, and the corresponding redshifts, as the M₁ and |$z$|₁ parameters entering our Brownian bridge excursions. In Section 5.2, we then use our constrained merger tree algorithm to generate the most likely distributions of growth histories for these high-redshift systems. Specifically, we calculate the distribution of |$z$|₀ = 0 descendant halo masses given each JWST galaxy candidate’s redshift and inferred halo mass, and we draw from these distributions to set M₀ in our constrained merger tree realizations.

5.1 Estimating high-z JWST Galaxy candidates’ halo masses

We consider the seven high-redshift JWST galaxy candidates reported by Labbe et al. (2022), which were selected based on the presence of Balmer breaks. We emphasize that these objects’ redshifts are not spectroscopically confirmed, and that several observational and modelling systematics may impact the interpretation of their stellar masses and redshifts; thus, we refer to them as galaxy candidates.⁸ Note that, we do not include the even more extreme candidate reported by Naidu et al. (2022a) in our calculations because its high-redshift solution is under scrutiny (e.g. Naidu et al. 2022b).

Our procedure for estimating halo masses is based on the extreme value statistics (EVS) calculations and code from Lovell et al. (2022).⁹ In particular, these authors present distributions of the maximum halo and stellar masses predicted to exist in a fiducial ΛCDM cosmology as a function of redshift. We calculate the quantile of the stellar mass EVS distribution, at the appropriate redshift, corresponding to each of the Labbe et al. (2022) galaxy candidates; these systems’ reported redshifts and stellar masses span 7.5 ≲ |$z$| ≲ 11 and 10¹⁰ M_⊙ ≲ M_* ≲ 10¹¹M_⊙, respectively. We confirm that the quantile of the least (most) extreme galaxy candidate relative to the stellar mass EVS distribution corresponds to a ≈2σ (4σ) outlier assuming Gaussian statistics, as shown by Lovell et al. (2022).

We then translate these statistics to halo masses by finding the halo mass that minimizes the difference between each galaxy’s stellar mass and halo mass EVS quantile, at the appropriate redshift. This yields a distribution of halo masses that spans 2.6 × 10¹¹M_⊙ ≲ M_halo ≲ 3.5 × 10¹² M_⊙.¹⁰ Specifically, for the seven Labbe et al. (2022) galaxy candidates with estimated redshifts and stellar masses of (⁠|$z$|_*, log (M_*/M_⊙) = (9.9, 11.0), (9.9,10.7), (10.8,10.4), (7.6,11.1), (8.9,10.8), (9.3,10.4), (8.5,10.2), we estimate halo masses of log (M_halo(⁠|$z$|_*)/M_⊙) = 12.5, 12.2, 11.9, 12.2, 11.9, 11.6, 11.4, respectively. We always order these galaxies in decreasing order of their average expected |$z$|₀ = 0 halo masses derived from our calculation in Appendix  C. In the following subsection, we use each galaxy candidate’s redshift and estimated halo mass as constraint parameters to generate a family of merger trees for the Labbe et al. (2022) sample.

Before moving on, we discuss the basic assumption underlying our halo mass estimates and areas for improvement. In essence, we have assumed that the haloes of JWST galaxy candidates are equal outliers in halo mass, relative to the underlying high-redshift halo population, as the galaxies themselves are in stellar mass, relative to the underlying high-redshift galaxy population. Although this assumption is plausible, it is almost certainly incorrect in detail. Empirically, high-redshift galaxy and halo populations are related via a non-linear galaxy–halo connection with potentially significant scatter, which implies that the one-to-one correspondence our calculation assumes likely breaks down (e.g. see Wechsler & Tinker 2018 for a review). More sophisticated abundance-matching techniques and empirical modelling based on halo populations in cosmological simulations (e.g. Behroozi & Silk 2015, 2018; Behroozi et al. 2020), or direct estimates of halo properties from hydrodynamic simulations (e.g. Lovell et al. 2022), would thus improve the accuracy of our estimates. Physically, correlations between secondary halo properties (e.g. accretion rate) and galaxy properties that correlate with observable quantities (e.g. star formation rate) imply that a one-dimensional correspondence based on halo and stellar mass is likely inadequate for explaining why some high-redshift haloes host extremely luminous galaxies (e.g. Mason et al. 2023; Mirocha & Furlanetto 2022). We therefore proceed with the caveat that our halo mass estimates are subject to improvement. The utility of our constrained merger tree realizations for predicting the growth histories of high-redshift systems is independent of these uncertainties.

5.2 Constrained excursion set predictions

We integrate our halo mass estimates for the high-redshift JWST galaxy candidates into constrained merger tree realizations as follows. For each galaxy from Labbe et al. (2022), we use our estimated halo mass and its reported redshift to set M₁ and |$z$|₁, respectively. We then calculate the most likely distribution of descendant halo masses, given M₁ and |$z$|₁, assuming the halo mass function at |$z$|₀ = 0. We derive this descendant mass distribution in Appendix  C, finding that the most likely |$z$| = 0 descendant masses for the seven galaxies listed above are log (M_halo(⁠|$z$|₀ = 0)/M_⊙) = 14.4, 14.2, 14.0, 13.9, 13.8, 13.5, and 13.3, respectively. We use the full descendant mass distributions derived in Appendix  C to draw values of M₀ for constrained merger tree realizations; thus, unlike the results presented up to this point, we now use a distribution of M₀ for each set of M₁ and |$z$|₁ (at fixed |$z$|₀ = 0) rather than a single value. Note that, because the JWST galaxy candidates are ≈2σ to 4σ outliers relative to the underlying halo population, our constrained merger tree realizations are on average ∼20 to 10⁴ times more efficient than a calculation based on down-sampling unconstrained merger trees.

We generate 1000 realizations of constrained merger trees for each of the seven Labbe et al. (2022) JWST galaxy candidates in this manner, and Fig. 9 shows the distribution of resulting constrained and main branch growth histories for each system. These growth histories behave similarly to the M₁ = 10¹² M_⊙, |$z$|₁ = 8 constraint with M₀ = 10¹⁴M_⊙ studied above, and typically reach halo masses of 10¹¹ M_⊙ (i.e. 10⁻³ times the final halo mass) by |$z$| ≈ 10. However, the spread of M₀ values predicted by our ePS descendant calculation for each JWST galaxy candidate yields a larger diversity of growth histories than our calculations with a single value of M₀; considering the Labbe et al. (2022) galaxy population as a whole, the spread of constraint parameters derived from the estimated galaxy properties further increases this scatter.

Figure 9.

Growth histories for the seven high-redshift JWST galaxy candidates from Labbe et al. (2022). Bands shows 68 per cent quantiles of the distribution of unconstrained main and constrained branch growth histories. The left-hand panel shows the main branches of the constrained merger trees, and the right-hand panel shows the constrained branches (see Section 3.3). Each constrained halo mass is derived from the corresponding galaxy’s estimated stellar mass and photometric redshift, as described in Section 5.1, and the constrained redshift is taken to be equal to the photometric redshift.

We predict a distribution of half-mass formation times for the main branches of the Labbe et al. (2022) galaxy candidates with a mean and standard deviation of |$z$|_1/2 = 2.4 ± 1.1, averaged over the seven systems; this result does not significantly change when we use the constrained branches instead. We also calculate the time, after each system’s observation epoch, until the first major merger (defined as any merger with a mass ratio above 1 to 3) along the constrained branch that represents each galaxy. For the seven galaxies listed above, we predict that major mergers occur t_merger = 0.8, 1.4, 2.0, 1.8, 2.5, 3.3, and 4.2 Gyr after the corresponding observation epochs, on average. These merger times systematically decrease as the typical final halo mass decreases because the constrained branches then merge onto (or are identical to) lower mass main branches.

Finally, we predict the fates of the high-redshift JWST galaxy candidates by calculating the fraction of these systems predicted to descend to satellites of the main-branch halo in their merger trees. In the context of our constrained merger tree realizations, these galaxies are predicted to remain isolated only if the constrained branch and main branches coincide for all |$z$| < |$z$|₁; otherwise, the constrained branch must eventually merge onto a larger halo and becomes a satellite of that system.¹¹ We find that the probability for each Labbe et al. (2022) galaxy candidate to coincide with its main branch, and thus remain a central galaxy at |$z$| = 0, is 33, 17, 17, 15, 14, 7, and 5 per cent, respectively. Interestingly, these average isolated fractions decrease as the mean expected |$z$|₀ = 0 halo mass decreases. For each galaxy candidate, these isolated fractions are fairly strong functions of M₀ and typically drop by a factor of ∼2 when only the upper half of each galaxy candidate’s M₀ distribution is considered, because the constrained halo is more likely to merge onto a more massive branch to achieve the final mass in such cases.

Our calculations imply that the JWST galaxy candidates likely become satellites of a larger system at |$z$| ≈ 2, although the most extreme high-redshift systems have a non-negligible chance to remain isolated and inhabit massive (≳ 10¹⁴ M_⊙) haloes at |$z$| = 0. Both evolutionary pathways are interesting observationally; mergers at |$z$| ≈ 2 coincide with the peak of the cosmic star formation rate density and may help quench massive clusters at these redshifts. Meanwhile, the fraction of extreme high-redshift systems predicted to remain isolated may be remarkably luminous today, if their star formation continues to be fueled and in the absence of feedback or merger-driven quenching; on the other hand, if they occupy isolated environments, they may quickly deplete their gas and appear as ‘fossils’ of the first galaxies today. As discussed below, we plan to study these constrained systems’ baryonic components in future work.

6 THE MILKY WAY–LARGE MAGELLANIC CLOUD MERGER

The constrained branches generated using our Brownian bridge excursions do not necessarily coincide with or merge directly onto the main branches of our constrained merger trees. In particular, our constrained solutions only guarantee that the constrained branch reaches a mass of precisely M₁ at redshift |$z$|₁; its position in the overall merger tree is a non-trivial prediction of our model. For example, the constrained halo can be a lower-mass progenitor along the main branch, or it can merge onto a secondary branch that eventually merges with the main branch. Meanwhile, observational constraints on systems’ accretion histories are often easiest to interpret as mergers with the main branch progenitor of the main halo. In future work, we therefore plan to consider more general classes of excursion set solutions, and to experiment with directly modifying merger trees during the construction process (e.g. by forcing mergers onto the main branch with desired properties to occur) in order to produce constrained growth histories that match specific sets of observational constraints in detail.

Here, we briefly illustrate that even our constrained Brownian bridge excursions can consistently yield mergers of the constrained halo directly onto the main branch, using the merger between the Milky Way and the Large Magellanic Cloud (LMC) as an example. The LMC is the Milky Way’s largest satellite galaxy, with a stellar mass of 3 × 10⁹ M_⊙ (van der Marel et al. 2002) and an estimated halo mass of ≈1 to 3 × 10¹¹M_⊙ (Erkal et al. 2019; Erkal & Belokurov 2020; Shipp et al. 2021). Kinematic measurements imply that the LMC fell into the Milky Way recently, within the last ≈2 Gyr (Kallivayalil et al. 2006, 2013). This makes the Milky Way–LMC merger fairly rare, in the sense that only |$\approx 10~{{\ \rm per\ cent}}$| of Milky Way-mass haloes host such a massive subhalo today; an even smaller fraction of these massive subhaloes accrete at very late times and reside near the host centre today like the LMC (e.g. Liu et al. 2011; Nadler et al. 2021).

The Milky Way–LMC merger provides a particularly clean example of a system for which the ‘constrained’ branch (i.e. the LMC) should merge directly onto the ‘main’ branch (i.e. the Milky Way progenitor at |$z$| ≈ 0.1, corresponding to a lookback time ≈1 Gyr). In particular, it is observationally disfavoured (and exceedingly unlikely in ePS theory) for a constrained branch that resembles the LMC by reaching a halo mass of M₁ ≈ 3 × 10¹¹ M_⊙ at |$z$|₁ ≈ 0.1 to merge with a non-main branch progenitor of the Milky Way. This would require another halo even more massive than the LMC to merge with the main branch of the Milky Way host halo at late times, which has a total mass of M₀ ≈ 10¹²M_⊙ at |$z$|₀ = 0 (e.g. see Bland-Hawthorn & Gerhard 2016 for a review). Thus, we assume that the constrained branch should represent the LMC.

In Fig. 10, we show that imposing a constraint of M₁ = 3 × 10¹¹ M_⊙ at redshift |$z$|₁ = 0.1 for a final halo of mass M₀ = 10¹²M_⊙ at |$z$|₀ = 0 consistently yields direct mergers of the constrained and main branches that resemble the infall of the LMC into the Milky Way. Interestingly, the main branches of Milky Way progenitors that experience late-time mergers with LMC analogues are biased towards later-forming growth histories relative to the unconstrained distribution at the ≈1σ to 2σ level. The ability to generate such growth histories without down-sampling represents a first step towards a controlled semi-analytic model of the Milky Way’s entire accretion history, including all of its most significant merger events. We intend to pursue this application of our constrained excursion set solutions in a future study.

Figure 10.

Growth histories for a final halo of mass M₀ = 10¹²M_⊙ at |$z$|₀ = 0, similar to the mass of the Milky Way halo. Dark (light) blue bands show 68 per cent (95 per cent) quantiles of the distribution of constrained main branch growth histories for a constraint of (M₁, |$z$|₁) = (3 × 10¹¹ M_⊙, 0.1), chosen such that the constrained branch resembles the Large Magellanic Cloud and merges onto the main branch. Five individual constrained growth histories are shown for the main branch (solid lines) and constrained branch (i.e. the Large Magellanic Cloud analogue; dashed lines). Dark (light) grey bands show 68 per cent (95 per cent) quantiles of the distribution of unconstrained main branch growth histories for the same final halo mass and redshift.

7 POTENTIAL CAVEATS

We have shown that our constrained merger trees satisfy the relevant Brownian bridge constraints, and we have described their qualitative and quantitative behaviour. However, we have not attempted to calibrate our predictions for merger rates and halo growth histories to analogous quantities measured in cosmological simulations. Importantly, we do not use the first-crossing distributions derived below to compute (conditional) halo mass functions. Instead, we only use the modified first-crossing rates to construct merger trees, as described in Section 3. A comparison between the (conditional) population statistics of haloes generated by our constrained solutions and cosmological simulations is beyond the scope of this study and represents an interesting avenue for future work.

There are multiple fronts along which calibration of merger rates and halo growth histories should proceed. First, because we modified standard ePS merger tree construction algorithms as described in Section 3, and particularly because we sampled progenitor masses from the full interval [M_res, M − M_res], even our predictions for unconstrained growth histories may require some recalibration. Concretely, such recalibration could follow the procedure in Benson (2017) by optimizing the parameters of our generalized modifier function (equations 22–25) to match merger rates and mass accretion histories measured in cosmological N-body simulations; we intend to address this in future work.

Second, it will be necessary to calibrate our predictions for constrained growth histories to simulation results. Rather than calibrating to mass accretion histories and merger rates for the full distribution of haloes (e.g. in some interval of final halo mass M₀ and redshift |$z$|₀), this will require calibration to the subset of simulated systems that satisfy a particular constraint (e.g. those with a progenitor that exceeds a mass M₁ by redshift |$z$|₁), performed for many sets of constraint parameters. In practice, this down-sampling procedure may yield very few suitable growth histories in simulations, depending on the simulation size and the extremity of the constraint; for example, haloes with masses comparable to our high-redshift JWST estimates are rare even in large-volume cosmological simulations (e.g. see Lovell et al. 2022). Techniques for ‘genetically modifying’ initial conditions to alter the accretion histories of particular haloes found in simulations (e.g. Porciani 2016; Roth, Pontzen & Peiris 2016; Rey & Pontzen 2018; Stopyra et al. 2021) may aid such comparisons. We anticipate that comparing the correlated nature of halo growth predicted by simulations and by our method in the presence of extreme constraints will be fruitful and may yield new insights into the impact of large-scale environment on smaller scale merger events. We intend to explore these questions in future work.

Finally, we note that, because our implementation of constrained Brownian bridge excursions assumes a constant barrier as a function of ΔS, generalizing our approach for non-constant barriers (e.g. for the warm dark matter barriers considered in Benson et al. 2013) will require dedicated study.

8 DISCUSSION AND CONCLUSIONS

We have presented a new class of constrained excursion set solutions to model halo growth histories. Specifically, we used Brownian bridges to model constrained excursions of the linear overdensity field in ePS theory. These constrained excursions are guaranteed to pass through a certain overdensity and mass variance, thereby fixing a point in the final halo’s growth history, and to accurately represent the distribution of unconstrained excursions that would have satisfied the desired constraint. We implemented these excursions in a modified merger tree construction algorithm to guarantee that a progenitor halo with a specific mass exists at a desired redshift, regardless of how rare the resulting growth histories are relative to the underlying unconstrained distribution.

Our main results focus on merger trees that satisfy extreme high-redshift constraints, and thus represent the tail of the ΛCDM halo population with the most rapid early growth. Specifically, we have quantified how halo growth histories (Fig. 6), formation time distributions (Fig. 7), and merger rates (Fig. 8) depend on the constraint parameters. As a proof of principle, we applied our method to predict the most probable halo growth histories of seven high-redshift galaxy candidates identified by Labbe et al. (2022) in early JWST data (Fig. 9). These realizations are at least thousands of times more efficient than calculations based on the subset of standard unconstrained merger trees that satisfy the relevant constraints. We predict that the descendants of high-redshift JWST galaxy candidates most likely merge with larger systems ≈2 Gyr after the observation epoch and descend to haloes with masses of ≈10¹⁴ M_⊙ today. We will explore whether these systems survive as distinct satellites or merge with the central galaxy in future modelling work that combines semi-analytic galaxy formation models implemented in galacticus with our constrained solutions.

We have not yet explored how constrained growth histories translate to the baryonic components of haloes, e.g. stellar mass, gas mass, and central black hole mass, and the other properties of the galaxies they host. There are significant complexities and uncertainties associated with semi-analytically modelling galaxy formation using ePS merger trees (e.g. see Benson 2010 for a review), and we plan to investigate how such astrophysical modelling interplays with our constrained solutions in future studies. Here, we describe some questions that will motivate this work.

First, are the galaxies that occupy rapidly growing haloes that satisfy high-redshift equally significant outliers, relative to the distribution of unconstrained galaxies that exist at the corresponding redshifts, as their dark matter haloes? Naively, the eventual merger of the constrained branch onto the main branch implies that the main halo’s accretion rate is increased as a result of the constraint, leading to increased star formation rates and overly luminous central galaxies, compared to unconstrained systems at the same redshift. Several effects may complicate this picture in detail. For example, as we have shown, the galaxy occupying the constrained-branch halo often becomes a satellite of the main-branch progenitor. Depending on the merger time-scale, these constrained-branch satellites can survive as distinct galaxies or merge onto the central; in the case of a major merger between the constrained and main-branch galaxies, star formation can either rapidly cease or be reignited (e.g. see Somerville & Davé 2015 for a review). Second, how do (potentially extreme) variations in haloes’ growth histories affect the star formation histories of the galaxies they host? As an illustrative example, consider the scatter in the stellar mass–halo mass relation for ultra-faint dwarf galaxies, most of which are thought to be quenched by reionization (e.g. see Simon 2019 for a review). Controlled simulations indicate that this scatter is set by how quickly stars formed before reionization. In particular, for a fixed final halo mass, haloes that formed more of their stars before reionization are expected to host more luminous dwarf galaxies (Rey et al. 2019); this effect may contribute to the large stellar mass–halo mass relation scatter predicted by hydrodynamic simulations at the faint end (e.g. Munshi et al. 2021). In turn, these upscattered systems are expected to accrete more rapidly at early times and can thus be sampled in a controlled manner using our Brownian bridge technique. We can therefore leverage our method – which enables us to generate thousands of constrained merger trees at the same computational expense as a non-cosmological simulation of comparable resolution – to explore the relationship between the early growth of low-mass haloes and the resulting properties of the ultra-faint galaxies they host.

Finally, we note that supermassive black holes (SMBHs) have been detected at redshifts |$z$| ≳ 6 in optical and near-infrared data (e.g. Onoue et al. 2019; Yang et al. 2020), and these systems are also likely to occupy rare, rapidly growing haloes (e.g. see Inayoshi, Visbal & Haiman 2020 for a review). Our method can thus be used to study potential connections between early galaxy and SMBH populations in the presence of appropriate seeding mechanisms (e.g. see Volonteri, Habouzit & Colpi 2023 for a recent study).

Our constrained excursion set solutions may in principle be generalized to model the growth histories of astrophysical systems subject to multiple constraints, rather than enforcing a single constraint on the constrained progenitor’s mass and redshift. For example, we have already demonstrated that certain regions of our Brownian bridge parameter space consistently yield constrained branches that directly merge onto the main branch, resembling the merger between the Milky Way and Large Magellanic Cloud (Fig. 10); we expect that similar techniques can likely reproduce any series of major mergers in a halo’s growth history. Enforcing multiple Brownian bridge constraints along a single excursion and implementing other classes of constrained random walks will allow desired growth histories to be reconstructed in increasing detail. In this context, it is interesting to consider the extent to which the detailed growth history of a system like the Milky Way encodes the precise excursion of the linear overdensity field from which it originated.

Acknowledgement

We thank Fangzhou Jiang for helpful discussions. Computing resources used in this work were made available by a generous grant from the Ahmanson Foundation. VG amd TD acknowledge the support from the National Science Foundation under grant no. PHY-2013951 and from NASA through the Astrophysics Theory Program, Award Number 21-ATP21-0135.

DATA AVAILABILITY

Our constrained merger tree code is publicly available at https://github.com/galacticusorg/galacticus.

Footnotes

1

These assumptions imply that standard ePS excursions are Wiener processes in the continuum limit; thus, the constrained solutions we present can be applied to any stochastic process y(x) with a suitable barrier.

2

We write M[S] as a functional because of the integral relation between these variables in equation (1).

3

In particular, every branching event adds two progenitors in the Cole et al. (2000) approach: one with mass M′ drawn from the progenitor mass distribution function, and the other with mass M − M′ to ensure mass conservation. The overall branching rate must be reduced by a factor 2 to counteract this. In the original Cole et al. (2000) method, this factor of 1/2 is implicitly included because the integral spans only half of the full mass range.

4

As discussed in Neistein & Dekel (2008), binary merger algorithms only approximate the true ePS progenitor mass distribution, which generally includes more than two progenitors per branching event. The binary Parkinson et al. (2008) algorithm and our modifications thereof have none the less been calibrated to match simulations well, and also predict progenitor mass functions and branching rates that are fairly self-consistent (e.g. see Jiang & van den Bosch 2014).

5

Russell (2014) used a similar technique to calculate void statistics in the ‘two-barrier’ excursion set formalism, following arguments from Marchal (2003); also see Sheth & van de Weygaert (2004); Russell (2013).

6

Such excursions are formally equivalent to Brownian meanders, i.e. Brownian bridges with a supremum of δ₁ at S₁ (also see Russell 2014).

7

Specifically, we use commit 5ea7873f6a3a1ef08f67d114710b77cca7d9529b.

8

As this paper was being prepared, a distinct sample of JWST galaxy candidates were spectroscopically confirmed (Curtis-Lake et al. 2022; Robertson et al. 2022); furthermore, Bouwens et al. (2022) reported that the Labbe et al. (2022) candidates’ photometrically derived properties are relatively robust based on consistency between independent analyses of these systems

9

See https://github.com/christopherlovell/evstats.

10

The halo mass function used to calculate the EVS distribution uses the M₂₀₀ mass definition (Lovell et al. 2022); we do not attempt to convert these masses to our ePS halo mass definition because this |$\mathcal {O}(10~{{\ \rm per\ cent}})$| difference is small compared to other uncertainties in our estimate.

11

In turn, this requirement implies that the constrained and main branches coincide at all redshifts, because we have defined the constrained branch as the main branch of the tree rooted at (M₁, |$z$|₁) for |$z$| > |$z$|₁.

12

We note that convergence properties are more subtle as a function of δ₁. As δ₁ → δ₀, the Brownian bridge drift term (equation (16)) vanishes, but the covariance (equation 18) differs from the unconstrained case. This follows because unconstrained trajectories may not pass through (S₁, δ₀) unless S₁ is sufficiently large. Thus, although the growth histories in the top-right-hand panel of Fig. 6 approach the unconstrained distribution as |$z$|₁ decreases, convergence as a function of |$z$|₁ (or δ₁) requires a sufficiently small M₁ (or large S₁).

APPENDIX A: INTEGRAL EQUATIONS FOR THE FIRST-CROSSING DISTRIBUTION

A1 Unconstrained excursions

For unconstrained excursions, the integral equation that determines f(S) follows from plugging equations (11) and (14) into equation (10), giving

$$ \begin{eqnarray*} 1 &=& \int _0^S f(\hat{{ S}})\mathrm{d}\hat{{ S}} + \int _{-\infty }^{\tilde{B}(S)} \left[\frac{1}{\sqrt{2\pi S}}\exp \left(-\frac{\delta ^2}{2S}\right)-\int _{0}^S f(\hat{{ S}})\frac{1}{\sqrt{2\pi (S-\hat{{ S}})}}\exp \left(-\frac{(\delta -\tilde{B}(\hat{S}))^2}{2(S-\hat{{ S}})}\right)\mathrm{d}\hat{{S}}\right]\mathrm{d}\delta \nonumber \\ &=& \int _0^S f(\hat{{ S}})\mathrm{d}\hat{{ S}} + \mathrm{erf}\left[\frac{\tilde{B}(S)}{\sqrt{2S}}\right] - \int _0^S f(\hat{{ S}}) \mathrm{erf}\left[\frac{\tilde{B}(S)-\tilde{B}(\hat{{ S}})}{\sqrt{2(S-\hat{S})}}\right]\mathrm{d}\hat{{ S}}. \end{eqnarray*} $$

(A1)

For a constant barrier |$\tilde{B}(S)\equiv B$| (relevant, for example, in ePS predictions for cold dark matter cosmologies), the final term vanishes because |$\tilde{B}(S)-\tilde{B}(\hat{\boldsymbol{ S}})=0$|⁠. Equation A1 can then be solved by differentiating with respect to S, which yields the standard first-crossing distribution

$$ \begin{eqnarray*} f_{\mathrm{constant\ barrier}}(S)=\frac{B}{\sqrt{2\pi S^3}}\exp \left(-\frac{B^2}{2S}\right). \end{eqnarray*} $$

(A2)

A2 Constrained Brownian bridge excursions

For our constrained Brownian bridge excursions, the the integral equation that determines f(S) follows from plugging equations (19) and (20) into equation (10). To write this explicitly, we first expand equation (19) and equation (20) according to

$$ \begin{eqnarray*} P^{\prime }_0(\delta ,S) = \mathcal {N}(\Delta \delta ^{\prime },\Delta S^{\prime }) \nonumber \\ &=& \mathcal {N}\left(\delta -\tilde{B}(\hat{{ S}})-\frac{S-S_0}{S_1-S_0}(\delta _1-\delta _0),\frac{(S_1-S)(S-S_0)}{S_1-S_0}\right) \nonumber \\ &=& \frac{1}{\sqrt{2\pi }}\sqrt{\frac{(S_1-S_0)}{ (S_1-S)(S-S_0)}}\exp \left(-\frac{\left[\delta -\tilde{B}(\hat{{ S}})-\frac{S-S_0}{S_1-S_0}(\delta _1-\delta _0)\right]^2}{2\frac{(S_1-S)(S-S_0)}{(S_1-S_0)}}\right), \end{eqnarray*} $$

(A3)

$$ \begin{eqnarray*} P^{\prime }_0[\delta ,S|\tilde{B}(\hat{\boldsymbol{ S}}),\hat{\boldsymbol{ S}}] &=& \mathcal {N}(\Delta \delta ^{\prime },\Delta S^{\prime }) \nonumber \\ &=& \mathcal {N}\left(\delta -\tilde{B}(\hat{\boldsymbol{ S}})-\frac{S-\hat{\boldsymbol{ S}}}{S_1-S_0}(\delta _1-\delta _0),\frac{(S_1-S)(S-\hat{\boldsymbol{ S}})}{S_1-S_0}\right) \nonumber \\ &=& \frac{1}{\sqrt{2\pi }}\sqrt{\frac{(S_1-S_0)}{ (S_1-S)(S-\hat{\boldsymbol{ S}})}}\exp \left(-\frac{\left[\delta -\tilde{B}(\hat{\boldsymbol{ S}})-\frac{S-\hat{\boldsymbol{ S}}}{S_1-S_0}(\delta _1-\delta _0)\right]^2}{2\frac{(S_1-S)(S-\hat{\boldsymbol{ S}})}{(S_1-S_0)}}\right). \end{eqnarray*} $$

(A4)

Finally, we plug these expressions into equation (10) to obtain

$$ \begin{eqnarray*} 1 &=& \int _0^S f(\hat{{ S}})\mathrm{d}\hat{S} + \int _{-\infty }^{\tilde{B}(S)} \Bigg [\frac{1}{\sqrt{2\pi }}\sqrt{\frac{(S_1-S_0)}{ (S_1-S)(S-S_0)}}\exp \left(-\frac{\left[\delta -\tilde{B}(\hat{{ S}})-\frac{S-S_0}{S_1-S_0}(\delta _1-\delta _0)\right]^2}{2\frac{(S_1-S)(S-S_0)}{(S_1-S_0)}}\right) \nonumber \\ &&-\,\int _{0}^S f(\hat{{ S}})\frac{1}{\sqrt{2\pi }}\sqrt{\frac{(S_1-S_0)}{ (S_1-S)(S-\hat{{ S}})}}\exp \left(-\frac{\left[\delta -\tilde{B}(\hat{{ S}})-\frac{S-\hat{{ S}}}{S_1-S_0}(\delta _1-\delta _0)\right]^2}{2\frac{(S_1-S)(S-\hat{{ S}})}{(S_1-S_0)}}\right)\mathrm{d}\hat{{ S}}\Bigg ]\mathrm{d}\delta \nonumber \\ &=& \int _0^S f(\hat{{ S}})\mathrm{d}\hat{{ S}} + \mathrm{erf}\left[\frac{\tilde{B}(S)-\tilde{B}(\hat{{ S}})-\frac{S-S_0}{S_1-S_0}(\delta _1-\delta _0)}{\sqrt{2\frac{(S_1-S)(S-\hat{{ S}})}{(S_1-S_0)}}}\right] - \int _0^S f(\hat{{ S}}) \mathrm{erf}\left[\frac{\tilde{B}(S)-\tilde{B}(\hat{{ S}})-\frac{S-\hat{{ S}}}{S_1-S_0}(\delta _1-\delta _0)}{\sqrt{2\frac{(S_1-S)(S-\hat{{ S}})}{(S_1-S_0)}}}\right]\mathrm{d}\hat{{ S}}. \end{eqnarray*} $$

(A5)

The arguments of the error functions in equation (A5) simplify for a constant barrier, but even in this case there is no analytic solution for f(S) in our constrained model because the argument of the error function inside the last integral depends on S.

A3 Numerical Solution Technique

To solve equations (A1) and (A5) for the unconstrained and constrained first-crossing distributions, respectively, we follow Benson et al. (2013) by discretizing variance into N logarithmically spaced intervals of width ΔS. We then use the midpoint integration method described in Du et al. (2017) to discretize the integrals; Du et al. (2017) show that this procedure outperforms previous solutions based on trapezoid integration. The solution is iterative: at each step, f(S_i) only depends on all f(S_j) (where j < i), and on the barrier. We implement this solution generally, such that the first-crossing distribution can be obtained for any (non-constant) barrier in both unconstrained and constrained cases.

APPENDIX B: CONVERGENCE TESTS OF CONSTRAINED SOLUTIONS

B1 Convergence to unconstrained solutions

As argued in Section 2.3 and indicated in the top-left-hand panel of Fig. 6, our constrained halo growth histories should converge to unconstrained distributions generated (using the same final halo mass and redshift) for appropriate choices of the constraint parameters. In particular, in the limit M₁ → 0 (S₁ → ∞), our Brownian bridge excursions precisely converge to unconstrained excursions.¹² We quantify how growth histories converge to the unconstrained case by comparing main branch progenitor halo mass distributions as a function of M₁, at fixed M₀, |$z$|₀, and |$z$|₁, at several redshifts |$z$|⁠. The left-hand panel of Fig. B1 shows this result at |$z$| = 8 for merger trees constructed using a final halo mass of M₀ = 10¹⁴M_⊙ at |$z$|₀ = 0, with |$z$|₁ = 8, and demonstrates the desired behaviour. In particular, two-sample Kolmogorov–Smirnov (KS) tests with the unconstrained distribution yield KS statistics of 0.77, 0.6, 0.25, and 0.05, for M₁ = 10¹², 5 × 10¹¹, 10¹¹, and 10¹⁰, respectively. The M₁ = 10¹⁰M_⊙ case is consistent with the unconstrained distribution; in all other cases, the KS test rejects the null hypothesis of distinguishable distributions with p ≪ 0.05. We repeat these tests for various redshifts and constraint parameters, finding that constrained merger trees with M₁ ≲ 10⁻⁴M₀ yield main branch progenitor mass distributions consistent with corresponding unconstrained cases.

Figure B1.

Left: The distribution of main branch progenitor masses at |$z$| = 8 for a final halo of mass M₀ = 10¹⁴ M_⊙ at |$z$|₀ = 0. Results are shown for 1000 realizations of constrained merger trees with M₁ = 10¹² M_⊙  (blue), M₁ = 5 × 10¹¹M_⊙ (pink), M₁ = 10¹¹ M_⊙ (purple), and M₁ = 10¹⁰M_⊙ (brown), all with |$z$|₁ = 8, and for unconstrained merger trees with the same final halo properties (grey). Right: Comparison of growth histories for constrained branches from constrained merger tree realizations with |$z$|₁ = 8 and M₁ = 10¹¹M_⊙ (purple) versus branches of unconstrained trees that satisfy this constraint (grey).

We also test whether the subset of unconstrained merger tree realizations that satisfy a given constraint have growth histories that are consistent with those generated directly by our constrained solutions. The right-hand panel of Fig. B1 demonstrates that the growth histories of specific branches of unconstrained trees that satisfy a |$z$|₁ = 8, M₁ = 10¹¹M _⊙ constraint agree well with the distribution of constrained branches from our constrained merger trees. We find similar agreement for all constraint parameters studied in this work; thus, our constrained solutions accurately sample the underlying unconstrained distribution as desired.

B2 Convergence as a function of numerical resolution

To test for convergence as a function of the merger tree mass resolution, M_res, we repeat the tests from the previous subsection using three sets of merger tree realizations with the same constraint parameters (M₁ = 10¹³M_⊙ and |$z$|₁ = 8) and final halo parameters (M₀ = 10¹⁴ M _⊙ and |$z$|₀ = 0), with M_res = 10¹¹, 10¹⁰, and 10⁹M_⊙, respectively, at several redshifts. We find that the M_res = 10¹¹M_⊙ progenitor mass distribution is inconsistent with our fiducial M_res = 10¹⁰M_⊙ results, with p < 0.05 from KS two-sample tests at all redshifts 1 < |$z$| ≤ 8 (for |$z$| ≲ 1, the distributions agree by construction because the merger trees all reach the same final halo mass). Meanwhile, the M_res = 10⁹M_⊙ progenitor mass distributions are consistent with the M_res = 10¹⁰ M_⊙ case at all redshifts we test, with KS statistics ≤0.1, implying that our fiducial results are well converged for these constraint parameters. We repeat these tests for all constraints used in this paper, finding good convergence for M_res ≤ 10⁻⁴M₀ in all cases. We note that our constraints also introduce another convergence criterion – namely, that M_res ≤ M₁ so that the constrained-branch merger is always resolved. In practice, we find that M_res ≤ 10⁻¹M₁ is sufficient to achieve converged constrained-branch progenitor mass distributions. We have also tested the effects of merger tree construction resolution parameters in galacticus that control the resolution of the first-crossing rate distribution as a function of variance and time. None of our results are sensitive to variations of our fiducial choices for these parameters, which are documented in our publicly available code.

APPENDIX C: DESCENDANT HALO MASS DISTRIBUTIONS

In our constrained solutions, we can view the final halo mass M₀ and redshift |$z$|₀ as free parameters, given a constraint that the merger tree exceeds a mass of M₁ at some redshift |$z$|₁. Thus, for the high-redshift constraints considered in this work, it is natural to consider what distribution of halo masses, p(M₀, |$z$|₀), these haloes descend to at redshift |$z$|₀, assuming they self-consistently evolve to produce the halo mass function at |$z$| = |$z$|₀ calculated from standard unconstrained ePS solutions. In the context of our excursion set solutions, the probability that a trajectory which passes through (S₁, δ₁) will descend to a halo of mass S₀ at some fixed δ₀ is

$$ \begin{eqnarray*} p(S_0|\delta _0,\delta _1,S_1) \propto p(\delta _1|S_0,S_1,\delta _0)p(S_0,\delta _0) \propto p(\delta _1|S_0,S_1,\delta _0)\frac{\mathrm{d}n(S_0,\delta _0)}{\mathrm{d}S_0} \propto p(\delta _1|S_0,S_1,\delta _0)\frac{\mathrm{d}n(M_0[S_0],z_0(\delta _0))}{\mathrm{d}M_0}\left|\frac{\mathrm{d}M_0[S_0]}{\mathrm{d}S_0} \right|, \end{eqnarray*} $$

(C1)

where p(δ₁|S₀, S₁, δ₀)dδ₁ is the probability that a trajectory starting from S₀ (at some fixed δ₀) will pass through a small region (S₁, δ₁ ± dδ₁), dn(M₀[S₀], |$z$|₀(δ₀))/dM₀ is the halo mass function evaluated at the redshift |$z$|₀ corresponding to δ₀, and the Jacobian |dM₀[S₀]/dS₀| follows from equation (1). Lacey & Cole (1993) provide an expression for p(δ₁|S₀, S₁, δ₀)dδ₁, assuming that any branch of the merger tree can satisfy the constraint (specifically, these are the ‘constrained’ branches in our implementation). In our notation, their equation (2.16) becomes

$$ \begin{eqnarray*} p_{\mathrm{constrained\ branch}}(\delta _1|S_0,S_1,\delta _0) \propto \frac{1}{\sqrt{2\pi }}\left[\frac{S_1}{S_0(S_1-S_0)}\right]^{3/2}\frac{\delta _0(\delta _1-\delta _0)}{\delta _1}\times \exp \left[-\frac{(\delta _0 S_1-\delta _1 S_0)^2}{2S_1S_0(S_1-S_0)}\right], \end{eqnarray*} $$

(C2)

for S₁ > S₀ and δ₁ > δ₀, and zero otherwise. We also calculate the probability that a single branch of the merger tree (which we assume to be the main branch) satisfies the constraint and descends to (S₀, δ₀). Since any single (unconstrained) trajectory is an uncorrelated random walk with Gaussian transition probabilities, we have

$$ \begin{eqnarray*} p_{\mathrm{main\ branch}}(\delta _1|S_0,S_1,\delta _0) \propto \frac{1}{\sqrt{2\pi }} \frac{1}{\sqrt{S_1-S_0}} \times \exp \left[-\frac{(\delta _1-\delta _0)^2}{2(S_1-S_0)}\right], \end{eqnarray*} $$

(C3)

which again holds for S₁ > S₀ and δ₁ > δ₀ and is zero otherwise. Note that, this expression can also be interpreted in terms of the conditional progenitor mass function for the final halo.

We calculate the constrained and main-branch forms of the descendant distribution by assuming |$z$|₀ = 0 and plugging equations (C2) and (C3) into equation (C1), respectively, using a |$z$|₀ = 0 halo mass function computed from our unconstrained solutions. We then change variables from (S, δ) to (M, |$z$|⁠) by tabulating numerical solutions for the mass–variance relation and halo mass function. Fig. C1 shows the resulting descendant halo mass distributions for several constraints considered throughout this work. These distributions are suppressed at both high- and low-descendant masses. At high descendant masses, the |$z$|₀ = 0 halo mass function cuts off exponentially, while at low-descendant masses, the probability of a growth history exceeding the desired mass by redshift |$z$|₁ is also exponentially suppressed; for M₀ < M₁, the probability vanishes. Note that, the constrained-branch distribution is shifted towards larger descendant halo masses relative to the main-branch distribution. This follows because a specific branch that exceeds the high-redshift constraint continues growing for |$z$| < |$z$|₁, while any branch that satisfies the constraint can later merge onto the main branch. The results in Fig. C1 imply that a halo that exceeds mass M₁ = 10¹² M_⊙ (10¹¹M_⊙) at redshift |$z$|₁ = 8 most likely descends to a halo of mass M₀ = 6 × 10¹³M_⊙ (6 × 10¹²M_⊙) at |$z$|₀ = 0. Decreasing the constraint redshift at fixed progenitor mass also shifts the |$z$|₀ = 0 descendant distribution towards lower masses, as expected. For example, in the (M₁ = 10¹²M_⊙|$z$|₁ = 6) case, the most likely descendant halo mass is M₀ = 3 × 10¹³ M_⊙.

Figure C1.

Left: The distribution of |$z$|₀ = 0 descendant halo masses for excursions that exceed a mass of M₁ = 10¹² M_⊙ at |$z$|₁ = 8 (blue), M₁ = 10¹²M_⊙ at |$z$|₁ = 6 (green), and M₁ = 10¹¹M_⊙ at |$z$|₁ = 8 (purple). Solid (dashed) lines correspond to the case where any branch (a single branch) satisfies the constraint. Right: Same as the left-hand panel, for the seven high-redshift galaxy candidates from Labbe et al. (2022), using the progenitor halo masses estimated in Section 5.1.