aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters/formation.py
blob: 524237c94231a6e2c353492839ed83fca04730f1 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
# Copyright (c) 2017-2018 Weitian LI <weitian@aaronly.me>
# MIT license

"""
Simulate cluster formation (i.e., merging history) using the extended
Press-Schechter formalism.

References
----------
.. [randall2002]
   Randall, Sarazin & Ricker 2002, ApJ, 577, 579
   http://adsabs.harvard.edu/abs/2002ApJ...577..579R

.. [cassano2005]
   Cassano & Brunetti 2005, MNRAS, 357, 1313
   http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
"""

import logging

import numpy as np
import scipy.integrate
import scipy.special
import scipy.optimize

from .mergertree import MergerTree
from ...share import COSMO


logger = logging.getLogger(__name__)


class ClusterFormation:
    """
    Simulate the cluster formation (i.e., merging history) using the extended
    Press-Schechter formalism by Monte Carlo methods.

    Parameters
    ----------
    M0 : float
        Cluster mass at redshift z0
        Unit: [Msun]
    z0 : float
        Redshift from where to simulate former merging history.
    zmax : float, optional
        The maximum redshift/age when to stop the formation trace.
        Default: 3.0 (i.e., looking back time ~11.5 Gyr)
    merger_mass_min : float, optional
        Minimum mass change to be regarded as a merger event instead of
        accretion.
        Unit: [Msun]

    Attributes
    ----------
    mtree : `~MergerTree`
        Merging history of this cluster.
    recent_major_merger : dict, or None
        An dictionary containing the properties of the found most recent
        major merger event, or ``None`` if not found.
    """
    def __init__(self, M0, z0, zmax=3.0, merger_mass_min=1e12):
        self.M0 = M0  # [Msun]
        self.z0 = z0
        self.zmax = zmax
        self.merger_mass_min = merger_mass_min

    @property
    def sigma_index(self):
        """
        The power-law spectral index assumed for the following density
        perturbations sigma(M).

        References: Ref.[randall2002],Eq.(2)
        """
        n = -7/5
        alpha = (n+3) / 6
        return alpha

    def f_sigma(self, mass):
        """
        Current rms density fluctuations within a sphere of specified
        mass (unit: Msun).

        It is generally sufficient to consider a power-law spectrum of
        density perturbations, which is consistent with the CDM models.

        References: Ref.[randall2002],Eq.(2)
        """
        alpha = self.sigma_index
        sigma = COSMO.sigma8 * (mass / COSMO.M8) ** (-alpha)
        return sigma

    def f_delta_c(self, z):
        """
        w = delta_c(z) is the critical linear overdensity for a region
        to collapse at redshift z.

        This is a monotone decreasing function.

        References: Ref.[randall2002],App.A,Eq.(A1)
        """
        return COSMO.overdensity_crit(z)

    def f_dw_max(self, mass):
        """
        Calculate the allowed maximum step size for tracing cluster
        formation, therefore, the adopted step size is chosen to be half
        of this maximum value.

        dw^2 ~< abs(d(ln(sigma(M)^2)) / d(ln(M))) * (dMc / M) * sigma(M)^2
              = 2 * alpha * sigma(M)^2 * dMc / M

        References: Ref.[randall2002],Sec.(3.1),Para.(1)
        """
        alpha = self.sigma_index
        dMc = self.merger_mass_min
        return np.sqrt(2 * alpha * self.f_sigma(mass)**2 * dMc / mass)

    def calc_z(self, delta_c):
        """
        Solve the redshift from the specified delta_c (a.k.a. w).
        """
        z = scipy.optimize.newton(
            lambda x: self.f_delta_c(x) - delta_c,
            x0=0, tol=1e-5)
        return z

    def calc_mass(self, S):
        """
        Calculate the mass corresponding to the given S.

        S = sigma(M)^2

        References: Ref.[randall2002],Sec.(3)
        """
        alpha = self.sigma_index
        mass = COSMO.M8 * (S / COSMO.sigma8**2)**(-1/(2*alpha))
        return mass

    @staticmethod
    def cdf_K(dS, dw):
        """
        The cumulative probability distribution function of sub-cluster
        masses.

        References: Ref.[randall2002],Eq.(5)
        """
        p = scipy.special.erfc(dw / np.sqrt(2*dS))
        return p

    @staticmethod
    def cdf_K_inv(p, dw):
        """
        Inverse function of the above ``cdf_K()``.
        """
        dS = 0.5 * (dw / scipy.special.erfcinv(p))**2
        return dS

    def gen_dS(self, dw, size=None):
        """
        Randomly generate values of dS by sampling the CDF ``cdf_K()``.
        """
        r = np.random.uniform(size=size)
        dS = self.cdf_K_inv(r, dw)
        return dS

    def simulate_mtree(self, main_only=True):
        """
        Simulate the merger tree of this cluster by tracing its formation
        using the PS formalism.

        Parameters
        ----------
        main_only : bool, optional
            Whether to only trace the forming history of the main
            halo/cluster.
            (default: True)

        References: Ref.[randall2002],Sec.(3.1)
        """
        logger.debug("Simulating cluster formation: " +
                     "M0=%.3e[Msun] " % self.M0 +
                     "from z={:.3f} back to z={zmax} ...".format(
                         self.z0, zmax=self.zmax))
        self.main_only = main_only
        if main_only:
            logger.debug("Only trace the formation of the *main* cluster ...")
            self.mtree = self._trace_main()
        else:
            logger.debug("Trace formations of both main and sub cluster ...")
            self.mtree = self._trace_formation(self.M0, _z=self.z0)
        logger.debug("Simulated cluster formation with merger tree.")
        return self.mtree

    def recent_major_merger(self, mtree=None, ratio_major=3.0):
        """
        Identify and return the most recent major merger event.

        Parameters
        ----------
        ratio_major : float, optional
            The mass ratio of the main and sub clusters to define whether
            the merger is a major event or a  minor one.
            If ``M_main/M_sub < ratio_major``, then it is a major merger.
            Default: 3.0

        Returns
        -------
        event : dict
            A dictionary with the properties of the found major event:
            ``{"M_main": M_main, "M_sub": M_sub, "R_mass": R_mass,
               "z": z, "age": age}``;
            ``None`` if no major event found.
        """
        if mtree is None:
            mtree = self.mtree

        for main, sub in mtree.itermain():
            if main["mass"] <= sub.get("mass", 0) * ratio_major:
                event = {"M_main": main["mass"],
                         "M_sub": sub["mass"],
                         "R_mass": main["mass"] / sub["mass"],
                         "z": main["z"],
                         "age": main["age"]}
                return event
        return None

    def maximum_merger(self, mtree=None):
        """
        The merger event corresponding to the biggest sub cluster, i.e.,
        the main cluster gains the most mass.

        NOTE
        ----
        Sometimes, the maximum merger event found here is not an major
        merger event.

        Returns
        -------
        event : dict
            Same as the above ``self.recent_major_event``.
            ``None`` if no mergers occurred during the traced period.
        """
        if mtree is None:
            mtree = self.mtree

        event_max = {"M_main": 0, "M_sub": 0, "R_mass": 0, "z": -1, "age": -1}
        for main, sub in mtree.itermain():
            if sub.get("mass", -1) > event_max["M_sub"]:
                event_max = {"M_main": main["mass"],
                             "M_sub": sub["mass"],
                             "R_mass": main["mass"] / sub["mass"],
                             "z": main["z"],
                             "age": main["age"]}

        if event_max["z"] <= 0:
            logger.warning("No mergers occurred at all!!")
            return None
        else:
            return event_max

    def _trace_main(self):
        """
        Iteratively trace the merger and accretion events of the
        main cluster/halo.
        """
        # Initial properties
        zc = self.z0
        Mc = self.M0
        mtree_root = MergerTree(data={"mass": Mc,
                                      "z": zc,
                                      "age": COSMO.age(zc)})
        logger.debug("[main] z=%.4f : mass=%g [Msun]" % (zc, Mc))

        mtree = mtree_root
        while True:
            # Whether to stop the trace
            if self.zmax is not None and zc > self.zmax:
                break
            if Mc <= self.merger_mass_min:
                break

            # Trace the formation by simulate a merger/accretion event
            # Notation: progenitor (*1) -> current (*2)

            # Current properties
            w2 = self.f_delta_c(z=zc)
            S2 = self.f_sigma(Mc) ** 2
            dw = 0.5 * self.f_dw_max(Mc)
            dS = self.gen_dS(dw)
            # Progenitor properties
            z1 = self.calc_z(w2 + dw)
            age1 = COSMO.age(z1)
            S1 = S2 + dS
            M1 = self.calc_mass(S1)
            dM = Mc - M1

            M_min = min(M1, dM)
            if M_min <= self.merger_mass_min:
                # Accretion
                M_main = Mc - M_min
                # NOTE: no sub node
            else:
                # Merger event
                M_main = max(M1, dM)
                M_sub = M_min
                mtree.sub = MergerTree(data={"mass": M_sub,
                                             "z": z1,
                                             "age": age1})
                logger.debug("[sub] z=%.4f : mass=%g [Msun]" % (z1, M_sub))

            # Update main cluster
            mtree.main = MergerTree(data={"mass": M_main,
                                          "z": z1,
                                          "age": age1})
            logger.debug("[main] z=%.4f : mass=%g [Msun]" % (z1, M_main))

            # Update for next iteration
            Mc = M_main
            zc = z1
            mtree = mtree.main

        return mtree_root

    def _trace_formation(self, M, _z=None, zmax=None):
        """
        Recursively trace the cluster formation and thus simulate its
        merger tree.
        """
        z = 0.0 if _z is None else _z
        node_data = {"mass": M, "z": z, "age": COSMO.age(z)}

        # Whether to stop the trace
        if self.zmax is not None and z > self.zmax:
            return MergerTree(data=node_data)
        if M <= self.merger_mass_min:
            return MergerTree(data=node_data)

        # Trace the formation by simulate a merger/accretion event
        # Notation: progenitor (*1) -> current (*2)

        # Current properties
        w2 = self.f_delta_c(z=z)
        S2 = self.f_sigma(M) ** 2
        dw = 0.5 * self.f_dw_max(M)
        dS = self.gen_dS(dw)
        # Progenitor properties
        z1 = self.calc_z(w2 + dw)
        S1 = S2 + dS
        M1 = self.calc_mass(S1)
        dM = M - M1

        M_min = min(M1, dM)
        if M_min <= self.merger_mass_min:
            # Accretion
            M_new = M - M_min
            return MergerTree(
                data=node_data,
                main=self._trace_formation(M_new, _z=z1),
                sub=None
            )
        else:
            # Merger event
            M_main = max(M1, dM)
            M_sub = M_min
            return MergerTree(
                data=node_data,
                main=self._trace_formation(M_main, _z=z1),
                sub=self._trace_formation(M_sub, _z=z1)
            )