continuum absorption accuracy

[ajwheeler]

There's lots to talk about in the figure comparing the continuum from korg, moog, and TS, but I think the jump in the Korg continuum at 3846\AA (the H2+ bf+ff cutoff, marked with a grey dashed line) is indicating that this source of absorption is overestimated.

Thoughts?

[mabruzzo]

So, the H₂⁺ formulation used by MARCS (and Turbospectrum) looks really easy to implement. I'll take a stab at that and see if it makes any difference. (The MARCS paper seems to suggest that this H₂⁺ source only includes ff absorption, and they don't list a separate source for bound-free absorption..., so maybe that's the answer?)

• Relatedly, it the future it may make sense to implement the coefficients from here. This is a more recent source than the source used for either our data or for Turbospectrum.

[ajwheeler]

Please don't do this unless you really, really want to. I'm happy to take it on myself, and I'm not sure that we should do something just because other code do. Is what we are doing a better or worse approximation?

[mabruzzo]

I don't know if what Turbospectrum is doing is better or worse (it's probably not better). It might be a different fit to the same data....

I didn't see you response until now... This is how far I got (feel free to pick up where I left off):

using Korg
using Interpolations: LinearInterpolation

const _C_interpolator, _D_interpolator = begin
    const table = [0.05    1.821 -0.01917;
                   0.10    2.556 -0.04473;
                   0.15    2.981 -0.07137;
                   0.20    3.275 -0.09884;
                   0.25    3.499 -0.1268;
                   0.30    3.684 -0.1544;
                   0.35    3.836 -0.1845;
                   0.40    3.967 -0.2138;
                   0.50    4.182 -0.2738;
                   0.60    4.356 -0.3354;
                   0.70    4.499 -0.3985;
                   0.80    4.611 -0.4626;
                   0.90    4.727 -0.5280;
                   1.0     4.817 -0.5946;
                   1.11111 4.915 -0.671;
                   1.25    5.013 -0.766;
                   1.42857 5.127 -0.893;
                   1.66667 5.243 -1.067;
                   2.0     5.395 -1.319;
                   2.5     5.551 -1.716;
                   3.33333 5.751 -2.424;
                   4.0     5.863 -3.020;
                   5.0     5.989 -3.976;
                   6.66666 6.099 -5.708;
                   10.0    6.289 -9.606;
                   20.0    6.235 -26.70]

    νs = (Korg.c_cgs.*(1.0e4 .* table[:,1])) # Hz
    println(Korg.c_cgs * 1e8 ./ νs)
    LinearInterpolation(νs, table[:,2]), LinearInterpolation(νs, table[:,3])
end
#println(size(table))

function _H2plus_mihalas(ν, T, nH_I, nH_II)
    # implicitly includes term for stimulated emission
    mass_H = Korg.atomic_masses[1]

    exponent = _C_interpolator(ν) + _D_interpolator(ν)*(5040.0/T)

    mass_H * 10.0^exponent * nH_I * nH_II / Korg.blackbody(T, Korg.c_cgs ./ ν)
end

I haven't actually tested it and it's unclear to me how many factors of mass_H are needed (it could be 0, 1, or 2). This came from Table A.I.2 of this paper.

[ajwheeler]

This goes away if you switch off H2+ bf+ff, confirming that this is, indeed the problem.

[ajwheeler]

I'm continuing to look into this. I'll definitely look at the MOOG/TS version, but I think a good way to go would be to update our prescription from Bates 1952 to Stancil 1994, as you suggested. This will also give us coverage down to 500 Å, although the temp lower bound will go up (perhaps, we could interpolate?). It's also 40 years more recent, and presumably much more accurate.

Also, I didn't realize that H2+ bf refers to H2 II + gamma -> H I + H II ! That should definitely be documented! (Sorry if I missed it.)

[ajwheeler]

To be clear, is TS using the data form the Mihalas paper? I can't actually find the part of the code where they do this...

[ajwheeler]

OK, so the new (Stancil) implementation reproduces the continuum offset with other codes, as expected, since it matches the Bates/Gray values. The good news is that we don't have the problematic discontinuity around 3850 \AA, where the old values stopped.

Manually turning it off gets a closer match, but I'm now pretty convinced we aren't actually doing anything wrong RE H2plus. It may be that something else is overestimated, or that the other codes are underestimating somthing.

[ajwheeler]

I've checked that the ff contribution from H2plus is very small and does not account for the discrepancy.

[ajwheeler]

Comparisons with Turbospectrum in other stars have suggested that at least one source of continuum opacity is significantly overestimated in Korg. (An alternative explanation is that there is something wrong with my radiative transfer. I'm not ruling this out. I suppose it could also be the molecular equilibrium but I wouldn't bet on it.)

Here's an underestimated (wrt TS) continuum for two giants (edit: ignore there legend here. This is using all continuum opacity sources in Korg):

This persists if H- bf is the only source turned on for Korg. This suggests that there is a problem with that particular source. TS is using the same table of cross-sections (Wishart+ 1979), so that's not the problem. There could be a disparity arising from how n(H-) is calculated. TS uses Mihalas 1967, which I haven't found a PDF for...

[mabruzzo]

This looks crazy!

[ajwheeler]

OK, so turbospec adjusts the ionization energy of H- for pressure, according to this paper: https://articles.adsabs.harvard.edu/pdf/1969tons.conf..143T. I'm not sure if this is a large enough effect to explain the difference, but I'm looking into it. It's probably something we should be aware of anyway.

[ajwheeler]

An overview of this effect is found on Hubeny and Mihalas p244-246, which states that ionization potentials can be adjusted by 3e-8 * Z * ne^0.5 T^-0.5 eV as a rough approximation.

julia> [3e-8  * sqrt(l.electron_number_density / l.temp) for l in atm.layers]
 3.338923987225036e-6
 3.98622586156897e-6
 4.746751058123049e-6
 5.628540558561197e-6
 6.649828953533077e-6
 7.826763746735217e-6
 9.165571895018155e-6
 1.0669809772082763e-5
 1.2337084775531654e-5
 1.4177154541496423e-5
 1.6225400209165468e-5
 1.7291895559478315e-5
 1.843716023129176e-5
 ⋮
 0.0009883337060364432
 0.0011845003364186083
 0.001384572744372603
 0.0015874316679815375
 0.0017938935198369107
 0.002004763014170044
 0.002221103664318355
 0.0026648894546756753
 0.0031523364249812463
 0.0036892123015088555
 0.004280991706155335
 0.00493208728816118

Too small to be the cause of this difference. (The ionization energy of H- is about 0.75 eV, for reference.)

[ajwheeler]

I've tracked down some problems with how we were interpolated model atmospheres, and am am going to switch Korg over to use different atmosphere params internally (see #97). This accounts for some discrepancy between Korg and Moog/TS. I still don't understand entirely what's going on.

In the meantime, here's a plot comparing alpha extracted from Moog and Korg for the sun at each layer in the atmosphere.

[ajwheeler]

With, #97, the continua look like this:

Note that our continuum absorption is still larger than Moog's, but it matters less with the renormalization to the model atmosphere tau, and it's not clear to me that Korg is less correct. I'll leave this issue open, since this remains an interesting thing to investigate, but it's no longer critical.

[mabruzzo]

So I have a naive question, why is there any disagreement at all at 5000 Angstroms? Aren't all of the code effectively using the same opacity at that wavelength?

Or do they all agree exactly when you consider the full spectra (when you additionally include the opacity from lines)?

[ajwheeler]

That is the naive expectation, yes. I think the differences come down to slight differences in how radiative transfer is done, but I could be wrong.

On Thu, Jan 6, 2022 at 4:52 PM Matthew Abruzzo @.***> wrote:

So I have a naive question, why is there any disagreement at all at 5000 Angstroms? Aren't all of the code effectively using the same opacity at that wavelength?

Or do they all agree exactly when you consider the full spectra (when you additionally include the line list).

— Reply to this email directly, view it on GitHub https://github.com/ajwheeler/Korg.jl/issues/83#issuecomment-1006960855, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAFN2GY367AIOXYE647JZQ3UUYFIPANCNFSM5GRJ6TQQ . You are receiving this because you authored the thread.Message ID: @.***>

[ajwheeler]

More confirmation that H2+ is really the culprit.

With H2+ opacity turned on:

With it turned off: (The missing opacity in the infrared is presumably metals.)

As far as I can tell, H2+ isn't in MOOG. I still can't see anything wrong with our implementation, so it really does seem that there's nothing wrong here. I'll return to this while/after updating radiative transfer and adding metal opacities.

[ajwheeler]

I got around to comparing our opacities to those form MARCS (which are the same as those from turbospectrum). Here's the same plot, indicating that we are doing pretty well. same thing with compressed colorbar:

[ajwheeler]

A view extending down to lower wavelengths:

[ajwheeler]

By the way, this is what bf metal absorption looks like. It's not at the right wavelength and tau (nor is it strong enough) to account for the discrepancy.

edit: colorbar units are cm^-1

[mabruzzo]

That is fascinating - although the y-axis is a little confusing. Could you remind me what that means again?

[ajwheeler]

It's log(optical depth), so the top (bottom) of the plot is the top (bottom) of the atmosphere.

[ajwheeler]

Looking at the MOOG output, it turns out that the difference is due to metal opacity, specifically Al I. Currently trying to figure out why we don't reproduce it.

[mabruzzo]

Fascinating. I have a faint recollection of seeing Al I opacity in the MOOG source code, and seem to recall that it was:

• overly simplistic
• somewhat unclear what publication it came from

[ajwheeler]

I believe all the sources in MOOG are form ATLAS.

From the ATLAS report:

MOOG implementation, which I think it from the ATLAS source.

      subroutine opacAl1
c******************************************************************************
c     This routine computes the bound-free absorption due to Al I.
c******************************************************************************

      implicit real*8 (a-h,o-z)
      include 'Atmos.com'
      include 'Kappa.com'

      do i=1,ntau
         if (freq .ge. 1.443d15) then
            aAl1(i) = 6.5d-17*(1.443d15/freq)**5*6.*
     .                numdens(5,1,i)/u(13,1,i)
         endif
      enddo

      return
      end

edit: this code doesn't quite match the screenshot. It has a steeper nu dependence and a larger constant, i.e. a "pointer" edge profile. I imagine this can from a subsequent version of ATLAS.

[ajwheeler]

Ok, so I've done some work on the metal bf opacities. Currently I'm only including just using neutral C, Mg, Al, and Si. These are what's in MOOG (minus Fe) and they are also the species that MARCS does an empirical energy level adjustment for. I'm also doing the adjustment, though I'm doing all energy levels, while MARCS only did those in the UV (?).

With that change, the delta from MARCS opacs looks like this: Our continuum looks like this:

There are two possible things going on with this.

• Our metal bf opacities are very spiky, possibly too much so. I did convolve them with a small doppler profile, to smear the edges a bit, but perhaps MARCS smears them more? We have more "features" in the UV than Turbospectrum does, and our alpha is overestimated relative to MARCS at the energy "edges", suggesting this is true. Another possibility is that I have some sort of bug which "compresses" the cross-section for each electron configuration in wavelength space, making it spikier (???)

and/or

• I just need to add more metals. More blanketing would make the UV less "spiky" maybe, and since we are using NORAD for some species, it's not surprising if our metal opacities are higher than MARCS's in some places.

[mabruzzo]

Our metal bf opacities are very spiky, possibly too much so. I did convolve them with a small doppler profile, to smear the edges a bit, but perhaps MARCS smears them more? We have more "features" in the UV than Turbospectrum does, and our alpha is overestimated relative to MARCS at the energy "edges", suggesting this is true. Another possibility is that I have some sort of bug which "compresses" the cross-section for each electron configuration in wavelength space, making it spikier (???)

I'm fairly confident that we aren't "compressing" the cross-sections in wavelength space. A while back I proved that we were doing things correctly by comparing OP opacities for hydrogenic species to the analytic formula. (If there is a bug, it's related to multi-electron atoms or it appeared after I performed that test)

I'm fairly confident the issue is related to smoothing. I seem to recall that the MARCS paper explicitly states that they smoothed opacities. I seem to recall that they were vague about exactly how they smoothed them (but I think they implied that they wanted to smear to remove resonance spikes).

It might be worth checking what happens if you artificially inflate the smoothing scale just to see if that helps explains the difference. Alternatively, it could also be worth checking what SME does here...

[ajwheeler]

This makes sense to me. I added in iron as well, and got this:

I can definitely imagine that adding back all the species in MARCS and smearing them a bit will get up pretty close to the MARCS opacs. (Though they do include some lines, so it's OK if they still have more absorption.)

[ajwheeler]

I'm going to close this since the topic of conversation has changed significantly since I started it and the original issue is solved. Future discussion of continuum absorption can happen in more specific issues.

[ajwheeler]

It's log(optical depth at 5000 Angstroms).

On Wed, Sep 7, 2022 at 11:00 AM Matthew Abruzzo @.***> wrote:

That is fascinating - although the y-axis is a little confusing. Could you remind me what that means again?

— Reply to this email directly, view it on GitHub https://github.com/ajwheeler/Korg.jl/issues/83#issuecomment-1239507011, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAFN2G2MPVWVMN37YEL2FDTV5CUX5ANCNFSM5GRJ6TQQ . You are receiving this because you authored the thread.Message ID: @.***>