Mathematica Notebook for "Exact, non-singular black holes from a phantom DBI Field as primordial dark matter" https://arxiv.org/abs/2511.14047
Abstract
This release contains the two mathematica codes that reproduces all the key results in the manuscript "Exact, non-singular black holes from a phantom DBI Field as primordial dark matter". See https://arxiv.org/abs/2511.14047
Full text
Black Hole Solution of the DBI Lagrangian S = ∫ d4x-gR 2κ2+Λ41-1+2X Λ4 X is the Kinetic Energy of the DBI field Λ is the intrinsic energy scale Our metric signature is (-,+,+,+) We are going to assume the following metric solution ds2 = -f(r)dt2 + dr2 f(r) + ρ2 (r) dθ2 + sin2 θ dϕ2 ) where f(r)= 1+ 3GM a2 (rr2+a2 a ArcTan ( a r ) ) ρ2 (r) = ( r2 + a2 ) This Mathematica notebook uses “DBI_file1.nb” as its input
X is the Kinetic Energy of the DBI field Λ is the intrinsic energy scale Our metric signature is (-,+,+,+) We are going to assume the following metric solution ds2 = -f(r)dt2 + dr2 f(r) + ρ2 (r) dθ2 + sin2 θ dϕ2 ) where f(r)= 1+ 3GM a2 (rr2+a2 a ArcTan ( a r ) ) ρ2 (r) = ( r2 + a2 ) This Mathematica notebook uses “DBI_file1.nb” as its input In the following section, we calculate the Einstein tensor components and the curvature invariants of the spacetime metric In[ ]:= $Assumptions =G∈PositiveReals && M ∈PositiveReals && a∈PositiveReals && B ∈PositiveReals && κ∈PositiveReals; Defining the coordinate system: In[ ]:= n=4; coord = {r, θ,ϕ, t}; f[r]=1+3*(G*M) a^2 *r-(r^2+a^2) a*ArcTana r; ρ[r] = Sqrt[r^2+a^2]; Defining the metric and its inverse: In[ ]:= metric = {{f[r]^(-1), 0, 0, 0},{0, ρ[r]^2, 0, 0}, {0, 0, ρ[r]^2 *Sin[θ]^2, 0},{0, 0, 0, -f[r]}}; inversemetric =Simplify[Inverse[metric]]; Defining the Christoffel symbols: In[ ]:= affine :=affine =Simplify[ Table[(1/2)*Sum[(inversemetric[[i, s]])*(D[metric[[s, j]], coord[[k]]]+ D[metric[[s, k]], coord[[j]]]-D[metric[[j, k]], coord[[s]]]), {s, 1, n}],{i, 1, n},{j, 1, n},{k, 1, n}]] 2 DBI_file2.nb
In[ ]:= listaffine :=Table[If [UnsameQ[affine [[i, j, k]], 0], {ToString[Γ[i, j, k]], affine[[ i, j, k]]}],{i, 1, n},{j, 1, n},{k, 1, j}] TableForm[Partition[DeleteCases[Flatten[listaffine], Null], 2], TableSpacing → {2, 2}]; Defining the Riemann and Ricci: In[ ]:= riemann :=riemann =Simplify[Table[ D[affine[[i, j, l]], coord[[k]] ]-D[affine[[i, j, k]], coord[[l]] ]+ Sum[ affine[[s, j, l]]×affine[[i, k, s]]-affine[[s, j, k]]×affine[[i, l, s]], {s, 1, n}], {i, 1, n},{j, 1, n},{k, 1, n},{l, 1, n}] ] listriemann :=Table[If[UnsameQ[riemann[[i, j, k, l]], 0], {ToString[R[i, j, k, l]], riemann[[i, j, k, l]]}] , {i, 1, n},{j, 1, n},{k, 1, n},{l, 1, k -1}] TableForm[Partition[DeleteCases[Flatten[listriemann], Null], 2], TableSpacing → {2, 2}]; In[ ]:= ricci :=ricci = Simplify[Table[Sum[riemann[[i, j, i, l]],{i, 1, n}],{j, 1, n},{l, 1, n}] ] listricci :=Table[If[UnsameQ[ricci[[j, l]], 0], {ToString[R[j, l]], ricci[[j, l]]}] ,{j, 1, n},{l, 1, j}] TableForm[Partition[DeleteCases[Flatten[listricci], Null], 2], TableSpacing → {2, 2}]; Ricciscalar =Simplify[Sum[inversemetric[[i, j]]×ricci[[i, j]], {i, 1, n},{j, 1, n}] ] // FullSimplify Out[ ]= -2a5+21 a3G M r +18 a G M r3-9GMa2+r2 a2+2 r2ArcTana r a3a2+r22 In[ ]:= (*TeXForm[Ricciscalar]*) Defining the Einstein Tensor: In[ ]:= einstein :=einstein =Simplify[ricci -(1/2)Ricciscalar *metric] listeinstein :=Table[If[UnsameQ[einstein[[j, l]], 0], {ToString[G[j, l]], einstein[[j, l]]}],{j, 1, n},{l, 1, j}] TableForm[Partition[DeleteCases[Flatten[listeinstein], Null], 2], TableSpacing → {2, 2}] // FullSimplify; DBI_file2.nb 3
Defining the Curvature Invariants (Kretschmann scalar, Ricci tensor squared, Ricci scalar) In[ ]:= Ricciscalar =Simplify[ Sum[inversemetric[[i, j]]×ricci[[i, j]],{i, 1, n},{j, 1, n}] ] // FullSimplify Out[ ]= -2a5+21 a3G M r +18 a G M r3-9GMa2+r2 a2+2 r2ArcTana r a3a2+r22 ◼ We calculate the Ricci scalar at r=0 Ricciscalaratzero = Ricciscalar /. ArcTana r → π 2-ArcTanr a /.r→0// FullSimplify Out[ ]= -2 a +9GMπ a3 ◼ We calculate the form of the Ricci scalar at r → ∞ In[ ]:= Series[Ricciscalar, {r, Infinity, 6}] // FullSimplify Out[ ]= -2 a2 r4+36 a2G M 5 r5+4 a4 r6+O1 r7 ◼ We calculate Rμν Rμν In[ ]:= Riccitensorsquared = Simplify[Sum[inversemetric[[m, i]]×inversemetric[[l, j]]×ricci[[i, j]]× ricci[[m, l]],{i, 1, n},{j, 1, n},{l, 1, n},{m, 1, n}]] // FullSimplify Out[ ]= 1 a6a2+r24 4 a2a8+15 a6G M r +189 a2G2M2r4+81 G2M2r6+9 a4G M r2(13 G M +r)+9GMa2+r22 ArcTana r -aa4+12 a2GMr+18GMr3+3GMa4+3 a2r2+3 r4ArcTana r ◼ We calculate Rμν Rμν at r = 0 In[ ]:= Riccitensorsquaredatzero = Riccitensorsquared /. ArcTana r → π 2-ArcTanr a /.r→0// FullSimplify Out[ ]= 4 a2-18aGMπ+27 G2M2π2 a6 ◼ We calculate the form of Rμν Rμν at r → ∞ In[ ]:= Series[Riccitensorsquared, {r, Infinity, 10}] Out[ ]= 4 a4 r8-96 a4G M 5 r9-16 25 a6-39 a4G2M2 25 r10 +O1 r11 ◼ We calculate the Kretschmann scalar Rμνρσ Rμνρσ 4 DBI_file2.nb
◼ The Riemann tensor components evaluated above are of the form Rα μνσ ◼ Hence, we lower the first index of the Riemann Tensor to have the components in the form R αμνσ to calculate further In[ ]:= riemanndown =Table[Sum[metric[[l, m]]×riemann[[m, b, c, d]],{m, 1, 4}], {l, 1, 4},{b, 1, 4},{c, 1, 4},{d, 1, 4}] // Simplify; In[ ]:= K=Sum[inversemetric[[p, l]]×inversemetric[[q, b]]×inversemetric[[s, c]]× inversemetric[[t, d]]*riemanndown[[p, q, s, t]]×riemanndown[[l, b, c, d]], {p, 1, 4},{l, 1, 4},{s, 1, 4},{c, 1, 4},{t, 1, 4}, {d, 1, 4},{q, 1, 4},{b, 1, 4}] // Simplify Out[ ]= 1 a6a2+r24 12 a2a8+8 a6G M r +42 a2G2M2r4+18 G2M2r6+a4G M r2(33 G M +2 r)- 2aGMa2+r2 2 a6+30 a2G M r3+18 G M r5+a4r(15 G M +r) ArcTana r+ 9 G2M2a2+r22a4+2 a2r2+2 r4ArcTana r2 ◼ We calculate Kretschmannscalar at r=0 In[ ]:= Kretschmannatzero =K/. ArcTana r → π 2-ArcTanr a /.r→0// FullSimplify Out[ ]= 34 a2-8aGMπ+9 G2M2π2 a6 ◼ We calculate the form of Kretschmannscalar at r → ∞ In[ ]:= Series[K, {r, Infinity, 8}] Out[ ]= 48 G2M2 r6+32 a2G M r7+12 a4-16 a2G2M2 r8+O1 r9 Defining the DBI part : In[ ]:= ϕ[r] = 2*a B*(κ^2) π 2-ArcTana r* 3*(G*M)*r a^2 +1 Out[ ]= 2 a π 2-1+3GMr a2ArcTana r Bκ2 ◼ We calculate ϕ at r = 0 In[ ]:= ϕatzero = ϕ[r] /. ArcTana r → π 2-ArcTanr a /.r→0// FullSimplify Out[ ]= 0 ◼ We calculate the form of ϕ at r → ∞ DBI_file2.nb 5
In[ ]:= Series[ϕ[r],{r, Infinity, 3}] Out[ ]= -6 G M +aπ Bκ2-2 a2 Bκ2r+2 a2G M Bκ2r2+2 a4 3 B κ2r3+O1 r4 We consider the BH Thermodynamics in the (2GM>>a) limit ◼ We expand f(r) and consider terms up to second order in a In[ ]:= Series[f[r],{a, 0, 3}] Out[ ]= 1-2 G M r+2 G M a2 5 r3+O[a]4 ◼ We locate the horizon for the above spacetime metric In[ ]:= Solve[(1-(2GM)/r)+(2 G M a^2)/(5 r^3) ⩵ 0, r]// Simplify Out[ ]= r→1 3 2GM+4×51/3G2M2 -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3+ -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3 51/3, r→2 G M 3-2ⅈ51/3-ⅈ+ 3G2M2 3-27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3+ ⅈ ⅈ+ 3 -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3 6×51/3, r→2 G M 3+2ⅈ51/3ⅈ+ 3G2M2 3-27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3ⅈ -ⅈ+ 3 -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3 6×51/3 ◼ We expand the real solution up to second order in a In[ ]:= Horizon =Series1 32GM+4×51/3G2M2 -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3+ -27 a2G M +40 G3M3+3 81 a4G2M2-240 a2G4M41/3 51/3, {a, 0, 2} // FullSimplify Out[ ]= 2 G M -a2 10 (G M)+O[a]3 6 DBI_file2.nb
◼ We obtain the location of the horizon at r = 2GM - a2 /10GM ◼ Surface gravity is defined as 1 2 * f'(rh ) where rh is the horizon and f(r) is the negative of the gtt component of the metric ◼ We calculate the surface gravity up to second order in a In[ ]:= SurfaceGravity =1 2*D[(1-(2GM)/r)+(2 G M a^2)/(5 r^3), r] /. r-> 2 G M -a^2/(10 (G M)) // FullSimplify Out[ ]= 100 G3M3a4-100 a2G2M2+400 G4M4 a2-20 G2M24 ◼ We write the final expression for surface gravity which is correct up to second order in a In[ ]:= finalsurfacegravity =Series[SurfaceGravity, {a, 0, 2}] Out[ ]= 1 4 G M -a2 80 G3M3+O[a]3 ◼ Hawking Temperature is defined as TH = 1 2π * (Surface gravity) ◼ We write the expression of the corresponding Hawking Temperature correct up to second order in a In[ ]:= HawkingTemperature =1 2*π *finalsurfacegravity Out[ ]= 1 8 G M π-a2 160 G3M3π +O[a]3 ◼ We write the expression for the area of the horizon In[ ]:= Horizonarea =4*π*((Normal[Horizon])^2 +a^2) Out[ ]= 4 a2+ - a2 10 G M +2 G M 2 π ◼ We write the expression for the evaporation rate P correct up to second order in a In[ ]:= P= (HawkingTemperature^4)*(Horizonarea)// FullSimplify Out[ ]= 1 256 G2M2π3-a2 5120 G4M4π3+O[a]3 ◼ The first term in the above expression is the corresponding evaporation rate for the Schwarzchild metric and the second term is the correction. DBI_file2.nb 7