Thick disc model

[phajy]

It would be good to produce a Gradus-based thick disc model that can be used to fit spectra.

The sorts of parameters that could be varied include:

• metric: multiple deformation parameters (or an arbitrary metric)
• spin
• inclination angle
• emissivity index or self-consistent emissivity from a lamp-post model
• inner and outer radii of the disc
• Eddington fraction of the disc

To start with it would be great to quantitatively reproduce the results of Abdikamalov et al. (2020). Let's start with their Figure 5.

Perhaps the easiest way to do a quantitative comparison is to compare with the output of RELLINE_NK.

[phajy]

See astro-group-bristol/SpectralFittingExtras.jl#7 regarding relline_nk.

[wiktoriatarnopolska]

Reproducing results from Abdikamalov et al. 2020

Now that the photon index q=3 was specified as in the Abdikamalov paper, the difference is quite striking. Code can be found here

[fjebaker]

The code that calculates emissivity profiles for different coronal spectra has been tested against that Gonzalez paper, right? I can't find any unit tests in Gradus.jl for it, so just wanted to be sure it's verified somewhere.

[wiktoriatarnopolska]

yess, we did that last summer. we tested the moving source and our results didn't agree with Gonzalez+2017 nor Wilkins+Fabian2012 BUT we did agree with Dauser+2013

[fjebaker]

I remember that, but for different corona powerlaw exponents? As in, we're definitely sure that bit of the code is working as intended?

[wiktoriatarnopolska]

yes, we did a smoke test to check if PowerLaw(2) is giving the same results as the default here and also reproduced Fanton figures with different emissivity indexes here

[wiktoriatarnopolska]

[fjebaker]

The Fanton paper has that bluring interpolation though, and also is normalized differently. Just looking at those plots side by side isn't convincing to me; are there not examples in Gonzalez for different power laws at better resolution that we can compare to?

[fjebaker]

This one.

[wiktoriatarnopolska]

uu I don't think we did that specifically but that can be done, I'll run the test

[wiktoriatarnopolska]

Recreated the Gonzalez figure:

also ran the test again to compare the PowerLawSpectrum() function with the default to verify if we actually use photon index = 2

[phajy]

Reproducing results from Abdikamalov et al. 2020

Now that the photon index q=3 was specified as in the Abdikamalov paper, the difference is quite striking. Code can be found here

I think the issue here might be that the emissivity index $q = 3$ rather than the photon index. They don't self-consistently calculate the illumination patter on the disc but yes an emissivity law $I_e \propto \rho^{-3}$ where $\rho$ is radial coordinate.

[wiktoriatarnopolska]

Reproducing results from Abdikamalov et al. 2020
Now that the photon index q=3 was specified as in the Abdikamalov paper, the difference is quite striking. Code can be found here

I think the issue here might be that the emissivity index q=3 rather than the photon index. They don't self-consistently calculate the illumination patter on the disc but yes an emissivity law Ie∝ρ−3 where ρ is radial coordinate.

PowerLawSpectrum defined as here: https://github.com/astro-group-bristol/Gradus.jl/blob/4bc304ab7a571334cfeded8b8c3c96ef089dbcde/src/corona/spectra.jl#L2

[wiktoriatarnopolska]

Following up on our meeting today, we actually do have the non-self-consistent version of the Abdikamalov recipe and the results from here do agree with Abdikamalov+ (2020) paper because they use the emissivity power law ~ power of -3 and so do we by default as per here. The reason our self-consistent figures did not agree with his results is because I mistakenly assumed emissivity and photon index are synonymous (mea culpa), however that prompted me to perform the same calculations using the self-consistent version of the code this time, replacing the PowerLaw argument with the default value of 2.
Running that model yielded different resulting line shapes.
Reproduced Abdikamalov fig. 5 (inclination angle of 20 degrees) calculated using the broken power law agrees with their results:

But the self-consistent version implementing the following
model = LampPostModel(h = 10.0)
spectrum = PowerLawSpectrum(2.0)
does not:

But this is expected since the self-consistent code calculates the illumination pattern of the disc returning the profile, while the broken power law uses the power of -3.
Abdikamalov mentions in his paper that his group used the broken power law profile, whilst Taylor & Reynolds (2018a) used the lamppost geometry to calculate theirs.

I think the self-consistent method might be worth exploring and we could try recreating Taylor & Reynolds figures to verify our approach first and then expand their work -- they only implemented their models for inclination angles 15 30 and 60 degrees and did not include the deformation parameter. Abdikamalov expanded their work, adding higher inclination angles and the deformation parameter but compromised the self-consistency of their approach. With Gradus we can do all of it -- cover the high inclination angles, explore the impact of the deformation parameter in the Johannsen metric and provide the self-consistency in the method using the bit of code I wrote over the summer that calculates spectra and the profile passing that as argument.

If Abdikamalov was expanding Taylor&Reynolds work by employing Johannsen's metric with the deformation parameters and higher inclination angles but used the power law instead of self-consistently calculated profiles, compromising this consistency, we could do all of those things -- therefore making a difference, providing a new approach and improving the models, right?

Let me know your thoughts if this idea is worth exploring!
Summoning everyone! @phajy @fjebaker @HannahLidgett

[phajy]

Yes! I think this is a great idea. We should do the self-consistent calculation with deformation parameters, going beyond the Abdikamalov et al. (2020) and Taylor & Reynolds (2018) papers.

[wiktoriatarnopolska]

Test Taylor&Reynolds:

• reproduce figs. 6, 7, 8 (in relation to Abdikamalov)
• reproduce fig. 3 (essentially the emissivity profile) <-- see Thick disc model #3 (comment)
• fig. 4
• fig. 5

Effects of the self-consistent illumination and metric dependency:

• figures varying height for chosen Abdikamalov's parameters (deformation parameter never done self-consistently)

[wiktoriatarnopolska]

Taylor&Reynolds tests -- results

Taylor&Reynolds2018a recipe for figures 6-8

note:
Differs for $a = 0.0$ at $h = 12$, specifically at edd. ratio = 0.
Differs for $a = 0.9$ at $h=3$ across inclinations, particularly end. ratio = 30%
Differs for $a = 0.99$ at $h=3$, most significantly at $inclination 15 and 30$

[wiktoriatarnopolska]

• reproduce fig. 3 (essentially the emissivity profile) <-- see [Thick disc model Thick disc model #3 (comment)]

so I am actually struggling with it a little bit. I am confused at how they are defining their inner and outer radii. so for now I'm simply not defining those.
the plots so far kind of agree but I am a bit confused. Perhaps defining inner and outer radii would help but if we could discuss this in the next meeting on Friday that would be great!

[wiktoriatarnopolska]

so far the code is here and this is the result:

[fjebaker]

Squeeze your disc together a bit: no point plotting out to 1000 $r_{g}$ for most of them, so maybe set limits between $0$ and $100$ $r_{g}$. The flat disc also seems to be extending within the ISCO!

[wiktoriatarnopolska]

fixed the flat disc extending into the ISCO by defining thin disc as dthin = ThinDisc(Gradus.isco(m), Inf) (yay!)
I want to also add the same line at $r = r_{ISCO}$ as Taylor & Reynolds did but I struggle -- I tried it by passing the metric as argument into plot_all function but I'm getting an error cannot convert Float64 to series data for plotting? I will try to work it out in a bit but so far so good:

Finding the best limits for r now.

[wiktoriatarnopolska]

Better limits for r and fixed the order so it matches T&R 😄 code file shall be updated soon

[fjebaker]

Looks good! Maybe the y axes need adjusting? Top row especially, there is a lot of empty space in the plots!

[wiktoriatarnopolska]

working on the fig. 4 code recipe -- map of the discs' redshift. results:

@fjebaker is there a condition for stopping integration of the photons inside the ISCO since we see the redshift from the photon ring? also, any way to give an outer edge to >0% Eddington ratios in Shakura--Sunyaev models? passing minrₑ and maxrₑ makes rendergeodesics unhappy (unrecognised keyword arguments found).

also I feel like at low Eddington ratios especially our BH shadow isn't the same shape..? might be due to heights/widths but that's just my guess..?

code available here :)

[wiktoriatarnopolska]

@fjebaker is there a condition for stopping integration of the photons inside the ISCO since we see the redshift from the photon ring? also, any way to give an outer edge to >0% Eddington ratios in Shakura--Sunyaev models? passing minrₑ and maxrₑ makes rendergeodesics unhappy (unrecognised keyword arguments found).

mainly because since our line profiles figures differ slightly from the source ones we suspect the photons from photon ring might travel back and hit the disc therefore affecting the shapes. unlikely, but would be neat to check just to be 100% sure

[fjebaker]

Yeah you can filter the geodesics. For the inner horizon, I already provide callback = domain_upper_hemisphere(), and for setting an outer radius, you can use a FilterPointFunction which filters the radii to within a certain r

r_filter = FilterPointFunction((m, gp, t) -> gp.r[2] < 100, NaN)
pf = ConstPointFunctions.redshift(m, x) ∘ r_filter ∘ ConstPointFunctions.filter_intersected()

[phajy]

I've written a routine that allows us to calibrate images so we can overlay Gradus (or any other) plot on top. I've pushed an example to the repository at overplot_example.jl which produces the following plot.

This is a comparison with Taylor & Reynolds (2018) Figure 7 bottom left panel ($a = 0.9$). Looks pretty good. The normalisation was adjusted by hand. There is still some minor disagreement with the blue curve, but that might not be significant.

The image calibration was done with Overlay.jl which is functional but perhaps not as polished as it could be.

[fjebaker]

Wow that is better agreement than I would first have thought! Well done all and very cool package @phajy!
Now to work out what their normalization is...

[wiktoriatarnopolska]

@phajy that looks amazing!! really nice agreement too!

[phajy]

Yes, it is good agreement! The normalisation I tuned by eye having set the normalisation of the $\dot M = 0$ to have a maximum of 1, as per the Taylor & Reynolds (2018) paper.

SelfConsistentIronLines/src/overplot_example.jl

Lines 90 to 93 in d41b09b

	m[1] = maximum(data.f[1])
	m[2] = 1.08*m[1]
	m[3] = 1.2*m[1]
	m[4] = 1.4*m[1]

I would prefer a different normalisation because the height of the sharp peak probably depends on the binning and might be subject to aliasing. Perhaps a normalisation based on integrated flux.

[wiktoriatarnopolska]

that's an eye for detail! I will have a deep dive into the Overlay.jl package, looks super interesting -- we can really neatly fit the data now. We can definitely improve the normalisation, but it's a good sanity check that our code does what it's supposed to and the results agree 🎉

[wiktoriatarnopolska]

Fixed up the filters on redshift maps!

[wiktoriatarnopolska]

his all I hope your Easter break is going great<3

I've been looking at and verifying results I revisited the reproduced figure from Fanton et al 1997 and I noticed the wrong parameter was varied -- we meant to vary the emissivity index and instead the photon index was varied.

I fixed the code but now our results differ majorly -- please see below.

I posted the code here. Please confirm if this looks correct now :( I can see they are calculating emissivity slightly differently -- we do ε(r) = r^(-3) and Fanton includes an extra factor of epsilon_0/4pi which we do not do. but could it impact the profile so much?? I beg for a sanity check.

summon @phajy @fjebaker

[wiktoriatarnopolska]

but then now that I think about it -- because of their normalisation they look very differently but actually, in fact I think the results now match better

[fjebaker]

plot(gs, flux, legend=true, label = "p = 2.0", colour =:hotpink2, title = "a = 0.0, r_in = 6.0, r_out = 20.0, θ = 30",  xlabel = L"\textrm{Observed~frequency~shift~} (\nu_o / \nu_e)", ylabel = L"\testrm{Flux~(arbitrary~units)")

line 34 in your code, you are plotting flux not flux1. If you correct it

[fjebaker]

You likely had flux still initialized from another parameter combination

[fjebaker]

Good on you for checking though! And well noticed regarding the parameters 👍