Accretion of DM onto a moving black hole recent results ( P.M & A.O arXiv:2009.13433 ) Andrzej Odrzywołek

(1)

Accretion of DM onto a moving black hole recent results ( P.M & A.O arXiv:2009.13433 )

Andrzej Odrzywołek

IFT UJ

Seminarium ZTWiA 10:15, Środa, 25 listopada 2020

β=10^-2 β=1 β=10² γ

2 4 6 8 10

M(v)/M(0)

(2)

Hoyle-Lyttleton accretion (cross section)

axially symmetric, stationary „accretion” on „point mass”

bound trajectories are assumed to be trapped v∞2

2 <GM R

where R is the Hoyle-Lyttleton radius (impact parameter) accretion via 1D column (black arrow) with rate ˙M

M = πR˙ ²v∞ρ∞= 4πG²M² v∞3 ρ∞.

Seminarium ZTWiA

(3)

Hoyle-Lyttleton accretion (comment)

M =˙ 4πG²M²

v∞3 ρ∞→ 4πG²M²

p

v∞2+ cs23ρ∞¬ 4πG²M² c³ ρ∞

model is nonsense, divergent at v∞→ 0

M can be regularized by „interpolation” between Bondi (gas with sound speed˙ cs) and H-L accretion formulas

two classes of trajectories:

(i) scattered (parabolic & hyperbolic) (ii) trapped in orbit (elliptical)

in Newton theory there are no true accretted trajectories (except direct hit at point mass)

above formulas are surprisingly resistant to time, frequently verified by numerical/Monte Carlo simulations

this is probably good example of dimensional analysis, giving little freedom on possible ˙M formulas

(4)

Vlasov gas model

parallel stream of uniform velocity particles is too simple approximation (ignore random/thermal motion)

motion in black hole (Schwarzschild) metric provides missing class of absorbed particles (see e.g: D. Hunik, 2014, praca lic.)

random motions can be incorporated using Vlasov gas formalism (earlier talks by P. Mach & A. Cieślik)

probability density distribution function f (x^µ, pν) is the solution of Vlasov equation:

∂H

∂pµ

∂f

∂x^µ− ∂H

∂x^ν

∂f

∂pν

= {H, f } = 0, H(x^µ, pν) =1

2g^µν(x^α)pµpν

ingenuious trick of Rioseco & Sarbach (2017) is to use canonical transformation using action-angle variables making above equation trivial

boundary (asymptotic) conditions are relativistic generalization of canonical ensemble (Maxwell-J¨uttner distribution) Lorentz - boosted with velocity v∞:

f (x^µ, pν) ∝ exp

_βp

t

m

→ α exp

n

_β

mγ

h

pt− v∞

cos θpr−sin θ r pθ

io

where β =^mc_kT²

then, most of the effort is to use various variables and use them to compute phase-space integrals . . . done by Patryk Mach

Seminarium ZTWiA

(5)

Vlasov gas model

parallel stream of uniform velocity particles is too simple approximation (ignore random/thermal motion)

motion in black hole (Schwarzschild) metric provides missing class of absorbed particles (see e.g: D. Hunik, 2014, praca lic.)

random motions can be incorporated using Vlasov gas formalism (earlier talks by P. Mach & A. Cieślik)

probability density distribution function f (x^µ, pν) is the solution of Vlasov equation:

∂H

∂pµ

∂f

∂x^µ− ∂H

∂x^ν

∂f

∂pν

= {H, f } = 0, H(x^µ, pν) =1

2g^µν(x^α)pµpν

ingenuious trick of Rioseco & Sarbach (2017) is to use canonical transformation using action-angle variables making above equation trivial

boundary (asymptotic) conditions are relativistic generalization of canonical ensemble (Maxwell-J¨uttner distribution) Lorentz - boosted with velocity v∞:

f (x^µ, pν) ∝ exp

_βp

t

m

→ α exp

n

_β

mγ

h

pt− v∞

cos θpr−sin θ r pθ

io

(6)

Random motion effect on Hoyle-Lyttleton accretion

Seminarium ZTWiA

(7)

Parallel stream of Schwarzschild geodesics

Trajectory classification

1 BLUE - accreted by a black hole

(8)

Parallel stream of Schwarzschild geodesics

Trajectory classification

1 BLUE - accreted by a black hole

2 GREEN - trapped in orbit

3 RED - scattered to infinity

Seminarium ZTWiA

(9)

0

dm pµf m³λ ξ²

p

ε²− U_λ(ξ)

Jµ^(scat)=

0

dm pµf m³λ ξ²

p

ε²− U_λ(ξ)

pµ= (pt, pr, pθ, pϕ) = m −ε,(1 − N η)ε + r

p

ε²− U_λ(ξ)

N , M λ cos χ, M λ sin θ sin χ

!

f ∝ F = α exp

n

βγ

h

−ε − v

p

ε²− 1 (rcos X(ξ, ε, λ) cos θ − sin X(ξ, ε, λ) sin θ cos χ)

io

(10)

Currents and density: ready to numerical evaluation

Jµ^(abs)= 2παm⁴ ξ²

∞

Z

1

dεe^−βγε

λ_c(ε)

Z

0

dλ λ exp(Y )

p

ε²− Uλ(ξ)







−εI₀(Z)

−εI0(Z) +

1+^λ2 ξ2 ε+√

ε²−U_λ(ξ)I0(Z) M λI1(Z)

0







Jµ^(scat)=4παm⁴ ξ²

∞

Z

ε_min(ξ)

dεe^−βγε

λmax(ξ,ε)

Z

λ_c(ε)

dλ λ cosh (Y )

p

ε²− U_λ(ξ)







−εI0(Z)

I₀(Z) N

(1 − N )ε −

p

ε²− U_λtgh Y

−M λI1(Z) 0







X(ξ, ε, λ) = λ

∞

Z

ξ

dξ⁰ ξ⁰²

p

ε²− Uλ(ξ⁰)

Y = βγv

p

ε²− 1 cos X(ξ, ε, λ) cos θ

Z = βγv

p

ε²− 1 sin X(ξ, ε, λ) sin θ

Surprisingly, above appeared to be impossible to evaluate numerically in reliable way . . .

Seminarium ZTWiA

(11)

Formula for mass accretion rate ˙ M

Fortunately, simple formula can be derived

M =˙ 4πG²M²ρ∞

c³

β 4K2(β)×

∞

Z

1

dε e^−βγελ²_c(ε)sinh βγv/c√ ε²− 1

βγv/c√ ε²− 1

where leading term is Hoyle-Lyttleton formula evaluated at v∞= c

MH-L(v = c) =4πG²M²ρ∞

c³ .

Numerical problems

Numerical evaluation of the above integral produced qualitatively contradictory results without any error messages (in Mathematica) depending on methods used.

Problem was pinpointed to a function λc(critical angular momentum)

λc(ε) =

s

12

1 − 4 δ −√

δ²+ δ

² ⁼

r

12

1 − 4δ − 8δ²+ 8δ√ δ²+ δ,

(12)

Numerical instability

Proper formula

Seminarium ZTWiA

(13)

Numerical instability

Proper formula

(14)

Numerical instability

Proper formula

λ²_c=



 

 

 

 

12 for ε =

p

8/9

12

1− ⁴

√3ε 9ε2 −8

+1

2

for

p

8/9 < ε < 10⁴

27ε²− 9 for ε 10⁴

Seminarium ZTWiA

(15)

Elliptic & inverse Weierstrass ℘ functions

Time-critical function

X(ξ, ε, λ) = λ

∞

Z

ξ

dξ⁰

ξ⁰²

r

ε²− 1 − ²

ξ⁰

1 + ^λ²

ξ^{0 2}

where dimensionless radius, energy and angular momentum are

r = M ξ, E = mε, l = M mλ

and roots of expression under√

classify orbits.

There are at least 3 methods to compute X:

1 direct numerical integration (slow, numerically clean)

2 use of inverse ℘ (fastest, numerically beyond control)

3 use of elliptic integrals and roots (obscure formulas, but working)

(16)

Inverse Weierstrass ℘ functions

X = −√ 2a

−^a

Z

6

a ξ−^a

6

dw

p

4w³− g2w − g3

g2= −2 λ²a²

2 −λ² 6

, g3= −2 λ²a³

ε²−2 3−λ²

54

z =

∞

Z

℘(z;g₂,g₃)

dw

p

4w³− g₂w − g3

Mathematica implementation

Above formulas can be computed very fast using Mathematica function InverseWeierstrassP

InverseWeierstrassP

InverseWeierstrassP. Unfortunately, in large-scale numerical calculations unstable behavior was observed, and location of 2D complex lattice out-of-control.

Seminarium ZTWiA

(17)

Elliptic functions

One can compute indefinite integral

I =

Z

dw

p

4w³−g2w−g3

=

2 (w − w3)

q

w−w1 w₃−w₁

q

w−w2 w₃−w₂

q

w−w₃ w2−w₃

p

4w³− g2w − g3

F

arc sin

r

w3−w w3−w2

|w2− w3

w1− w3

+C

where w1, w2, w3denote (possibly complex) roots of the polynomial

4w³− g2w − g3= 0, and F (φ|m) is Legendre’s elliptic integral of the first kind. X is then computed as

X = −√

2a [I(−a/6) − I(a/ξ − a/6)] .

Mathematica implementation

Above formulas can be computed, although slowly, using Mathematica functions EllipticF

EllipticF

EllipticF and RootRootRoot. To recall, in every point of 2D space, we must compute numerically six double integrals involving X.

(18)

M >˙ 27G²M² c³ mn∞γ

4 for ultra-relativistic velocities flow is complicated and a bit unexpected

5 computation of the X(ξ, ε, λ) elliptic/Weierstrass ℘ related function is now time-critical operation; possible room for performance improvement (now days for 100-core machine)

6 moving black hole might be subject of the runaway accretion with M → ∞ in finite time — possible growth of primordial black holes

7 . . .

Seminarium ZTWiA

(53)

(54)

v = 0.95c, mc²/kT = 1

Seminarium ZTWiA

(69)

Density profiles (polar)

v = 0.99c, mc²/kT = 1

Accretion of DM onto a moving black hole recent results ( P.M & A.O arXiv:2009.13433 ) Andrzej Odrzywołek