Added the third-order bias parameter from the peak-background split

Author: komatsu5147

README.md

@@ -45,11 +45,11 @@ These multiplicity functions are assumed to be "universal", in the sense that th
 
 ## Halo bias parameters from the peak-background split approximation
 
-The halo bias parameters can be calculated from a halo multiplicity function using the so-called peak-background split (PBS) approximation. See [Desjacques, Jeong & Schmdit, Phys. Rept., 733, 1 (2018)](https://www.sciencedirect.com/science/article/pii/S0370157317304192).
+The halo bias parameters can be calculated from a halo multiplicity function using the so-called peak-background split (PBS) approximation. See Equation (2.16) and Section 3.3 of [Desjacques, Jeong & Schmdit, Phys. Rept., 733, 1 (2018)](https://www.sciencedirect.com/science/article/pii/S0370157317304192).
 
 The package contains 
 - `pbsBias(lnν, MF)` and  `pbsBias1(lnν, MF)`: The linear bias parameter from the PBS.
-- `pbsBias2(lnν, MF)`: The second-order bias parameter from the PBS.
+- `pbsBias2(lnν, MF)` and `pbsBias3(lnν, MF)`: The second- and third-order bias parameters from the PBS.
 
 Here, `MF(lnν)` is any of the halo multiplicity functions from the above list. For example, ``pbsBias(lnν, stMF)`` and ``pbsBias(lnν, x -> tinker10MF(x, z, Δm))``.
 

src/HaloMF.jl

@@ -5,7 +5,7 @@ export bocquetMFhy, bocquetMFdm
 export psMF, stMF, jenkinsMF
 export stBias, stBias1, stBias2, stBias3
 export tinker10Bias
-export pbsBias, pbsBias1, pbsBias2
+export pbsBias, pbsBias1, pbsBias2, pbsBias3
 include("tinkerMF.jl")
 include("bocquetMF.jl")
 include("classicMF.jl")

src/pbsBias.jl

@@ -4,7 +4,7 @@
 
 Linear halo bias, b1(lnν), from the peak-background split (PBS) approximation.
 
-*Reference*: Section 3.3 of Desjacques, Jeong & Schmidt, Phys. Rept., 733, 1 (2018)
+*Reference*: Equation (2.14) and Section 3.3 of Desjacques, Jeong & Schmidt, Phys. Rept., 733, 1 (2018)
 
 # Arguments
 - `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R,z)]^2. Here, δc = 1.6865 and σ(R,z) is the r.m.s. mass fluctuation within
@@ -17,7 +17,8 @@ The linear PBS bias is normalized as
 """
 function pbsBias(lnν::Real, MF)
     δc = 1.6865
-    b1 = 1 - ForwardDiff.derivative(MF, lnν) * 2 / δc / MF(lnν)
+    b1L = -ForwardDiff.derivative(MF, lnν) * 2 / δc / MF(lnν)
+    b1 = 1 + b1L
     return b1
 end
 const pbsBias1 = pbsBias
@@ -27,7 +28,7 @@ const pbsBias1 = pbsBias
 
 Second-order PBS halo bias, b2(lnν), from the peak-background split (PBS) approximation.
 
-*Reference*: Section 3.3 of Desjacques, Jeong & Schmidt, Phys. Rept., 733, 1 (2018)
+*Reference*: Equation (2.15) and Section 3.3 of Desjacques, Jeong & Schmidt, Phys. Rept., 733, 1 (2018)
 
 # Arguments
 - `lnν::Real`: natural logarithm of a threshold, ν (see `pbsBias`).
@@ -41,10 +42,38 @@ function pbsBias2(lnν::Real, MF)
     δc = 1.6865
     a2 = -17 / 21
     dMFdlnν(x) = ForwardDiff.derivative(MF, x)
-    b2 = (
-        -2 * (1 + a2) * dMFdlnν(lnν) * 2 / δc / MF(lnν) +
-        ForwardDiff.derivative(dMFdlnν, lnν) * 4 / δc^2 / MF(lnν) -
-        dMFdlnν(lnν) * 2 / δc^2 / MF(lnν)
-    )
+    b1L = -dMFdlnν(lnν) * 2 / δc / MF(lnν)
+    b2L = (ForwardDiff.derivative(dMFdlnν, lnν) * 2 - dMFdlnν(lnν)) * 2 / δc^2 / MF(lnν)
+    b2 = 2 * (1 + a2) * b1L + b2L
     return b2
-end
+end
+
+"""
+    pbsBias3(lnν, MF)
+
+Third-order PBS halo bias, b3(lnν), from the peak-background split (PBS) approximation.
+
+*Reference*: Equation (2.16) and Section 3.3 of Desjacques, Jeong & Schmidt, Phys. Rept., 733, 1 (2018)
+
+# Arguments
+- `lnν::Real`: natural logarithm of a threshold, ν (see `pbsBias`).
+- `MF`(lnν): a function which returns a halo multiplicity function with the argument lnν.
+
+The third-order PBS bias satisfies
+
+``∫_-∞^∞ dlnν MF(lnν) pbsBias3(MF) = 0``
+"""
+function pbsBias3(lnν::Real, MF)
+    δc = 1.6865
+    a2, a3 = -17 / 21, 2815 / 3969
+    dMFdlnν(x) = ForwardDiff.derivative(MF, x)
+    d2MFdlnν2(x) = ForwardDiff.derivative(dMFdlnν, x)
+    b1L = -dMFdlnν(lnν) * 2 / δc / MF(lnν)
+    b2L = (d2MFdlnν2(lnν) * 2 - dMFdlnν(lnν)) * 2 / δc^2 / MF(lnν)
+    b3L = -(
+        (ForwardDiff.derivative(d2MFdlnν2, lnν) * 2 - d2MFdlnν2(lnν) * 3 + dMFdlnν(lnν)) *
+        4 / δc^3 / MF(lnν)
+    )
+    b3 = 6 * (a2 + a3) * b1L + 3 * (1 + 2 * a2) * b2L + b3L
+    return b3
+end

src/stBias.jl

@@ -77,7 +77,7 @@ function stBias3(lnν::Real)
     ν = exp(lnν)
     qν = _st_q * ν
     a2 = -17 / 21
-    a3 = 341 / 567
+    a3 = 2815 / 3969
     ϵ1 = (qν - 1) / _st_δc
     ϵ2 = qν * (qν - 3) / _st_δc^2
     ϵ3 = qν * (qν^2 - 6 * qν + 3) / _st_δc^3