ABSTRACT (original article): Assuming a sufficiently general form for the matter distribution function of a galactic halo, we have derived solutions to the Einstein-Maxwell equations describing a charged black hole embedded in such a halo, while also allowing for a nonzero cosmological constant. These solutions generalize our earlier results for neutral black holes in asymptotically flat spacetime. As specific realizations of the general distribution, we consider the Hernquist, Navarro-Frenk-White, Burkert, Taylor-Silk, and Moore halo profiles, thereby capturing a broad range of astrophysically motivated scenarios. CITATION (original article): R.A. Konoplya, Z. Stuchlík, A. Zhidenko, Charged black hole surrounded by a galactic halo in a de Sitter universe, Phys. Rev. D 112, 083014 – 2025. https://doi.org/10.1103/l33g-6hlr

Main function

Here is the definition of FindNumericalSolution, which solves the Einstein equations using NDSolve for the charged spherically symmetric black hole surrounded by a galactic halo with given density profile. It finds numerically the mass function, the metric function B, and the electric potential V, which are matched to Reissner–Nordström-de Sitter asymptotic. The resulting piecewise functions are defined for any value of the radial coordinate.
FindNumericalSolution::usage="FindNumericalSolution[density profile,{coordinate variable,electric charge,event horizon,halo mass,halo radius,cosmological horizon}] solves Einstein equation for the spherically symmetric black hole surrounded by a galactic halo in a de Sitter universe and returns a solution as described in \!\(\*ButtonBox[\"R.A. Konoplya, Z. Stuchlík, A. Zhidenko, Phys. Rev. D 112 (2025) 8, 083014\",BaseStyle->\"Hyperlink\",ButtonData->{URL[\"https://inspirehep.net/record/2966363\"],None},ButtonNote->\"Click to open the INSPIRE page.\"]\).";​​FindNumericalSolution::largehalomass="Mass of the halo is too large.";​​FindNumericalSolution::wrongmass="Unable to solve numerically the equation for the mass function.";​​FindNumericalSolution::wrongequations="Unable to solve numerically the equations for the metric function.";​​FindNumericalSolution::nonvanishingdensity="Density function does not vanish at the event horizon `1`=`2`.";​​FindNumericalSolution[density_,{r_Symbol,Qv_?NumericQ,rv_/;rv>0,Mv_/;Mv≥0,sv_/;sv>0,cv_/;cv>0},opt:OptionsPattern[]]:=Withbasemassfunc=IfNumericQ[cv](*<Infinity*),
2
Qv
2rv
+
rv
2
-
2
Qv
2#
+
3
#
-
3
rv
6
3(
2
cv
rv-cv
2
Qv
-2cvMvrv+
2
Qv
rv-cv
2
rv
)
cvrv(
3
cv
-
3
rv
)
&,
2
Qv
2rv
+
rv
2
-
2
Qv
2#
&,Λv=IfNumericQ[cv](*<Infinity*),
3(
2
cv
rv-cv
2
Qv
-2cvMvrv+
2
Qv
rv-cv
2
rv
)
cvrv(
3
cv
-
3
rv
)
,0,​​If[Mv+basemassfunc[sv]≥sv/2,Message[FindNumericalSolution::wrongmass];Return[{}]];​​If[!PossibleZeroQ[density/.rrv],Message[FindNumericalSolution::nonvanishingdensity,ToString[r],ToString[rv,InputForm]]];​​With{FSol=NDSolve[{m'[r]density*r^2,m[rv]0},m,{r,rv,sv},opt]},​​If[!ListQ[FSol]||Length[FSol]≠1,Message[FindNumericalSolution::wrongmass];Return[{}]];​​With{FN=m/.First[FSol]},​​With{densitynorm=Mv/FN[sv]},​​WithBSol=QuietNDSolveB'[r]densitynorm*
r^2*B[r]*density
r-2*(basemassfunc[r]+densitynorm*FN[r])
,V'[r]-
Qv
r^2
*B[r],V[sv]
Qv
sv
,B[sv]1,{B,V},{r,rv,sv},opt,​​If[!ListQ[BSol]||Length[BSol]≠1,Message[FindNumericalSolution::wrongequations];Return[{}]];​​With{BN=B/.First[BSol],VN=V/.First[BSol]},​​With{rm=BN["Domain"]//Min},​​With{Brv=BN[rv]//Quiet,Brm=BN[rm],Vrv=VN[rv]//Quiet,Vrm=VN[rm]},​​QQv,r0rv,r1rm,MMv,ssv,rccv,ΛΛv,ρ[r]Piecewise
densitynorm
4π
*density,rv≤r≤sv,0,m[r]Piecewise[{{Mv,r>sv},{densitynorm*FN[r],rv≤r≤sv}},0]+basemassfunc[r],B[r]Piecewise
Brm(r-rv)+Brv(rm-r)
rm-rv
,rv<r<rm,{Brv,r≤rv},{1,r>sv},{BN[r],rm≤r≤sv},0,V[r]Piecewise
Vrm(r-rv)+Vrv(rm-r)
rm-rv
,rv<r<rm,Brv
Qv
r
-
Qv
rv
+Vrm,r≤rv,
Qv
r
,r>sv,{VN[r],rv≤r≤sv},0​​(*EndWithBrv,Brm,Vrv,Vrm*)​​(*EndWithrm*)​​(*EndWithBN,VN*)​​(*EndWithBSol*)​​(*EndWithdensitynorm*)​​(*EndWithFN*)​​/;Λv≥0​​/;rv<sv<cv;
Lapse funcion definition
f[r_]=1-
2m[r]
r
B[r]^2;

Example for the Dehnen-like profile, which vanishes at the horizon

Dehnen-like density profile depending on the radial coordinate r
densityprofile=
1+n
1-
horizon
r
α
a
r
α-γ
k
1
2
+
1
2
k
r
a
;
Solution for the fixed parameters (chorizon=∞ corresponds to the asymptotically flat spacetime):
BHSolution=BlockMHalo=
1
2
,sHalo=20,BHcharge=
1
10
,horizon=1,chorizon=∞,a=10,α=1,γ=4,k=1,n=0,​​FindNumericalSolution[densityprofile,{r,BHcharge,horizon,MHalo,sHalo,chorizon}]​​;
Lapse function plot
Plot[f[r]/.BHSolution,{r,1,50},AxesLabel{"r","f(r)"},PlotRangeAll]
For r>20 the solution coincides with Reissner-Nordström. Inside the horizon (r<1) the mass function is negative due to the electric charge.
Plot[m[r]/.BHSolution,{r,1/10,30},AxesLabel{"r","m(r)"},PlotRangeAll]
The density is nonzero for 1<r<20. The normalization coefficient is fixed to match the asymptotic mass.
Plot[ρ[r]/.BHSolution,{r,0,25},AxesLabel{"r","ρ(r)"},PlotRangeAll]

Calculation of the black-hole shadow radius

Einstein equations can be used to express derivatives in terms of the functions
FunctionDerivatives=
′
m
[r]
2
Q
+
4
r
Λ+8π
4
r
ρ[r]
2
2
r
,
′
B
[r]
4π
2
r
B[r]ρ[r]
r-2m[r]
;
We use r for the radial variable here
CalculateShadowRadius::usage="CalculateShadowRadius[solution] calculates the black-hole shadow radius for the solution obtained with FindNumericalSolution. It uses FindRoot in order to find the circular photon orbit position.";​​CalculateShadowRadius[L:{__Rule},opt:OptionsPattern[]]:=Withphotonorbitequation=D
f[r]
r^2
,r/.FunctionDerivatives/.L,​​Withphotonorbit=FindRootphotonorbitequation0,r,
3/2
1/r0+1/rc
/.L,opt,​​IfListQ[photonorbit]&&Length[photonorbit]1,
r
f[r]
/.L/.photonorbit​​;
We solve the Einstein equations with enhanced numerical precision. Then we calculate shadow radius for each solution. Since the black hole horizon r0=1, M is measured in units of r0.

Figures

Figure 3 from arXiv:2509.03301
Figure 4 from arXiv:2509.03301:

CITE THIS NOTEBOOK

Charged black hole surrounded by a galactic halo in a de Sitter universe​
by Roman A. Konoplya, Zdeněk Stuchlík & Alexander Zhidenko​
Wolfram Community, STAFF PICKS, December 16, 2025
​https://community.wolfram.com/groups/-/m/t/3592181