Excursion Set First Crossing Statistics

Class providing the first-crossing distribution in the extended Press-Schechter excursion set framework—the probability that a random walk in the density field, starting from variance \(\sigma^2_0\) at time \(t\), first crosses the collapse barrier at a larger variance \(\sigma^2 > \sigma^2_0\). This distribution gives the halo mass function (via the unconditional first-crossing probability) and the conditional mass function of progenitors (via the transition rate between variances). Implementations depend on the choice of barrier shape and random walk statistics, with a flat barrier and Markovian walks recovering the analytic Press-Schechter result.

Default implementation: excursionSetFirstCrossingLinearBarrier

Methods

probability → double precision

Return the probability for a trajectory to make its first crossing of the barrier at the given variance and time.

  • double precision , intent(in ) :: variance, time

  • type (treeNode), intent(inout) :: node

rate → double precision

Return the rate of first crossing for excursion sets beginning at the given variance and time to transition to a first crossing at the given varianceProgenitor.

  • double precision , intent(in ) :: variance, varianceProgenitor, time

  • type (treeNode), intent(inout) :: node

rateNonCrossing → double precision

Return the rate of non-crossing for excursion sets beginning at the given variance and time.

  • double precision , intent(in ) :: variance, massMinimum, time

  • type (treeNode), intent(inout) :: node

coordinatedMPI → void

Sets the state of coordination under MPI. If set to true then the object can assume that any calculations it performs are being performed identically by all other MPI processes. This permits, for example, coordinated tabulation of results across MPI processes.

  • logical, intent(in ) :: state

excursionSetFirstCrossingFarahi

An excursion set first crossing statistics class using the algorithm of Benson et al. (2013). For trajectories originating from a point \((S_1,\delta_1)\), the distribution of first crossings of a barrier \(B(S)\), \(f(S)\), is obtained by finding the solution to the integral equation:

(5)\[ 1 = \int_0^S f(\tilde{S})\mathrm{d}\tilde{S} + \int_{-\infty}^{B(S)} P(\delta,S) \mathrm{d} \delta,\]

where \(P(\delta,S) \mathrm{d} \delta\) is the probability for a trajectory to lie between \(\delta\) and \(\delta + \mathrm{d} \delta\) at variance \(S\), having originated from the point \((S_1,\delta_1)\) having not crossed the barrier at any smaller \(\tilde{S} < S\). In the absence of a barrier, \(P(\delta,S)\) would be equal to \(P_0(\delta,S)\). However, since some trajectories will have crossed the barrier at \(\tilde{S} < S\) we must subtract off their contribution to \(P_0(\delta,S)\). Writing the distribution of \(\delta\) at \(S\) for trajectories originating at some \((\tilde{\delta},\tilde{S})\) as \(P^\prime(\delta|S,\tilde{\delta},\tilde{S})\), we can therefore write

\[P(\delta,S) = P_0(\delta,S) - \int_{0}^{S} f(\tilde{S}) P^\prime(\delta|S,\tilde{\delta},\tilde{S}) \mathrm{d}\tilde{S},\]

For an interesting class of cases, both \(P_0(\delta,S)\) and \(P^\prime(\delta|S,\tilde{\delta},\tilde{S})\) are normal distributions, and we can write

\[P_0(\delta,S) = \frac{1}{\sqrt{2 \pi S}} \exp\left\{-{\Delta \delta^2 [\delta,\delta_1,S,S_1] \over 2 \Delta S [S,S_1]}\right\},\]

and

\[P^\prime(\delta|S,\tilde{\delta},\tilde{S}) = \frac{1}{\sqrt{2 \pi \Delta S[S,\tilde{S}]}} \exp\left\{-{\Delta \delta^2 [\delta,\tilde{\delta},S,\tilde{S}] \over 2 \Delta S [S,\tilde{S}]}\right\},\]

where we refer to \(\Delta \delta[\delta,\tilde{\delta},S,\tilde{S}]\) as the “effective offset”, and to \(\Delta S[S,\tilde{S}] = \mathrm{Var}(S) - \mathrm{Cov}(S,\tilde{S})\) as the “residual variance”. Note that \(\mathrm{Cov}(S,S_1) = 0\) (since all trajectories pass through \(\delta_1\) at \(S_1\)), and so \(\Delta S[S,S_1] = \mathrm{Var}(S)\).

For a standard Weiner process (such as applies to the standard case considered in excursion set theory, namely uncorrelated and unconstrained steps), we have trivially that

\[\Delta \delta[\delta,\tilde{\delta},S,\tilde{S}] = \delta - \tilde{\delta},\]

and

\[\Delta S[S,\tilde{S}] = S - \tilde{S},\]

since the Weiner process is invariant under translations of the starting point.

Using the above results, we can rewrite eqn. ((5)):

\[1 = \int_0^S f(\tilde{S})\mathrm{d}\tilde{S} + \int_{-\infty}^{B(S)} \left[ P_0(\delta,S) - \int_{0}^{S} f(\tilde{S}) P^\prime(\delta|S,\tilde{\delta},\tilde{S}) \mathrm{d} \delta \right] ,\]

in general and, for the case of a Gaussian distribution:

\[1 = \int_0^S f(\tilde{S})\mathrm{d}\tilde{S} + \int_{-\infty}^{B(S)} \left[ \frac{1}{\sqrt{2 \pi \Delta S[S,S_1]}} \exp\left(-\frac{\Delta \delta^2[\delta,\delta_1,S,S_1]}{2 \Delta S[S,S_1]}\right) - \int_{0}^{S} f(\tilde{S}) \frac{1}{\sqrt{2 \pi \Delta S [S,\tilde{S}]}} \exp\left(-\frac{\Delta \delta^2[\delta,B(\tilde{S}),S,\tilde{S}]}{2 \Delta S [S,\tilde{S}]}\right)\mathrm{d}\tilde{S} \right] \mathrm{d} \delta .\]

The integral over \(\mathrm{d}\delta\) can be carried out analytically to give:

(6)\[ 1 = \int_0^S f(\tilde{S})\mathrm{d}\tilde{S}+ \hbox{erf}\left\{\frac{\Delta \delta [B(S),\delta_1,S,S_1]}{\sqrt{2\Delta S[S,S_1]}}\right\} - \int_{0}^{S} f(\tilde{S}) \hbox{erf}\left\{\frac{\Delta \delta [B(S),B(\tilde{S}),S,\tilde{S}]}{\sqrt{2 \Delta S [S,\tilde{S}]}}\right\} \mathrm{d}S^{\prime\prime}.\]

We now discretize eqn. ((6)). Specifically, we divide the \(S\) space into \(N\) intervals defined by the points:

\[\begin{split}S_i = \left\{ \begin{array}{ll} 0 & \hbox{if } i=0 \\ \sum_{j=0}^{i-1} \Delta S_j & \hbox{if } i > 1. \end{array} \right.\end{split}\]

Note that \(f(0)=0\) by definition, so \(f(S_0)=0\) always. We choose \(\Delta S_i = S_\mathrm{max}/N\) (i.e. uniform spacing in \(S\)) when computing first crossing distributions, and \(\Delta S_i \propto S_i\) (i.e. uniform spacing in \(\log(S)\)) when computing first crossing rates.

Discretizing the integrals in eqn. ((6)) gives:

(7)\[\int_0^{S_i} f(\tilde{S})\d \tilde{S} = \sum_{j=0}^{i-1} \frac{f(S_j) + f(S_{j+1})}{2} \Delta S_j\]

and:

(8)\[\int_{0}^{S_i} f(\tilde{S}) \hbox{erf}\left\{\frac{\Delta \delta [B(S),B(\tilde{S}),S,\tilde{S}]}{\sqrt{2 \Delta S[S,\tilde{S}]}}\right\} \d \tilde{S} = \sum_{j=0}^{i-1} \frac{1}{2} \left(f(S_j) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_j), S_i, S_j]}{\sqrt{2 \Delta S[S_i,S_j]}}\right\} + f(S_{j+1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_{j+1}),S_i,S_{j+1}]}{\sqrt{2 \Delta S[S_i,S_{j+1}]}}\right\} \right) \Delta S_j.\]

We can now rewrite eqn. ((6)) in discretized form:

(9)\[1 = \sum_{j=0}^{i-1} \frac{f(S_j) + f(S_{j+1})}{2} \Delta S_j + \hbox{erf}\left\{\frac{\Delta \delta [B(S_i),\delta_1,S_i,S_1]}{\sqrt{2 \Delta S[S_i,S_1]}}\right\} - \frac{1}{2} \sum_{j=0}^{i-1} \left( f(S_j) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_j),S_i,S_j]}{\sqrt{2 \Delta S[S_i,S_j]}}\right\} + f(S_{j+1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_{j+1}),S_i,S_{j+1}]}{\sqrt{2 \Delta S[S_i,S_{j+1}]}}\right\} \right) \Delta S_j.\]

Solving eqn. ((9)) for \(f(S_i)\):

(10)\[\begin{split}\left( \frac{1}{2} - \frac{1}{2} \hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_i), S_i, S_i]}{\sqrt{2 \Delta S[S_i,S_i]}}\right\} \right) \Delta S_{i-1} f(S_i) & = 1 - \sum_{j=0}^{i-2} \frac{f(S_j) + f(S_{j+1})}{2} \Delta S_j - \frac{f(S_{i-1})}{2} \Delta S_{i-1} - \hbox{erf}\left\{\frac{\Delta \delta [B(S_i),\delta_1,S_i,S_1]}{\sqrt{2 \Delta S[S_i,S_1]}}\right\} \nonumber \\ & + \frac{1}{2} \sum_{j=0}^{i-2} \left( f(S_j) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_j),S_i,S_j]}{\sqrt{2 \Delta S [S_i,S_j]}}\right\} + f(S_{j+1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_{j+1}),S_i,S_{j+1}]}{\sqrt{2 \Delta S[S_i,S_{j+1}]}}\right\} \right)\Delta S_j \nonumber \\ & + \frac{1}{2} f(S_{i-1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_{i-1}),S_i,S_{i-1}]}{\sqrt{2 \Delta S [S_i,S_{i-1}]}}\right\} \Delta S_{i-1}.\end{split}\]

For all barriers that we consider:

\[\hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_i),S_i,S_i]}{\sqrt{2 \Delta S[S_i,S_i]}}\right\} = 0.\]

We can then simplify eqn. ((10)):

(11)\[\begin{split}f(S_i) & = {2 \over \Delta S_{i-1}}\left[1 - \sum_{j=0}^{i-2} \frac{f(S_j) + f(S_{j+1})}{2} \Delta S_j - \frac{f(S_{i-1})}{2} \Delta S_{i-1} - \hbox{erf}\left\{\frac{\Delta \delta [B(S_i),\delta_1,S_i,S_1]}{\sqrt{2 \Delta S [S_i,S_1] }}\right\} \right. \nonumber \\ & + \frac{1}{2} \sum_{j=0}^{i-2} \left( f(S_j) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_j),S_i,S_j]}{\sqrt{2 \Delta S [S_i,S_j]}}\right\} + f(S_{j+1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_{j+1}),S_i,S_{j+1}]}{\sqrt{2 \Delta S [S_i,S_{j+1}]}}\right\} \right)\Delta S_j \nonumber \\ & \left. + \frac{1}{2} f(S_{i-1}) \hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_{i-1}),S_i,S_{i-1}]}{\sqrt{2 \Delta S [S_i,S_{i-1}]}}\right\} \Delta S_{i-1}\right].\end{split}\]

Consolidating terms in the summations:

(12)\[ f(S_i) = {2 \over \Delta S_{i-1}}\left[1 - \hbox{erf}\left\{\frac{\Delta \delta [B(S_i),\delta_1,S_i,S_1]}{\sqrt{2\Delta S[S_i,S_1]}}\right\} - \sum_{j=0}^{i-1} \left( 1-\hbox{erf}\left\{\frac{\Delta \delta [B(S_i) , B(S_j),S_i,S_j]}{\sqrt{2 \Delta S [S_i,S_j]}}\right\} \right) f(S_j) {\Delta S_{j-1} + \Delta S_j \over 2} \right].\]

In the case of constant \(\Delta S_j(=\Delta S)\) this can be simplified further:

(13)\[ f(S_i) = {2 \over \Delta S}\left[1 - \hbox{erf}\left\{\frac{\Delta \delta [B(S_i),\delta_1,S_i,S_1]}{\sqrt{2\Delta S [S_i,S_1]}}\right\}\right] - 2 \sum_{j=0}^{i-1} \left(1- \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(S_j),S_i,S_j]}{\sqrt{2 \Delta S[S_i,S_j]}}\right\} \right) f(S_j).\]

In either case (i.e. eqns. (12) and (13)) solution proceeds recursively: \(f(S_0)=0\) by definition, \(f(S_1)\) depends only on the known barrier and \(f(S_0)\), \(f(S_i)\) depends only on the known barrier and \(f(S_{<i})\).

The first crossing rate is computed using the same method but with an effective barrier which is offset by the position of the progenitor in the \((\delta,S)\) plane, plus a small shift in time. The non-crossing rate, \(g(S_\mathrm{max})\), defined as the rate at which trajectories reach the maximum possible variance, \(S_\mathrm{max}\), without ever crossing the barrier—is computed directly by integrating over the first crossing rate distribution, i.e. \(g(S_\mathrm{max}) = 1 -\int_0^{S_\mathrm{max}} f(\tilde{S}) \mathrm{d}\tilde{S}\). Note that since the numerical integration occurs only up to a finite maximum \(S_\mathrm{max}\), a non-zero non-crossing rate will be computed for CDM-like barriers even though in reality they should have zero non-crossing rate. As such, use of this method for such barriers is not recommended.

Methods

  • tabulate — Tabulate the virial density contrast as a function of mass and time.

  • restoreTable — Restore a tabulated solution from file.

  • storeTable — Store a tabulated solution to file.

Parameters

  • [fileName] (string; default none) — The name of the file to/from which tabulations of barrier first crossing probabilities should be written/read. If set to “none” tables will not be stored.

  • [fractionalTimeStep] (real; default 0.01d0) — The fractional time step used when computing barrier crossing rates (i.e. the step used in finite difference calculations).

  • [varianceNumberPerUnitProbability] (integer; default 1000) — The number of points to tabulate per unit variance for first crossing probabilities.

  • [varianceNumberPerUnit] (integer; default 40) — The number of tabulation points per unit of \(\sigma^2\) used when building the rate look-up table for the Farahi excursion-set first-crossing distribution; higher values improve interpolation accuracy at the cost of memory and initialization time.

  • [varianceNumberPerDecade] (integer; default 400) — The number of points to tabulate per decade of progenitor variance for first crossing rates.

  • [varianceNumberPerDecadeNonCrossing] (integer; default 40) — The number of points to tabulate per decade of progenitor variance for non-crossing rates.

  • [timeNumberPerDecade] (integer; default 10) — The number of tabulation points per decade of cosmic time used when building the first-crossing rate look-up table as a function of time; higher values improve temporal interpolation accuracy for rapidly evolving cosmologies.

  • [varianceIsUnlimited] (boolean; default .false.) — If true, the variance is assumed to have no upper limit (e.g. as in the case of CDM). This allows the tabulated solutions to be extended arbitrarily. Otherwise, tables are extended to encompass just the range of variance requested.

excursionSetFirstCrossingFarahiMidpoint

An excursion set first crossing statistics class using the algorithm of Benson et al. (2013), but using a midpoint method to perform the integrations (Du et al., 2017).

Specifically, in the method used by Benson et al. (2013) the integral equation for \(f(S)\) (see equation (5)) is solved using the trapezoidal rule (see excursionSetFirstCrossingFarahi for complete details):

\[\int_0^{S_j} f(\tilde{S})K(S_j,\tilde{S}) \mathrm{d}\tilde{S} = \sum_{j=0}^{i-1} \frac{f(S_j)K(S_i,S_j)+f(S_{j+1})K(S_i,S_{j+1})}{2} \Delta S_j,\]

where

\[K(S_i,\tilde{S}) = \hbox{erf}\left\{\frac{\Delta \delta [B(S_i), B(\tilde{S}), S_i, \tilde{S}]}{\sqrt{2 \Delta S[S_i,\tilde{S}]}}\right\},\]

is the kernel of the integral equation for \(f(S)\) which is being solved.

The method used in this class increases the precision of the solver by replacing the trapezoidal integration rule in the above with a mid-point integration rule:

\[\int_0^{S_j} f(\tilde{S})K(S_j,\tilde{S}) \mathrm{d}\tilde{S} = \sum_{j=0}^{i-1} f(S_{j+1/2})K(S_i,S_{j+1/2}) \Delta S_j,\]

such that the first crossing distribution at \(S_{i-1/2}\) is given by

\[f(S_{i-1/2}) = \frac{1}{K(S_i,S_{i-1/2})\Delta S_{i-1/2}} \left( \hbox{erfc}\left\{ \frac{\Delta \delta [B(S_i),B(S_1),S_i,S_1]}{\sqrt{2 \Delta S[S_i,S_1]}} \right\} - \sum_{j=0}^{i-2} f(S_{j+1/2} K(S_i,S_{j+1/2}) \Delta S_j) \right).\]

Methods

  • tabulate — Tabulate the virial density contrast as a function of mass and time.

  • restoreTable — Restore a tabulated solution from file.

  • storeTable — Store a tabulated solution to file.

Parameters

  • [fileName] (string; default none) — The name of the file to/from which tabulations of barrier first crossing probabilities should be written/read. If set to “none” tables will not be stored. (inherited from excursionSetFirstCrossingFarahi)

  • [fractionalTimeStep] (real; default 0.01d0) — The fractional time step used when computing barrier crossing rates (i.e. the step used in finite difference calculations). (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerUnitProbability] (integer; default 1000) — The number of points to tabulate per unit variance for first crossing probabilities. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerUnit] (integer; default 40) — The number of tabulation points per unit of \(\sigma^2\) used when building the rate look-up table for the Farahi excursion-set first-crossing distribution; higher values improve interpolation accuracy at the cost of memory and initialization time. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerDecade] (integer; default 400) — The number of points to tabulate per decade of progenitor variance for first crossing rates. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerDecadeNonCrossing] (integer; default 40) — The number of points to tabulate per decade of progenitor variance for non-crossing rates. (inherited from excursionSetFirstCrossingFarahi)

  • [timeNumberPerDecade] (integer; default 10) — The number of tabulation points per decade of cosmic time used when building the first-crossing rate look-up table as a function of time; higher values improve temporal interpolation accuracy for rapidly evolving cosmologies. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceIsUnlimited] (boolean; default .false.) — If true, the variance is assumed to have no upper limit (e.g. as in the case of CDM). This allows the tabulated solutions to be extended arbitrarily. Otherwise, tables are extended to encompass just the range of variance requested. (inherited from excursionSetFirstCrossingFarahi)

excursionSetFirstCrossingFarahiMidpointBrownianBridge

An excursion set first crossing statistics class using the algorithm of Benson et al. (2013), but using a midpoint method to perform the integrations (Du et al., 2017), and with a Brownian bridge constraint.

Specifically, the trajectories are constrained to pass through a point \((S_2,\delta_2)\) (specified by the parameter [varianceConstrained] and [criticalOverdensityConstrained], or equivalently by the parameters [massConstrained] and [timeConstrained]—note that \((S_2,\delta_2)\) here follow the convention in excursion set literature that \(S_2\) is the variance evaluated at the present day, while \(\delta_2\) the the critical overdensity for collapse divided by the linear growth factor), and, of course, always pass through the initial point \((S_1,\delta_1)\) corresponding to the current halo.

For a Brownian bridge the distribution of \(\delta\) at some \(S\) (where \(S_1 \le S \le S_2\)), \(P_0(\delta|S)\), is given by a normal distribution with mean

\[\mu(S) = \delta_1 + \frac{\delta_2-\delta_1}{S_2-S_1}(S - S_1)\]

and variance

\[\mathrm{Var}(S) = \frac{(S_2-S)(S-S_1)}{S_2-S_1},\]

and the covariance between two points \(S_\mathrm{a}\) and \(S_\mathrm{b} > S_\mathrm{a}\) is given by,

\[\mathrm{Cov}(S_\mathrm{a},S_\mathrm{b}) = \frac{(S_2-S_\mathrm{b})(S_\mathrm{a}-S_1)}{S_2-S_1}.\]

Therefore, the same approach to solving for the first crossing distribution as was utilized by Benson et al. (2013) and improved by (Du et al., 2017) can be used (see excursionSetFirstCrossingFarahi and excursionSetFirstCrossingFarahiMidpoint for details), with just the appropriate change in the effective offset, \(\Delta \delta\), and residual variance, \(\Delta S\).

Considering two points \((S,\delta)\) and \((\tilde{S},\tilde{\delta})\) the effective offset is just the difference in their offsets relative to their local means:

\[\Delta \delta = [ \delta - \mu(S) ] - [ \tilde{\delta} - \mu(\tilde{S}) ] = \delta - \tilde{\delta} - \frac{S-\tilde{S}}{S_2-S_1}(\delta_2-\delta_1),\]

while the residual variance is, as always, just the variance at \((S,\delta)\) minus the covariance between the two points:

\[\Delta S = \mathrm{Var}(S) - \mathrm{Cov}(B({\tilde{S}}),\delta) = \frac{(S_2-S)(S-S_1)}{S_2-S_1} - \frac{(S_2-S)(\tilde{S}-S_1)}{S_2-S_1}.\]

which simplifies to

\[\Delta S = \frac{(S_2-S)(S-\tilde{S})}{S_2-S_1}.\]

Note that, in solving for the first crossing distribution we must also evaluate terms of the form

\[\int_{-\infty}^{\delta} P_{0}(\delta^\prime,S) \mathrm{d}\delta^\prime = \mathrm{erf}\left( \frac{\delta^\prime - \mu(S)}{\sqrt{2 \mathrm{Var}(S)}}\right).\]

In these cases we still use the residual variance since \(\Delta S \rightarrow \mathrm{Var}(S)\) as \(\tilde{S} \rightarrow S_1\).

When computing the distribution, \(p(\delta,s)\), of trajectories at variance \(S\), given that the first crossed the barrier, \(B(\tilde{S})\), at some smaller variance, \(\tilde{S}\), (equation A2 of Benson et al. 2013) we must condition both the residual variance and drift term on the intermediate point, \(\tilde{S},B(\tilde{S})\). Fortunately, given any Brownian random walk (including Brownian bridges) for which we know two points, the distribution of trajectories between those points is simply another Brownian bridge. Therefore, we can write:

\[\Delta \delta = \delta - \tilde{\delta} - \frac{S-\tilde{S}}{S_2-\tilde{S}}(\delta_2-\tilde{\delta}),\]

and:

\[\Delta S = \frac{(S_2-S)(S-\tilde{S})}{S_2-\tilde{S}}.\]

This class provides functions implementing these modified effective offset and residual variance.

Parameters

  • [criticalOverdensityConstrained] (real) — The linear theory critical overdensity \(\delta_\mathrm{c}\) (extrapolated to the present epoch) that defines the constrained end-point of the Brownian bridge in excursion-set space; used together with varianceConstrained to pin the random walk to a specific progenitor halo.

  • [varianceConstrained] (real) — The mass variance \(\sigma^2(M)\) corresponding to the constrained end-point mass of the Brownian bridge; together with criticalOverdensityConstrained it specifies the progenitor mass to which the excursion-set random walk is conditioned.

  • [redshiftConstrained] (real) — The redshift of the progenitor epoch that defines the constrained end-point of the Brownian bridge; converted internally to a cosmic time and then to a linear overdensity threshold via the critical overdensity at that epoch.

  • [massConstrained] (real) — The halo mass (\(\mathrm{M}_\odot\)) of the constrained progenitor at the end of the Brownian bridge; converted internally to a mass variance \(\sigma^2(M)\) that pins the excursion-set random walk to a specific progenitor scale.

  • [fileName] (string; default none) — The name of the file to/from which tabulations of barrier first crossing probabilities should be written/read. If set to “none” tables will not be stored. (inherited from excursionSetFirstCrossingFarahi)

  • [fractionalTimeStep] (real; default 0.01d0) — The fractional time step used when computing barrier crossing rates (i.e. the step used in finite difference calculations). (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerUnitProbability] (integer; default 1000) — The number of points to tabulate per unit variance for first crossing probabilities. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerUnit] (integer; default 40) — The number of tabulation points per unit of \(\sigma^2\) used when building the rate look-up table for the Farahi excursion-set first-crossing distribution; higher values improve interpolation accuracy at the cost of memory and initialization time. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerDecade] (integer; default 400) — The number of points to tabulate per decade of progenitor variance for first crossing rates. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceNumberPerDecadeNonCrossing] (integer; default 40) — The number of points to tabulate per decade of progenitor variance for non-crossing rates. (inherited from excursionSetFirstCrossingFarahi)

  • [timeNumberPerDecade] (integer; default 10) — The number of tabulation points per decade of cosmic time used when building the first-crossing rate look-up table as a function of time; higher values improve temporal interpolation accuracy for rapidly evolving cosmologies. (inherited from excursionSetFirstCrossingFarahi)

  • [varianceIsUnlimited] (boolean; default .false.) — If true, the variance is assumed to have no upper limit (e.g. as in the case of CDM). This allows the tabulated solutions to be extended arbitrarily. Otherwise, tables are extended to encompass just the range of variance requested. (inherited from excursionSetFirstCrossingFarahi)

excursionSetFirstCrossingLinearBarrier

An excursion set first crossing statistics class for linear barriers. Specifically, the first crossing distribution is

\[f(S,t) = B(0,t) \exp(- B(S,t)^2/2S)/S/\sqrt{2 pi S},\]

where \(B(S,t)\) is the (assumed-to-be-linear-in-\(S\)) barrier at time \(t\) and variance \(S\). The first crossing rate is computed using a finite difference approximation between two closely-spaced times. The non-crossing rate is zero.

(Default implementation)

Methods

  • tabulate — Tabulate the virial density contrast as a function of mass and time.

  • restoreTable — Restore a tabulated solution from file.

  • storeTable — Store a tabulated solution to file.

Parameters

  • [fractionalTimeStep] (real; default 0.01d0) — The fractional time step used when computing barrier crossing rates (i.e. the step used in finite difference calculations).

excursionSetFirstCrossingLinearBarrierBrownianBridge

An excursion set first crossing statistics class for linear barriers and where the trajectories are constrained to follow a Brownian bridge.

If we consider the Brownian bridge to originate from \((0,0)\) (i.e. we apply the usual shift of coordinates to move our starting point to the origin), and to end at \((\delta_2,S_2)\) then we can transform this Brownian bridge into the standard bridge with non-zero drift through the transformations:

\[\begin{split}\tau & = \frac{S}{S_2}, \\ b & = \frac{\delta_2}{\sqrt{S_2}}.\end{split}\]

To find the first crossing time distribution we then follow the general approach outlined by Kiwiakos (2014), but with an important difference that we will detail below.

The standard Brownian bridge (with no drift), \(Y_0\), can be written in terms of a standard Weiner process, \(W\), through a change of variables

\[Y_0(t) = (1-t) W\left(\frac{t}{1-t}\right).\]

The first crossing time distribution for our Brownian bridge can therefore be expressed as:

\[\tau_{Y}(B) = \mathrm{inf}\left\{ t : Y(t) = B(t)\right\} = \mathrm{inf}\left\{ \mu(t) + (1-t) W\left(\frac{t}{1-t}\right) \right\} = B(t) = \mathrm{inf}\left\{ t : W\left(\frac{t}{1-t}\right) = B(t) - \mu(t) \right\},\]

where \(B(t)\) is our barrier, and \(\mu(t)\) is the drift term in the Brownian bridge.

As can be seen from the above, for the case of a linear barrier, a Brownian bridge with non-zero drift effectively results in a new linear barrier equal to the original one minus the drift term, i.e.:

\[B^\prime(t) \rightarrow B(t) - \mu(t),\]

where \(B(S)\) is the barrier and \(\mu(S)\) is the Brownian bridge term. This means that the first-crossing time of the Brownian bridge is just the hitting time of this time changed Weiner process. That is, is

\[\tau_W(B) = \mathrm{inf}\left\{ s : W(s) = B\right\},\]

then

\[\frac{\tau_Y(B)}{1 - \tau_Y(B)} = \tau_W(B) \implies \tau_Y(B) = \frac{\tau_W(B)}{1+\tau_W(B)}.\]

Here is where the solution presented by Kiwiakos (2014) is slightly wrong. We must use the first crossing time solution for the standard Weiner process, but with a linear barrier (because, even if the actual barrier is constant, the effective barrier is linear due to the Brownian bridge drift term). Therefore (e.g. Zhang and Hui 2006):

\[f_{W}(\tau_{W}) = B(0) \exp\left( - \frac{B(\tau_{W})^2}{2\ tau_{W}} \right) / \sqrt{2 pi \tau_{W}^3}.\]

We then have that

\[f_{Y}(\tau_{Y}) \mathrm{d}\tau_{Y} = f_{W}(\tau_{W}) \mathrm{d}\tau_{W},\]

such that

\[f_{Y}(\tau_{Y}) = \frac{B(0)}{\sqrt{2 \pi \tau_{Y}^3 (1 - \tau_{Y})}} \exp\left( \frac{B^{\prime 2}(\tau_{Y})}{2 \tau_{Y} (1-\tau_{Y})} \right),\]

or, expressed in our usual variables

\[f(S) = \frac{S(0)}{\sqrt{2 \pi S^3 (1 - S/S_2)}} \exp\left( \frac{[B(S)-\mu(S)]^2}{2 S (1-S/S_2)} \right).\]

Parameters

  • [fractionalTimeStep] (real; default 0.01d0) — The fractional time step used when computing barrier crossing rates (i.e. the step used in finite difference calculations).

  • [criticalOverdensityConstrained] (real) — The linear theory critical overdensity \(\delta_\mathrm{c}\) (extrapolated to the present epoch) that defines the constrained end-point of the Brownian bridge in excursion-set space; used together with varianceConstrained to pin the random walk to a specific progenitor halo.

  • [varianceConstrained] (real) — The mass variance \(\sigma^2(M)\) corresponding to the constrained end-point mass of the Brownian bridge; together with criticalOverdensityConstrained it specifies the progenitor mass to which the excursion-set random walk is conditioned.

  • [redshiftConstrained] (real) — The redshift of the progenitor epoch that defines the constrained end-point of the Brownian bridge; converted internally to a cosmic time and then to a linear overdensity threshold via the critical overdensity at that epoch.

  • [massConstrained] (real) — The halo mass (\(\mathrm{M}_\odot\)) of the constrained progenitor at the end of the Brownian bridge; converted internally to a mass variance \(\sigma^2(M)\) that pins the excursion-set random walk to a specific progenitor scale.

excursionSetFirstCrossingZhangHui

An excursion set first crossing statistics class utilizing the algorithm of Zhang and Hui (2006). First crossing (and non-crossing) rates are not supported by this method.

Methods

excursionSetFirstCrossingZhangHuiHighOrder

An excursion set first crossing statistics class utilizing a higher order generalization of the algorithm of Zhang and Hui (2006). First crossing (and non-crossing) rates are not supported by this method.

In this method we discretize the first-crossing distribution function, and use a closed Newton-Cotes method to perform integrations. First, we mesh the \(S\) space using uniform spacing in \(S\):

\[S_i = i \times \Delta S , i = 0, 1, \dots, N, \Delta S = \frac{S}{N}.\]

Then we discretize the integral equation by Boole’s rule. The integral equation becomes a set of linear algebraic equations:

\[\begin{split}f(S_i) & = g_1(S_i) + \frac{ \Delta S}{90} \sum_{j=0:4}^{i-4} \left\{ 7 f(S_j) g_2(S_i, S_j) + 32 f(S_{j+1}) g_2(S_i, S_{j+1}) \right.\nonumber \\ & \qquad \left. {} + 12 f(S_{j+2}) g_2(S_i, S_{j+2}) + 32 f(S_{j+3}) g_2(S_i, S_{j+3}) + 7 f(S_{j+4}) g_2(S_i, S_{j+4}) \right\}.\end{split}\]

Since \(g_2(S, S^\prime)\) approaches infinity when \(S\) approaches \(S^\prime\), one needs to define \(g_2(S_i, S_i)\) carefully when \(j = i\). We can rewrite the equation:

\[\begin{split}f(S_i) & = g_1(S_i) + \frac{ 4 \Delta S}{90} \sum_{j=0:4}^{i-8} \left\{ 7 f(S_j) g_2(S_i, S_j) + 32 f(S_{j+1}) g_2(S_i, S_{j+1}) \right.\nonumber \\ & \qquad \left. {} + 12 f(S_{j+2}) g_2(S_i, S_{j+2}) + 32 f(S_{j+3}) g_2(S_i, S_{j+3}) + 7 f(S_{j+4}) g_2(S_i, S_{j+4}) \right\} \nonumber \\ & \qquad {} + \frac{ 4 \bar{g}_2(S_i) \Delta S}{90} \Big( 7 f(S_{i-4}) + 32 f(S_{i-3}) + 12 f(S_{i-2}) + 32 f(S_{i-1}) + 7 f(S_i) \Big).\end{split}\]

For \(\bar{g}_2(S_i)\) we have:

\[\bar{g}_2(S_i) = \frac{1}{\delta S} \int_{S - \delta S}^S g_2(S,S^\prime)\mathrm{d}S^\prime.\]

In the above, \(\delta S\) depends on the range of the previous integral part. Generally, \(\delta S\) is equal to \(4 \Delta S\). The above equation can be solved for \(f(S_i)\), giving:

\[\begin{split}f(S_i) & = \left( g_1(S_i) + \frac{ 4 \Delta S}{90} \sum_{j=0:4}^{i-8} \left\{ 7 f(S_j) g_2(S_i, S_j) + 32 f(S_{j+1}) g_2(S_i, S_{j+1}) \right.\right.\nonumber \\ & \qquad \left.\left. {} + 12 f(S_{j+2}) g_2(S_i, S_{j+2}) + 32 f(S_{j+3}) g_2(S_i, S_{j+3}) + 7 f(S_{j+4}) g_2(S_i, S_{j+4}) \right\} \right. \nonumber \\ & \qquad \left. {} + \frac{4 \bar{g}_2(S_i) \Delta S }{90} \Big( 7 f(S_{i-4}) + 32 f(S_{i-3}) + 12 f(S_{i-2}) + 32 f(S_{i-1}) \Big) \right) \nonumber \\ & \qquad {} \times \left( 1 - \frac{ 4 \bar{g}_2(S_i) \Delta S}{90} \right)^{-1}.\end{split}\]

Not all of the \(i\)’s are divisible by \(4\). So, for the first \(m^\mathrm{th}\) spaces, we need to calculate the integral part, separately, where \(m\) is the remainder of \(i\) by the modulus of \(4\). It is a good approximation to calculate the first part linearly. Consequently, the final formula for the general problem is:

\[\begin{split}f(S_i) & = \left( g_1(S_i) + \frac{ \Delta S}{2} \sum_{j=0}^{m-1} \Big( f(S_j) g_2(S_i, S_j) + f(S_{j+1}) g_2(S_i, S_{j+1}) \Big) \right.\nonumber \\ & \qquad \left. {} + \frac{ 4 \Delta S}{90} \sum_{\frac{j-m}{4}=0}^{i-8} \left\{ 7 f(S_j) g_2(S_i, S_j) + 32 f(S_{j+1}) g_2(S_i, S_{j+1}) \right.\right.\nonumber \\ & \qquad \left.\left. {} + 12 f(S_{j+2}) g_2(S_i, S_{j+2}) + 32 f(S_{j+3}) g_2(S_i, S_{j+3}) + 7 f(S_{j+4}) g_2(S_i, S_{j+4}) \right\} \right. \nonumber \\ & \qquad \left. {} + \frac{4 \bar{g}_2(S_i) \Delta S}{90} \Big( 7 f(S_{i-4}) + 32 f(S_{i-3}) + 12 f(S_{i-2}) + 32 f(S_{i-1}) \Big) \right) \nonumber \\ & \qquad {} \times \left( 1 - \frac{ 4 \bar{g}_2(S_i) \Delta S }{90} \right)^{-1}.\end{split}\]

Finally, we can solve the first-crossing distribution function, giving:

\[\begin{split}f(S_0) & = g_1(S_0) = 0 \\ f(S_1) & = g_1(S_1) \left( 1 - \frac{\bar{g}_2(S_1) \Delta S}{2} \right)^{-1} \\ f(S_2) & = \Big( g_1(S_2) + \frac{ \Delta S}{2} 2 \bar{g}_2(S_2) f(S_1) \Big) \left( 1 - \frac{ \bar{g}_2(S_2) \Delta S}{2} \right)^{-1} \\ f(S_3) & = \Big( g_1(S_2) + \frac{ \Delta S}{2} 2 \bar{g}_2(S_3) (f(S_1) + f(S_2)) \Big) \left( 1 - \frac{ \bar{g}_2(S_3) \Delta S}{2} \right)^{-1} \\ f(S_i) & = \left( g_1(S_i) + \frac{ \Delta S}{2} \sum_{j=0}^{m-1} \Big( f(S_j) g_2(S_i, S_j) + f(S_{j+1}) g_2(S_i, S_{j+1}) \Big) \right.\nonumber \\ & \qquad \left. {} + \frac{ 4 \Delta S}{90} \sum_{\frac{j-m}{4}=0}^{i-8} \left\{ 7 f(S_j) g_2(S_i, S_j) + 32 f(S_{j+1}) g_2(S_i, S_{j+1}) \right.\right.\nonumber \\ & \qquad \left.\left. {} + 12 f(S_{j+2}) g_2(S_i, S_{j+2}) + 32 f(S_{j+3}) g_2(S_i, S_{j+3}) + 7 f(S_{j+4}) g_2(S_i, S_{j+4}) \right\} \right. \nonumber \\ & \qquad \left. {} + \frac{4 \bar{g}_2(S_i) \Delta S}{90} \Big( 7 f(S_{i-4}) + 32 f(S_{i-3})+ 12 f(S_{i-2}) + 32 f(S_{i-1}) \Big) \right) \nonumber \\ & \qquad {} \times \left( 1 - \frac{ 4 \bar{g}_2(S_i) \Delta S}{90} \right)^{-1},\end{split}\]

where \(\bar{g}_2(S_1)\) is defined as:

\[\begin{split}\bar{g}_2(S_1) & = \frac{1}{\Delta S} \int_{0}^{S_1} g_2(S,S^\prime)\mathrm{d}S^\prime \hbox{ if } i = 1 \\ \bar{g}_2(S_2) & = \frac{1}{2 \Delta S} \int_{0}^{S_2} g_2(S,S^\prime)\mathrm{d}S^\prime \hbox{ if } i = 0 \\ \bar{g}_2(S_3) & = \frac{1}{3 \Delta S} \int_{0}^{S_3} g_2(S,S^\prime)\mathrm{d}S^\prime \hbox{ if } i = 0 \\ \bar{g}_2(S_i) & = \frac{1}{4 \Delta S} \int_{S_i - 4 \Delta S}^{S_i} g_2(S,S^\prime)\mathrm{d}S^\prime \hbox{ if } i > 3.\end{split}\]

The error term for this method of discretization is:

\[\epsilon = \frac{(\delta S)^7}{1935360} f^{(6)}(\xi),\]

where \(f^{(6)}(\xi)\) is the absolute maximum of the sixth derivative in the range of \([S_i , S_i + \delta S]\). For this problem, \(\delta S = 4 \Delta S\), so:

\[\epsilon \leq \frac{(\Delta S)^7}{118.125} f^{(6)}(\xi).\]

Since there are \(N/4\) intervals, the maximum deviation from the real value of the function is:

\[\begin{split}\epsilon \leq \sum_{i=1}{\frac{N}{4}}\frac{(\Delta S)^7}{118.125} f^{(6)}(\xi _i) \nonumber\\ \leq N4\frac{(\Delta S)^7}{472.5} f^{(6)}(\xi),\end{split}\]

where \(f^{(6)}(\xi)\) is the absolute maximum of the sixth derivative in the domain, \([0 , S]\).

Methods

  • initialize — Initialize the high order Zhang and Hui (2006) class.