Performing a FindRoot from Numerical integration

Question

I am trying to evaluate:

FindRoot[Xfnew[γ, w, Hi, mx]*92*10^23*γ^(4/(1 + w)) == 1198/10000, {Hi, hs}, 
    MaxIterations -> 10000]

for some starting value hs, and γ and w are known numbers listed below in the definitions. I have to evaluate this FindRoot in a range of mx, and plot it, something like:

LogLogPlot[Hi /. FindRoot[Xfnew[γ, w, Hi, mx]*9.2*10^24*γ^(4/(1 + w)) == 
    0.1198, {Hi, 1.6946*10^-8}], {mx, 10^(-4), 10^(-10)}, PlotStyle -> Green]

I have successfully evaluated this in the range of mx (0.1 to 0.001), and (10^-10 to 10^-19), with starting values 0.0046 and 0.1 respectively.

But in the range 10^(-4) to 10^(-10), I am not able to. I tried using a "for loop" and break this range to form a list {mx,Hi} with small steps of mx(and doing a ListLogLogPlot afterwards), and perturbing the initial starting of FindRoot, like this:

For[i = 1, i < 2000, i++,h[0] = hs; h[i] = Hi /. Check[FindRoot[
    Xfnew[γ, w, Hi, mxs - i*step]*92.*10^23*γ^(4/(
    1 + w)) == 1198./10000, {Hi, h[i - 1]}, MaxIterations -> 10000], 
    {Hi -> 25}]; If [h[i] == 25, Break[], AppendTo[hin, {mxs - i*step, h[i]}]]]

where hs and mxs are some starting values of Hi and mx, and I break the loop if FindRoot generates error messages. And "step" is some stepsize I have chosen.

But, currently my step size for this loop has decreased so much, that it would take me ages to cover the entire range.I decrease my stepsize, whenever I get the error:

Machine precision is insufficient to achieve the accuracy

Is there a faster way to cover this range? I have tried increasing Working precision, it still doesnot help.

Function definitions

Xfnew is a result of a complicated list of functions listed below, and the values of gamma and w are listed below:: (Xfnew and the values of gamma and w are listed at the End) (I apologize for the Greek letters)

g = 61;
k2 = 25/100;    
k1 = 2/10;
ai = 1;
Trh[γ_?NumericQ, Hi_?NumericQ] := k2*γ*(( Hi)^(1/2))    
arh[γ_?NumericQ, w_?NumericQ] := ai*γ^(-4/(3 (1 + w)))
Hrh[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Hi (a/ai)^(-3 (1 + w)/2)
H[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Hi γ^2 (a/arh[γ, w])^-2
Trh11[a_?NumericQ, γ_?NumericQ, w_?NumericQ, 
    Hi_?NumericQ] := ((k1 (γ Hi)^(1/2))/((1 + 3/5 w)^((1/
    4)))) ((a^(-3 (1 - w)/2) - a^-4)^((1/4))) 
T[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ] := 
    Trh[γ, Hi] arh[γ, w]/a
nXeq[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ, T_Symbol] := g/(2 π^2)
    mx^2 T[a, γ, w, Hi] BesselK[2, mx/T[a, γ, w, Hi]]
σv0[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
    mx_?NumericQ, T_Symbol] := 3 π ((mx^2)/
    8) (3 BesselK[1, mx/T[a, γ, w, Hi]]^2/(5 BesselK[2, mx/T[a, γ, w, Hi]]^2) + 2/5 + 
    4/5 T[a, γ, w, Hi]/mx BesselK[1, mx/T[a, γ, w, Hi]]/
    BesselK[2, mx/T[a, γ, w, Hi]] + 8/5 T[a, γ, w, Hi]^2/mx^2) 
σv1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ,
    mx_?NumericQ, T_Symbol] := π (
    mx^2) (11/20 BesselK[1, mx/T[a, γ, w, Hi]]^2/
    BesselK[2, mx/T[a, γ, w, Hi]]^2 + 9/20 + 13/20 T[a, γ, w, Hi]/
    mx BesselK[1, mx/T[a, γ, w, Hi]]/BesselK[2, mx/T[a, γ, w, Hi]] + 
    13/10 T[a, γ, w, Hi]^2/mx^2)   
σvtotal[a_?NumericQ, γ_?NumericQ, w_?NumericQ, 
    Hi_?NumericQ, mx_?NumericQ, T_Symbol] := 
    4 σv0[a, γ, w, Hi, mx, T] + (45 + 12) σv1[a, γ, w, Hi, mx, T]
g1[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ] := (mx*a^2)/(Hrh[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
    a, γ, w, Hi, mx, Trh11] nXeq[a, γ, w, Hi, mx, Trh11]^2; 
g2[a_?NumericQ, γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, 
    mx_?NumericQ] := (mx*a^2)/(H[a, γ, w, Hi]*Trh[γ, Hi]^3) σvtotal[
    a, γ, w, Hi, mx, T] nXeq[a, γ, w, Hi, mx, T]^2; 
Xfnew[γ_?NumericQ, w_?NumericQ, Hi_?NumericQ, mx_?NumericQ] := 
    NIntegrate[g1[a, γ, w, Hi, mx], {a, 1, arh[γ, w]}] + 
    NIntegrate[g2[a, γ, w, Hi, mx], {a, arh[γ, w], ∞}]
γ = 1/10;
w = 5/10;

Answer 1

Experimentation indicates that a larger value of WorkingPrecision indeed is necessary. However, too large a value is computationally expensive and also sometimes causes the computation to crash. WorkingPrecision -> 30 seems to be a good compromise. In addition, computations for mx in the range of {10^-6, 10^-4} do not converge well without very good initial guesses. An initial guess as a function of mx fit to values of Hi for mx in the vicinity of 10^-2 and of 10^-7 and again a bit of experimentation suggests 0.2 mx^1.7 for mx > 10^-6. Then,

ParallelTable[{mx = N[10^i], 
    Hi /. FindRoot[Xfnew[γ, w, Hi, mx]*92*10^23*γ^(4/(1 + w)) == 1198/10000, 
    {Hi, Max[10^-11, .2 mx^(1.7)]}, WorkingPrecision -> 30, PrecisionGoal -> 6]}, 
    {i, -19, -1, .1}];
ListLogLogPlot[%, Joined -> True, PlotStyle -> Green, ImageSize -> Large, 
    AxesLabel -> Automatic, LabelStyle -> {Black, Bold, 14}]