Mismatch of agnpy's absorption with Finke (2016) reference: BLR

[cosimoNigro]

As suggested by @jsitarek in issue #50, I am opening a different issue for each of the absorption crosschecks.
This one regards the absorption on the photon field of the BLR.
I obtained the opacity vs energy computed at several distances by Finke - before I was using values I had fetched with webplotdigitizer from the paper's figures.




I have made a comparison for all the absorption on the single Lyman alpha line for the 4 distances represented in Figure 14 of Finke (2016): r=10^(-1) R(Ly alpha), r=R(Ly alpha), r=10^(0.5) R(Ly alpha), r=10 R(Ly alpha).
The agreement is good only at small distances (within the BLR) and deteriorates at larger ones (outside the BLR).

[cosimoNigro]

Copying @jsitarek 's last comment from issue # 50.

Some more updates on the BLR absorption.
Somehow the energy threshold that I'm getting does not agree with what you get.
The line in the example is 10.2 eV, the threshold E_thr = 2 m_e^2 / (epsilon * (1-cos theta)). For head-head collisions you get cos theta=-1 and you end up with E_thr = 2.56e10 eV, i.e. nu = 6.2 e24 Hz. In the frame of the observer you would get a factor of 1+z lower value i.e. nu = 3.35e24 Hz, while Finke's and agnpy show 6.6e24 Hz.

I did also full calculations and assymptotically at high energies I get the values consistent with agnpy, but at energies close to the threshold, consistently with what I wrote above I get higher values. See the attached notebook, in particular [10]
check_tau_blr_js.ipynb.gz

[jsitarek]

The plots that you did clearly show that there is something wrong with agnpy code, when you put the emission region at the distance of the RLy alpha the angles change considerably compared to the 0.1 R Ly alpha case, but the threshold still stays the same for agnpy.
I suspect that there can be something wrong with the integration over the angles, especially when the distance gets larger the actual integration shrinks to a very small angular size with strong dependence and this is numerically tricky

[jsitarek]

I think I can at least partially explain the difference that you got in those tests. The formulae 123 and 126 in Finke's paper use mu_star in the place where you use cos psi, i.e. he assumes that the source is moving along the jet axis. This is becoming very important if the emission region is located far from the BLR, please see an example

Can you rerun the plots assuming that theta_s is ~0 (i.e. Delta =~2 * Gamma)

There is another thing that worries me if you want to make it really 3D with absorption being angle dependent.
Look at Fig 9 in Finke's paper, The vertical line (along which r is measured) might seem to go along the jet axis, but this is the case only for the moment of the emission of the gamma ray. When it is moving towards the observer mu_s will be getting smaller that will also affect the cos_phi angle (looking at cos_psi function I do not think you take this into account).

[cosimoNigro]

Thank you @jsitarek, replying to the first part of your message...
I am checking it now, the mistake on my side was introducing the dependency on the entire blob when initialising the absorption class Absorption(blob, target, r=r). That way I never know or check the actual value of the parameter mu_s, the angle the blob forms with the line of sight, but I always fetch it from the blob object, i.e. I do self.blob.mu_s somewhere in the code.

I think it would be good to make it explicit in the absorption definition and to remove the dependency on the entire blob. In the end the only two properties of the blob we use are blob.mu_s and blob.z.
So I propose to change the class initialiser to
Absorption(blob, target, r=r) to Absorption(target, r, mu_s, z)
Where we can set by default mu_s=0 and z=0, for example (that would be absorption for a blob aligned with the line of sight and in the galaxy reference frame). What do you think?

[jsitarek]

yes, I think it makes sense, the absorption is not connected with the emission region (unless we add at some point absorption in synchrotron radiation). so specyfying just the starting r, mu_s and z is the best solution.

[cosimoNigro]

After removing the dependency on the Blob from the Absorption (PR #76) things have improved a little, not too much though, there is still some mismatch when going outside the BLR. In the agnpy comparison now mu_s=1, before mu_s was taken from the blob.



(for larger distances the absorption values provided by Finke are null at all energies).
Again, the fact that outside the BLR there is some disagreement, but in the most extreme case we can successfully match the point-source approximation is really strange.

Might there be a problem with the paper's calculations?
We use the same cross section and the same integration method for all the absorption cases. So it's a bit difficult to figure out what's going wrong.
Also it is strange that alredy for 1e2 * R_Ly_alpha = 1e19 cm Finke gets already null values at all energies (check the sampled_taus in agnpy/data) while we can go to 1e20 cm with no problems and match the point-source case.

[jsitarek]

please note that the comparisons with point source are very tricky, because for large distances if you assume that the observer is along the jet axis, for a point like source the angle between photons is exactly 0, which means that tau=0.
It makes more sense to compare BLR vs DT, for those large distances, there will be still some difference in angles that would affect the calculated taus, but a large effect like this seen in those plots we should be able to still see

[jsitarek]

Hi @cosimoNigro

I did some tests comparing the DT and BLR absorption (taking T_DT = epsilon_lin * mc^2 / k_B)

1. if distance/radius = 0.1 and mu_s =1 tau_BLR = tau_DT *pi/2 - this might come from different geometry
   EDIT: thinking more about this, it does not look right. BLR, has also part in front of the photon, so the peak should be shifted by something close to a factor of 2 to the left

2. if distance/radius = 1.1 and mu_s =1 tau_BLR is a bit higher then tau_DT and goes down to lower energies - this is reasonable because the tangential to the BLR gives you better absorption angle then the angle to DT, and you still have a large fraction of BLR with such good angles

3. if distance/radius = 10.1 and mu_s =1 tau_BLR is less then tau_DT and also shifted a bit to higher energies - also understandable, this time all the angles are very unfavourable for the absorption, and most of the BLR is more close to the jet axis then DT

4. if distance/radius = 10.1 and mu_s =0.9 both absorptions give the same, in this case both DT and BLR are effectively a point-like source. I find it only somewhat strange that it is enought to change mu_s so little, I would expect that at acos(mu_s) ~ asin (radius/distance) you should still have some effect from the angles. I will check this part of the code carefully

[jsitarek]

I've rerun the dist / radius = 0.1 case with the code of PR #79

and now it looks much more reasonable, because of the part of the shell before the blob the absorption is stronger and also the threshold energy more or less agrees with the theoretical one (6.2e24 Hz for BLR, 1.2e25 Hz for DT)

[jsitarek]

Hi @cosimoNigro , after this PR #79 change I also checked the actual comparison with Finke's points.

I plot them as they would be calculated with z=0, since this is how it looks like they are presented in the paper (I will write about this also in the DT issue, because there is something really fishy there).

The agreement is very nice, for the case of 0.1 R_line it is definitely nicer than before the fix (the fix did not affect the case of r>R_line).
For higher values there are some differences but they are not large, can be related with the accuracy of the computations.

I think we can close this issue, but probably it would be worth to open another issue with a comparison of the BLR over all the lines, since @pawel21 is working on implementing this and the curves for all the lines are also in the Finke's paper (however I do not see them in the reference_tau directory)

[cosimoNigro]

@jsitarek can we close after #89?

[jsitarek]

closing