Max observed flux vs Luminet's paper

[cesmoak]

In Luminet's 1979 paper, Fig 10 notes that the max observed flux for the 10 degree image is 2.62. Running your code, I get max observed flux of 2.81. I have made a separate implementation of Luminet (ray tracing forward from the viewer) that gets 2.18. I can't seem to figure out how to things to line up with 2.62.

Do you know if 2.62 number is incorrect due to some of the typos in the paper, or are you aware of what changes are necessary to match this number?

Cheers!

[bgmeulem]

It's very odd, but I can't seem to reproduce either one of these numbers.
Luminet for sure has some typo's, I can't be sure if I weeded all of them out. That being said, I have also made plenty errors reproducing his results before. I can't be sure I weeded all of those out either.

I'll try to shed some light on these values below, hopefully we can work it out together.

The radius of max intrinsic flux is said to be about 9.55, but if we want a more accurate value, we can numerically solve Equation 15:

# Define optimization function for M=1 and acc_rate = 1. These have no impact on the result for r
def opt_func(r): return - bhmath.calc_flux_intrinsic_swarzschild(1., r, 1.)
from scipy.optimize import minimize as opt_minimize
r_max_flux = opt_minimize(opt_func, 10., tol=1e-12).x[0]
f_s_max = bhmath.calc_flux_intrinsic_swarzschild(1., r_max_flux, 1.)

This gives us more precise numbers:

• radius for the max flux: r_max_flux = 9.550928098667045
• a maximum intrinsic flux of f_s_max = 1.3678110819829872e-05

Plotting the intrinsic flux along an isoradial of radius r_max_flux in units of F_smax shows:

ir = Isoradial(radius=r_max_flux, incl=np.pi*(80)/180, bh_mass=1., angular_resolution=500)
f_norm = calc_flux_observed(r=r_max_flux, acc=1, bh_mass=1, redshift_factor=ir.calc_redshift_factors()) / f_s_max

plt.plot(f_norm)
plt.scatter(np.argmax(f_norm), np.max(f_norm), label="$F_{smax}=$"+f"{max(f_norm)}", c='red')
plt.axvline(np.argmax(f_norm), c='k', label="$\\alpha=$"+f"{np.argmax(f_norm)/500}")
plt.xlabel("$\\alpha$ (radians)")
plt.ylabel("$F_o$")
plt.legend()
plt.show()

So the maximum observed flux in units of max intrinsic flux should be about $F_o/F_{smax}=2.768$. Increasing angular precision to 50k bumps up the flux estimate to about $2.76930$. This is probably close to the numerical solution I didn't bother to implement.

Running the code with the normal approximate values (which should be the implementation on master):

• $R \approx 9.55$
• $F_{smax} \approx 1.146e-4 * 3/8\pi$

gives me $F=2.7170399411715374$.
Yet another number.
All ballpark close, but all too far away to be just roundoff error.

In any case, it should be noted that the maximum flux is very steeply peaked around this value. Any small deviation in angle or radius will quickly lead to a significant underestimation of the max flux. Slight deviations for $F_{smax}$ will normalize them slightly wrong as well.

Could you show me how you got to 2.81 with my code?
Could you tell me which values you get with your code if you calculate using the more exact numerical solutions?
Regarding Luminet's 2.62, I have no clue sadly... Maybe it's a typo. Maybe calculating this value was more of an afterthought and they didn't spend too much numerical precision on it. Maybe I made a mistake somewhere in my implementation.