.. _physics-excursionSetFirstCrossing: 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 :math:`\sigma^2_0` at time :math:`t`, first crosses the collapse barrier at a larger variance :math:`\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`` .. _physics-excursionSetFirstCrossingFarahi: ``excursionSetFirstCrossingFarahi`` ----------------------------------- An excursion set first crossing statistics class using the algorithm of :cite:t:`benson_dark_2012`. For trajectories originating from a point :math:`(S_1,\delta_1)`, the distribution of first crossings of a barrier :math:`B(S)`, :math:`f(S)`, is obtained by finding the solution to the integral equation: .. math:: :label: eq-OldExcursionMethod 1 = \int_0^S f(\tilde{S})\mathrm{d}\tilde{S} + \int_{-\infty}^{B(S)} P(\delta,S) \mathrm{d} \delta, where :math:`P(\delta,S) \mathrm{d} \delta` is the probability for a trajectory to lie between :math:`\delta` and :math:`\delta + \mathrm{d} \delta` at variance :math:`S`, having originated from the point :math:`(S_1,\delta_1)` having not crossed the barrier at any smaller :math:`\tilde{S} < S`. In the absence of a barrier, :math:`P(\delta,S)` would be equal to :math:`P_0(\delta,S)`. However, since some trajectories will have crossed the barrier at :math:`\tilde{S} < S` we must subtract off their contribution to :math:`P_0(\delta,S)`. Writing the distribution of :math:`\delta` at :math:`S` for trajectories originating at some :math:`(\tilde{\delta},\tilde{S})` as :math:`P^\prime(\delta|S,\tilde{\delta},\tilde{S})`, we can therefore write .. math:: 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 :math:`P_0(\delta,S)` and :math:`P^\prime(\delta|S,\tilde{\delta},\tilde{S})` are normal distributions, and we can write .. math:: 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 .. math:: 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 :math:`\Delta \delta[\delta,\tilde{\delta},S,\tilde{S}]` as the "effective offset", and to :math:`\Delta S[S,\tilde{S}] = \mathrm{Var}(S) - \mathrm{Cov}(S,\tilde{S})` as the "residual variance". Note that :math:`\mathrm{Cov}(S,S_1) = 0` (since all trajectories pass through :math:`\delta_1` at :math:`S_1`), and so :math:`\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 .. math:: \Delta \delta[\delta,\tilde{\delta},S,\tilde{S}] = \delta - \tilde{\delta}, and .. math:: \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. (:eq:`eq-OldExcursionMethod`): .. math:: 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: .. math:: 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 :math:`\mathrm{d}\delta` can be carried out analytically to give: .. math:: :label: eq-NewExcursionMethod 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. (:eq:`eq-NewExcursionMethod`). Specifically, we divide the :math:`S` space into :math:`N` intervals defined by the points: .. math:: 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. Note that :math:`f(0)=0` by definition, so :math:`f(S_0)=0` always. We choose :math:`\Delta S_i = S_\mathrm{max}/N` (i.e. uniform spacing in :math:`S`) when computing first crossing distributions, and :math:`\Delta S_i \propto S_i` (i.e. uniform spacing in :math:`\log(S)`) when computing first crossing rates. Discretizing the integrals in eqn. (:eq:`eq-NewExcursionMethod`) gives: .. math:: :label: eq-Des1 \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: .. math:: :label: eq-Des2 \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. (:eq:`eq-NewExcursionMethod`) in discretized form: .. math:: :label: eq-DesFinal1 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. (:eq:`eq-DesFinal1`) for :math:`f(S_i)`: .. math:: :label: eq-DesFinal11 \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}. For all barriers that we consider: .. math:: \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. (:eq:`eq-DesFinal11`): .. math:: :label: eq-DesFinal2 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]. Consolidating terms in the summations: .. math:: :label: eq-DesFinal2a 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 :math:`\Delta S_j(=\Delta S)` this can be simplified further: .. math:: :label: eq-DesFinal3 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. :eq:`eq-DesFinal2a` and :eq:`eq-DesFinal3`) solution proceeds recursively: :math:`f(S_0)=0` by definition, :math:`f(S_1)` depends only on the known barrier and :math:`f(S_0)`, :math:`f(S_i)` depends only on the known barrier and :math:`f(S_{`_ constraint. Specifically, the trajectories are constrained to pass through a point :math:`(S_2,\delta_2)` (specified by the parameter ``[varianceConstrained]`` and ``[criticalOverdensityConstrained]``, or equivalently by the parameters ``[massConstrained]`` and ``[timeConstrained]``---note that :math:`(S_2,\delta_2)` here follow the convention in excursion set literature that :math:`S_2` is the variance evaluated at the present day, while :math:`\delta_2` the the critical overdensity for collapse divided by the linear growth factor), and, of course, always pass through the initial point :math:`(S_1,\delta_1)` corresponding to the current halo. For a Brownian bridge the distribution of :math:`\delta` at some :math:`S` (where :math:`S_1 \le S \le S_2`), :math:`P_0(\delta|S)`, is given by a normal distribution with mean .. math:: \mu(S) = \delta_1 + \frac{\delta_2-\delta_1}{S_2-S_1}(S - S_1) and variance .. math:: \mathrm{Var}(S) = \frac{(S_2-S)(S-S_1)}{S_2-S_1}, and the covariance between two points :math:`S_\mathrm{a}` and :math:`S_\mathrm{b} > S_\mathrm{a}` is given by, .. math:: \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 :cite:t:`benson_dark_2012` and improved by :cite:p:`du_substructure_2017` can be used (see :galacticus-class:`excursionSetFirstCrossingFarahi` and :galacticus-class:`excursionSetFirstCrossingFarahiMidpoint` for details), with just the appropriate change in the effective offset, :math:`\Delta \delta`, and residual variance, :math:`\Delta S`. Considering two points :math:`(S,\delta)` and :math:`(\tilde{S},\tilde{\delta})` the effective offset is just the difference in their offsets relative to their local means: .. math:: \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 :math:`(S,\delta)` minus the covariance between the two points: .. math:: \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 .. math:: \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 .. math:: \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 :math:`\Delta S \rightarrow \mathrm{Var}(S)` as :math:`\tilde{S} \rightarrow S_1`. When computing the distribution, :math:`p(\delta,s)`, of trajectories at variance :math:`S`, given that the first crossed the barrier, :math:`B(\tilde{S})`, at some smaller variance, :math:`\tilde{S}`, (equation A2 of :cite:author:`benson_dark_2012` :cite:year:`benson_dark_2012`) we must condition *both* the residual variance and drift term on the intermediate point, :math:`\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: .. math:: \Delta \delta = \delta - \tilde{\delta} - \frac{S-\tilde{S}}{S_2-\tilde{S}}(\delta_2-\tilde{\delta}), and: .. math:: \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 :math:`\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 :math:`\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 (:math:`\mathrm{M}_\odot`) of the constrained progenitor at the end of the Brownian bridge; converted internally to a mass variance :math:`\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 :math:`\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 :term:`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``\ *)* .. _physics-excursionSetFirstCrossingLinearBarrier: ``excursionSetFirstCrossingLinearBarrier`` ------------------------------------------ An excursion set first crossing statistics class for linear barriers. Specifically, the first crossing distribution is .. math:: f(S,t) = B(0,t) \exp(- B(S,t)^2/2S)/S/\sqrt{2 pi S}, where :math:`B(S,t)` is the (assumed-to-be-linear-in-:math:`S`) barrier at time :math:`t` and variance :math:`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). .. _physics-excursionSetFirstCrossingLinearBarrierBrownianBridge: ``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 :math:`(0,0)` (i.e. we apply the usual shift of coordinates to move our starting point to the origin), and to end at :math:`(\delta_2,S_2)` then we can transform this Brownian bridge into the standard bridge with non-zero drift through the transformations: .. math:: \tau & = \frac{S}{S_2}, \\ b & = \frac{\delta_2}{\sqrt{S_2}}. To find the first crossing time distribution we then follow the general approach outlined by :cite:t:`kiwiakos_answer_2014`, but with an important difference that we will detail below. The standard Brownian bridge (with no drift), :math:`Y_0`, can be written in terms of a standard Weiner process, :math:`W`, through a change of variables .. math:: 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: .. math:: \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 :math:`B(t)` is our barrier, and :math:`\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.: .. math:: B^\prime(t) \rightarrow B(t) - \mu(t), where :math:`B(S)` is the barrier and :math:`\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 .. math:: \tau_W(B) = \mathrm{inf}\left\{ s : W(s) = B\right\}, then .. math:: \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 :cite:t:`kiwiakos_answer_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. :cite:author:`zhang_random_2006` :cite:year:`zhang_random_2006`): .. math:: 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 .. math:: f_{Y}(\tau_{Y}) \mathrm{d}\tau_{Y} = f_{W}(\tau_{W}) \mathrm{d}\tau_{W}, such that .. math:: 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 .. math:: 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 :math:`\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 :math:`\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 (:math:`\mathrm{M}_\odot`) of the constrained progenitor at the end of the Brownian bridge; converted internally to a mass variance :math:`\sigma^2(M)` that pins the excursion-set random walk to a specific progenitor scale. .. _physics-excursionSetFirstCrossingZhangHui: ``excursionSetFirstCrossingZhangHui`` ------------------------------------- An excursion set first crossing statistics class utilizing the algorithm of :cite:t:`zhang_random_2006`. First crossing (and non-crossing) rates are not supported by this method. **Methods** * ``g1`` — Returns the function :math:`g_1(S)` :cite:p:`zhang_random_2006`. * ``g2`` — Returns the function :math:`g_2(S,S^\prime)` :cite:p:`zhang_random_2006`. * ``g2Integrated`` — Returns the function :math:`g_2(S,S^\prime)` integrated over a range :math:`\Delta S` :cite:p:`zhang_random_2006`. * ``delta`` — Returns the function :math:`g_2(S,S^\prime)` integrated over a range :math:`\Delta S` :cite:p:`zhang_random_2006`. .. _physics-excursionSetFirstCrossingZhangHuiHighOrder: ``excursionSetFirstCrossingZhangHuiHighOrder`` ---------------------------------------------- An excursion set first crossing statistics class utilizing a higher order generalization of the algorithm of :cite:t:`zhang_random_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 :math:`S` space using uniform spacing in :math:`S`: .. math:: 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: .. math:: 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\}. Since :math:`g_2(S, S^\prime)` approaches infinity when :math:`S` approaches :math:`S^\prime`, one needs to define :math:`g_2(S_i, S_i)` carefully when :math:`j = i`. We can rewrite the equation: .. math:: 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). For :math:`\bar{g}_2(S_i)` we have: .. math:: \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, :math:`\delta S` depends on the range of the previous integral part. Generally, :math:`\delta S` is equal to :math:`4 \Delta S`. The above equation can be solved for :math:`f(S_i)`, giving: .. math:: 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}. Not all of the :math:`i`'s are divisible by :math:`4`. So, for the first :math:`m^\mathrm{th}` spaces, we need to calculate the integral part, separately, where :math:`m` is the remainder of :math:`i` by the modulus of :math:`4`. It is a good approximation to calculate the first part linearly. Consequently, the final formula for the general problem is: .. math:: 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}. Finally, we can solve the first-crossing distribution function, giving: .. math:: 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}, where :math:`\bar{g}_2(S_1)` is defined as: .. math:: \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. The error term for this method of discretization is: .. math:: \epsilon = \frac{(\delta S)^7}{1935360} f^{(6)}(\xi), where :math:`f^{(6)}(\xi)` is the absolute maximum of the sixth derivative in the range of :math:`[S_i , S_i + \delta S]`. For this problem, :math:`\delta S = 4 \Delta S`, so: .. math:: \epsilon \leq \frac{(\Delta S)^7}{118.125} f^{(6)}(\xi). Since there are :math:`N/4` intervals, the maximum deviation from the real value of the function is: .. math:: \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), where :math:`f^{(6)}(\xi)` is the absolute maximum of the sixth derivative in the domain, :math:`[0 , S]`. **Methods** * ``initialize`` — Initialize the high order :cite:t:`zhang_random_2006` class.