1 """Routines related to the reionization history of the IGM."""
2
3 import math
4
5 import numpy
6 import scipy
7 import scipy.integrate as si
8
9 import perturbation as cp
10 import distance as cd
11 import constants as cc
12 import density as cden
13 import utils as cu
14
16 """The Lyman-alpha wavelength shift given light-travel distance.
17
18 Wavelengths are in Angstroms.
19
20 Returns lambda(z), lambda(z - Deltaz), z, z - Deltaz
21
22 """
23
24 dl = cd.light_travel_distance(z, **cosmo)[0]
25 dl0 = dl - delta_dl
26 z0 = cd.redshift_d_light(dl0, **cosmo)
27
28 wavelength = cc.lambda_Lya_0 * (1+z)
29 wavelength0 = cc.lambda_Lya_0 * (1+z0)
30 return wavelength, wavelength0, z, z0
31
33 """Recombination rate coefficients for HII, HeII and HeIII.
34
35 Parameters
36 ----------
37
38 temp is the temperature in K
39
40 species is 'H', 'He0', or 'He1'.
41
42 case is 'A' or 'B'.
43
44 Notes
45 -----
46
47 From Hui and Gnedin (1997MNRAS.292...27H).
48
49 Valid for He0 for temperatures between 5e3 and 5e5 K.
50
51 """
52
53 if not((species == 'H') or (species == 'He0') or (species == 'He1')):
54 raise(exception.ValueError("Bad species: " + species))
55
56 if case == 'A':
57 case_N = 0
58 elif case == 'B':
59 case_N = 1
60 else:
61 raise(exception.ValueError("Case must be 'A' or 'B': " + case))
62
63
64 T_TR = {'H' : 157807.,
65 'He0' : 285335.,
66 'He1' : 631515.
67 }
68
69
70 a = {'H' : [1.269e-13, 2.753e-14],
71 'He0' : [3e-14, 1.26e-14],
72 'He1' : [2. * 1.269e-13, 2. * 2.753e-14]
73 }
74
75 p0 = {'H' : [1.503, 1.500],
76 'He0' : [0.654, 0.750],
77 'He1' : [1.503, 1.500]
78 }
79
80 p1 = {'H' : [0.470, 0.407],
81 'He0' : [0, 0],
82 'He1' : [0.470, 0.407]
83 }
84
85 p2 = {'H' : [1.923,2.242],
86 'He0' : [0, 0],
87 'He1' : [1.923,2.242]
88 }
89
90 con = {'H' : [0.522, 2.740],
91 'He0' : [1, 1],
92 'He1' : [0.522, 2.740],
93 }
94
95 lam = 2. * T_TR[species] / temp
96 return (a[species][case_N] * lam ** p0[species][case_N] /
97 (1.0 + (lam / con[species][case_N]) ** p1[species][case_N])
98 ** p2[species][case_N])
99
102 """The ionized fraction of the universe using perturbation theory.
103
104 Parameters
105 ----------
106
107 z:
108
109 Redshift values at which to evaluate the ionized fraction.
110
111 coeff_ion:
112
113 Coefficient giving the ratio between collapse fraction and
114 ionized fraction (neglecting recombinations and assuming all
115 photons are instantly absorbed).
116
117 temp_min:
118
119 Either the minimum Virial temperature (in Kelvin) or minimum
120 mass of halos (in solar masses) contributing to reionization.
121
122 passed_temp_min: Boolean
123
124 Set this to True if you pass a minimum mass, False (default) if
125 you pass a minimum Virial temperature.
126
127 cosmo: dict
128
129 Cosmological parameters.
130
131 Notes
132 -----
133
134 See Furlanetto et al. (2004ApJ...613....1F).
135
136 """
137 sd = cp.sig_del(temp_min, z, passed_min_mass=passed_min_mass, **cosmo)
138 cf = cp.collapse_fraction(*sd)
139 w = cf * coeff_ion
140 return w
141
142 -def quick_ion_col_function(coeff_ion, temp_min, passed_min_mass = False,
143 zmax = 20., zmin = 0., zstep = 0.1, **cosmo):
144 """Return a function giving ionization_from_collapse as a function
145 of redshift (based on interpolation).
146
147 Calling the resulting function is much faster than evaluating
148 ionization_from_collapse.
149 """
150 z = numpy.arange(zmin, zmax, zstep)
151 w = ionization_from_collapse(z, coeff_ion, temp_min,
152 passed_min_mass, **cosmo)
153 return scipy.interpolate.interp1d(z, w)
154
156 """Clumping factor as a function of redshift used by Bagla et al. 2009.
157
158 See Bagla, Kulkarni & Padmanabhan (2009MNRAS.397..971B).
159 """
160 return numpy.sqrt(26.2917 * numpy.exp(-0.1822 * z + 0.003505 * z**2.))
161
163 """Clumping factor as a function of redshift used by Haiman & Bryan (2006).
164
165 See Haiman & Bryan (2006ApJ...650....7H).
166 """
167 return 1 + 9 * ((1 + z)/7)**(-beta)
168
170 """Clumping factor as a function of redshift estimated from Chary (2008)
171
172 Chary, R.-R. 2008, ApJ, 680, 32 (2008ApJ...680...32C) shows a nice
173 plot (Figure 2a) of clumping factor for neutral and ionized gas
174 with and without halos included and adopts the clumping factor for
175 ionized gas without source halos (but with other halos), which
176 rises (apparently, from the graph) as a constant powerlaw from ~2
177 and z=15 to 6 at z=8, steepens to reach 8 at z=7, and ~17 at
178 z=5.
179
180 This function returns the values of a piecewise powerlaw (as a
181 function of redshift) interpolated/extrapolated through those
182 points.
183 """
184 _zclumpChary = numpy.array([15, 8, 7, 5])
185 _clumpChary = numpy.array([2, 6, 8, 17])
186 _logclumpChary = numpy.log10(_clumpChary)
187 _logczfunc = cu.Extrapolate1d(_zclumpChary, _logclumpChary,
188 bounds_behavior=['extrapolate',
189 'extrapolate'],
190 slopes=[None, None],
191 npoints = [2, 2])
192
193 return 10.0**_logczfunc(z)
194
195 -def _udot(u, t, coeff_rec_func, redshift_func, ion_func, bubble=True):
196 """du/dt where u = x - f_* f_esc,gamma N_gamma F
197
198 Parameters
199 ----------
200
201 u: integral of du/dt as defined below
202
203 t: cosmic age in s
204
205 redshift_func: function returning redshift given t
206
207 ion_func: function returing ionized fraction neglecting recombinations
208
209 coeff_rec_func: function returning clumping_factor alpha_B n_H_0 (1+z)^3
210
211 bubble: If True, assume ionized gas is in fully-ionized bubbles
212 and other gas is fully neutral. If False, asssume gas is
213 uniformly fractionally ionized.
214
215 Notes
216 -----
217
218 This is implemented as a reformulation of the normal ODE
219 describing ionization and recombination (see, e.g. Bagla, Kulkarni
220 & Padmanabhan (2009MNRAS.397..971B).
221
222 The original ODE is:
223
224 dx/dt = -alpha_B C n_H x + f_* f_esc,gamma N_gamma dF/dt
225
226 If we let u = x - w, where w = f_* f_esc,gamma N_gamma F(t) then
227
228 du/dt = dx/dt - dw/dt
229
230 which gives
231
232 du/dt = -alpha_B C n_H x = -alpha_B C n_H (u + w)
233
234 We have an analytical expression for w, so we can numerically
235 integrate the ODE to give us u(t) or x(t) = u(t) + w(t).
236
237 """
238 z = redshift_func(t)
239 w = ion_func(z)
240 crf = coeff_rec_func(z)
241
242
243
244
245 x = u + w
246 x = x * (x <= 1.) + 1.0 * (x > 1.)
247 if bubble:
248 udot = -1. * crf * x
249 else:
250 udot = -1. * crf * x**2
251
252 if (False):
253 print ("z=%.3f; t=%.1g; c=%.2g; udot=%.2g; w,x,u = %.2g, %.2g, %.2g" %
254 (z, t, crf, udot, w, x, u))
255 return udot
256
257 -def integrate_ion_recomb(z,
258 ion_func,
259 clump_fact_func,
260 xHe=1.0,
261 temp_gas=1e4,
262 alpha_B=None,
263 bubble=True,
264 **cosmo):
265 """Integrate IGM ionization and recombination given an ionization function.
266
267 Parameters:
268
269 z: array
270
271 The redshift values at which to calculate the ionized
272 fraction. This array should be in reverse numerical order. The
273 first redshift specified should be early enough that the
274 universe is still completely neutral.
275
276 ion_func:
277
278 A function giving the ratio of the total density of emitted
279 ionizing photons to the density hydrogen atoms (or hydrogen
280 plus helium, if you prefer) as a function of redshift.
281
282 temp_gas:
283
284 Gas temperature used to calculate the recombination coefficient
285 if alpha_b is not specified.
286
287 alpha_B:
288
289 Optional recombination coefficient in units of cm^3
290 s^-1. In alpha_B=None, it is calculated from temp_gas.
291
292 clump_fact_func: function
293
294 Function returning the clumping factor when given a redshift,
295 defined as <n_HII^2>/<n_HII>^2.
296
297 cosmo: dict
298
299 Dictionary specifying the cosmological parameters.
300
301 Notes:
302
303 We only track recombination of hydrogen, but if xHe > 0, then the
304 density is boosted by the addition of xHe * nHe. This is
305 eqiuvalent to assuming the the ionized fraction of helium is
306 always proportional to the ionized fraction of hydrogen. If
307 xHe=1.0, then helium is singly ionized in the same proportion as
308 hydrogen. If xHe=2.0, then helium is fully ionized in the same
309 proportion as hydrogen.
310
311 We assume, as is fairly standard, that the ionized
312 fraction is contained in fully ionized bubbles surrounded by a
313 fully neutral IGM. The output is therefore the volume filling
314 factor of ionized regions, not the ionized fraction of a
315 uniformly-ionized IGM.
316
317 I have also made the standard assumption that all ionized photons
318 are immediately absorbed, which allows the two differential
319 equations (one for ionization-recombination and one for
320 emission-photoionizaion) to be combined into a single ODE.
321
322 """
323
324
325 if alpha_B is None:
326 alpha_B_cm = recomb_rate_coeff_HG(temp_gas, 'H', 'B')
327 else:
328 alpha_B_cm = alpha_B
329 alpha_B = alpha_B_cm * cc.Gyr_s / (cc.Mpc_cm**3.)
330 print ("Recombination rate alpha_B = %.4g (Mpc^3 Gyr^-1) = %.4g (cm^3 s^-1)"
331 % (alpha_B, alpha_B_cm))
332
333
334 if 'deltaSqr' not in cosmo:
335 cosmo['deltaSqr'] = cp.norm_power(**cosmo)
336
337
338 rho_crit, rho_0, n_He_0, n_H_0 = cden.baryon_densities(**cosmo)
339
340
341 nn = (n_H_0 + xHe * n_He_0)
342
343
344
345 coeff_rec_func = lambda z1: (clump_fact_func(z1) *
346 alpha_B *
347 nn * (1.+z1)**3.)
348
349
350 red_func = cd.quick_redshift_age_function(zmax = 1.1 * numpy.max(z),
351 zmin = -0.0,
352 dz = 0.01,
353 **cosmo)
354
355 ref_func_Gyr = lambda t1: red_func(t1 * cc.Gyr_s)
356
357
358 t = cd.age(z, **cosmo)[0]/cc.Gyr_s
359
360
361 u = si.odeint(_udot, y0=0.0, t=t,
362 args=(coeff_rec_func, ref_func_Gyr, ion_func, bubble))
363 u = u.flatten()
364
365 w = ion_func(z)
366 x = u + w
367
368 return x, w, t
369
377
378 """IGM ionization state with recombinations from halo collapse
379 fraction. Integrates an ODE describing IGM ionization and
380 recombination rates.
381
382 z: array
383
384 The redshift values at which to calculate the ionized
385 fraction. This array should be in reverse numerical order. The
386 first redshift specified should be early enough that the
387 universe is still completely neutral.
388
389 coeff_ion:
390
391 The coefficient converting the collapse fraction to ionized
392 fraction, neglecting recombinations. Equivalent to the product
393 (f_star * f_esc_gamma * N_gamma) in the BKP paper.
394
395
396 temp_min:
397
398 See docs for ionization_from_collapse. Either the minimum virial
399 temperature or minimum mass of halos contributing to
400 reionization.
401
402 passed_temp_min:
403
404 See documentation for ionization_from_collapse.
405
406 temp_gas:
407
408 Gas temperature used to calculate the recombination coefficient
409 if alpha_b is not specified.
410
411 alpha_B:
412
413 Optional recombination coefficient in units of cm^3
414 s^-1. In alpha_B=None, it is calculated from temp_gas.
415
416 clump_fact_func: function
417
418 Function returning the clumping factor when given a redshift.
419
420 cosmo: dict
421
422 Dictionary specifying the cosmological parameters.
423
424 We assume, as is fairly standard, that the ionized
425 fraction is contained in fully ionized bubbles surrounded by a
426 fully neutral IGM. The output is therefore the volume filling
427 factor of ionized regions, not the ionized fraction of a
428 uniformly-ionized IGM.
429
430 I have also made the standard assumption that all ionized photons
431 are immediately absorbed, which allows the two differential
432 equations (one for ionization-recombination and one for
433 emission-photoionizaion) to be combined into a single ODE.
434
435 """
436
437
438 if alpha_B is None:
439 alpha_B_cm = recomb_rate_coeff_HG(temp_gas, 'H', 'B')
440 else:
441 alpha_B_cm = alpha_B
442 alpha_B = alpha_B_cm / (cc.Mpc_cm**3.)
443 print ("Recombination rate alpha_B = %.4g (Mpc^3 s^-1) = %.4g (cm^3 s^-1)"
444 % (alpha_B, alpha_B_cm))
445
446
447 if 'deltaSqr' not in cosmo:
448 cosmo['deltaSqr'] = cp.norm_power(**cosmo)
449
450
451 rho_crit, rho_0, n_He_0, n_H_0 = cden.baryon_densities(**cosmo)
452
453
454
455 coeff_rec_func = lambda z: (clump_fact_func(z)**2. *
456 alpha_B *
457 n_H_0 * (1.+z)**3.)
458
459
460 redfunc = cd.quick_redshift_age_function(zmax = 1.1 * numpy.max(z),
461 zmin = -0.0,
462 **cosmo)
463
464 ionfunc = quick_ion_col_function(coeff_ion,
465 temp_min,
466 passed_min_mass = passed_min_mass,
467 zmax = 1.1 * numpy.max(z),
468 zmin = -0.0,
469 zstep = 0.1, **cosmo)
470
471
472 t = cd.age(z, **cosmo)
473
474
475 u = si.odeint(_udot, y0=0.0, t=t,
476 args=(coeff_rec_func, redfunc, ionfunc))
477 u = u.flatten()
478
479 w = ionization_from_collapse(z, coeff_ion, temp_min,
480 passed_min_mass = passed_min_mass,
481 **cosmo)
482 x = u + w
483 x[x > 1.0] = 1.0
484 return x, w, t
485
486 -def ionization_from_luminosity(z, ratedensityfunc, xHe=1.0,
487 rate_is_tfunc = False,
488 ratedensityfunc_args = (),
489 method = 'romberg',
490 **cosmo):
491 """Integrate the ionization history given an ionizing luminosity
492 function, ignoring recombinations.
493
494 Parameters
495 ----------
496
497 ratedensityfunc: callable
498 function giving comoving ionizing photon emission rate
499 density, or ionizing emissivity (photons s^-1 Mpc^-3) as a
500 function of redshift (or time).
501
502 rate_is_tfunc: boolean
503 Set to true if ratedensityfunc is a function of time rather than z.
504
505 Notes
506 -----
507
508 Ignores recombinations.
509
510 The ionization rate is computed as ratedensity / nn, where nn = nH
511 + xHe * nHe. So if xHe is 1.0, we are assuming that helium becomes
512 singly ionized at proportionally the same rate as hydrogen. If xHe
513 is 2.0, we are assuming helium becomes fully ionizing at
514 proportionally the same rate as hydrogen.
515
516 The returened x is therefore the ionized fraction of hydrogen, and
517 the ionized fraction of helium is xHe * x.
518
519 """
520
521 cosmo = cd.set_omega_k_0(cosmo)
522 rhoc, rho0, nHe, nH = cden.baryon_densities(**cosmo)
523 nn = (nH + xHe * nHe)
524 if rate_is_tfunc:
525 t = cd.age(z, **cosmo)[0]
526 def dx_dt(t1):
527 return numpy.nan_to_num(ratedensityfunc(t1, *ratedensityfunc_args) /
528 nn)
529 sorti = numpy.argsort(t)
530 x = numpy.empty(t.shape)
531 x[sorti] = cu.integrate_piecewise(dx_dt, t[sorti], method = method)
532 return x
533 else:
534 dt_dz = lambda z1: cd.lookback_integrand(z1, **cosmo)
535 def dx_dz(z1):
536 z1 = numpy.abs(z1)
537 return numpy.nan_to_num(dt_dz(z1) *
538 ratedensityfunc(z1, *ratedensityfunc_args) /
539 nn)
540 sorti = numpy.argsort(-z)
541 x = numpy.empty(z.shape)
542 x[sorti] = cu.integrate_piecewise(dx_dz, -z[sorti], method = method)
543 return x
544
546 """The electron scattering optical depth given ionized filling
547 factor vs. redshift.
548
549 Parameters
550 ----------
551
552 x_ionH: array
553
554 Ionized fraction of hydrogen as a function of z. Should be [0,1].
555
556 x_ionHe: array
557
558 Set x_ionHE to X_HeII + 2 * X_HeIII, where X_HeII is the
559 fraction of helium that is singly ionized, and X_HeII is the
560 fraction of helium that is doubly ionized. See Notes below.
561
562 z: array
563 Redshift values at which the filling factor is specified.
564
565 cosmo: cosmological parameters
566
567 uses: 'X_H' and/or 'Y_He', plus parameters needed for hubble_z
568
569 Returns
570 -------
571
572 tau: array
573 The optical depth as a function of z.
574
575 Notes
576 -----
577
578 The precision of your result depends on the spacing of the input
579 arrays. When in doubt, try doubling your z resolution and see if
580 the optical depth values have converged.
581
582 100% singly ionized helium means x_ionHe = 1.0, 100% doubly
583 ionized helium means x_ionHe = 2.0
584
585 If you want helium to be singly ionized at the same rate as
586 hydrogen, set x_ionHe = x_ionH.
587
588 If you want helium to be doubly ionized at the same rate as
589 hydrogen is ionized, set x_ionHe = 2 * x_ionH.
590
591 """
592
593 rho_crit, rho_0, n_He_0, n_H_0 = cden.baryon_densities(**cosmo)
594
595
596 n_p = n_H_0 + 2. * n_He_0
597
598
599 n_e = n_H_0 * x_ionH + n_He_0 * x_ionHe
600
601
602 x = n_e / n_p
603
604 H_0 = cc.H100_s * cosmo['h']
605
606
607 tau_star = cc.c_light_Mpc_s * cc.sigma_T_Mpc * n_p
608
609
610 H_z = cd.hubble_z(z, **cosmo)
611
612
613 integrand = -1. * tau_star * x * ((1. + z)**2.) / H_z
614
615 integral = numpy.empty(integrand.shape)
616 integral[...,1:] = si.cumtrapz(integrand, z)
617 integral[...,0] = 0.0
618 return numpy.abs(integral)
619
620 -def optical_depth_instant(z_r, x_ionH=1.0, x_ionHe=1.0, z_rHe = None,
621 return_tau_star=False, verbose=0, **cosmo):
622 """Optical depth assuming instantaneous reionization and a flat
623 universe.
624
625 Calculates the optical depth due to Thompson scattering off free
626 electrons in the IGM.
627
628 Parameters
629 ----------
630
631 z_r:
632 Redshift of instantaneos reionization.
633
634 x_ionH:
635 Ionized fraction of hydrogen after reionization.
636
637 x_ionHe:
638 Set to 2.0 for fully ionized helium. Set to 1.0 for singly
639 ionized helium. Set to 0.0 for neutral helium. This value
640 equals X_HeII + 2 * X_HeIII after z_r (where X_HeII is the
641 fraction of helium that is singly ionized, and X_HeII is the
642 fraction of helium that is doubly ionized).
643
644 z_rHe (optional):
645 Redshift of instantaneos Helium reionization, i.e. when helium
646 becomes doubly ionized. z_rHe should be less than z_r.
647
648 return_tau_star: Boolean
649 whether or not to return the value of tau_star, as defined by
650 Griffiths et al. (arxiv:astro-ph/9812125v3)
651
652 cosmo: cosmological parameters
653
654 Returns
655 -------
656
657 tau: array
658 optical depth to election
659
660 tau_star: array or scalar
661
662 Notes
663 -----
664
665 See, e.g. Griffiths et al. (arxiv:astro-ph/9812125v3, note that
666 the published version [ 1999MNRAS.308..854G] has typos)
667
668 """
669
670 if numpy.any(cden.get_omega_k_0(**cosmo) != 0):
671 raise ValueError, "Not valid for non-flat (omega_k_0 !=0) cosmology."
672
673
674 if z_rHe is not None:
675
676 tau_short_all = optical_depth_instant(z_rHe, x_ionH, x_ionHe=2.0,
677 **cosmo)
678
679
680 tau_short_H = optical_depth_instant(z_rHe, x_ionH, x_ionHe, **cosmo)
681
682
683 tau_short_He = tau_short_all - tau_short_H
684 if(verbose > 0) :
685 print "tau_short_He = ", tau_short_He
686
687 rho_crit, rho_0, n_He_0, n_H_0 = cden.baryon_densities(**cosmo)
688
689
690 n_p = n_H_0 + 2. * n_He_0
691
692
693 n_e = (n_H_0 * x_ionH) + (n_He_0 * x_ionHe)
694
695
696 x = n_e / n_p
697
698 if(verbose > 0) :
699 print "n_He/n_H = ", n_He_0 / n_H_0
700 print "x = ne/np = ", x
701 print "n_e/n_H_0 = ", n_e/n_H_0
702
703 H_0 = cc.H100_s * cosmo['h']
704
705
706 tau_star = cc.c_light_Mpc_s * cc.sigma_T_Mpc * n_p * x / H_0
707
708
709
710
711 e_z_reion = cd.e_z(z_r, **cosmo)
712
713 tau = 2. * tau_star * (e_z_reion - 1.0) / (3. * cosmo['omega_M_0'])
714
715 if z_rHe is not None:
716
717 tau += tau_short_He
718
719 if return_tau_star:
720 return tau, tau_star
721 else:
722 return tau
723
725 """Recombination rate density from Madau, Haardt, & Rees 1999.
726
727 Assumes hydrogen is fully ionized.
728
729 Units are s^-1 coMpc^-3.
730
731 """
732 return 1e50 * clump_fact * ((1. + z)/7.)**3
733