ajwheeler / Korg.jl

fast 1D LTE stellar spectral synthesis
https://ajwheeler.github.io/Korg.jl/
BSD 3-Clause "New" or "Revised" License
44 stars 7 forks source link

Behaviour next the balmer break for hotter stars #149

Closed segasai closed 10 months ago

segasai commented 1 year ago

Hi,

While looking at synthesised spectra of hotter stars I'm seeing weird behaviour next to the Balmer break like here (spectrum rebinned by a factor 100)

image

It looks like high order Balmer lines don't produce enough opacity here. As far as I can see the linelist has lines down to 3648.26000, so that's not the problem.

This was the atmosphere: p7500_g+5.0_m0.0_t01_st_z-0.25_a+0.10_c+0.00_n+0.00_o+0.10_r+0.00_s+0.00.mod

I looked at the code looking for some kind of truncation on Hydrogen lines, but couldn't immediately identify it.

ajwheeler commented 1 year ago

Korg actually ignores the hydrogen line in the linelist (which is why there's a hydrogen_lines switch kwarg for synthesize). It looks like this is happening because the stark broadening profiles for the hydrogen lines don't exist beyond particular electron number densities.

The Balmer break for these stellar parameters according to Korg and PHOENIX, with the locations of the H lines marked. (This was for my own sake--I realize it's clear to you.) image

And the $\tau_{5000}$ in this atmosphere at which the Stark-broadening profile for each line ceases to exist and Korg gives up: image

Pretty bad! Korg is "officially" only for FGK stars, and it's largely untested for hotter and colder stars, which was a decision we made just for the sake of getting the paper out without unmanageable delay. To be honest I would not be surprised if you encounter other problems with hot stars, but you are welcome to persevere and I will do my best to address them. I'm going to see what I can do about this by using regular voigt profiles outside the currently supported regime.

segasai commented 1 year ago

Thanks for looking into this. I guess the F star boundary is technically at ~7000K :) and the issue exists there as well, that's why I felt okay to report. I think I also saw something similar in the Paschen break for not so hot stars. I'll try to see if I can understand what's going there to suggest a fix (something that comes to mind -- not the drop the line but use the last point in the interpolation grid )

ajwheeler commented 1 year ago

This is good point, and you should feel free to report issues at all temperatures anyway. I will check the Paschen break as well.

something that comes to mind -- not the drop the line but use the last point in the interpolation grid

This is a good idea.

On Mon, Nov 21, 2022 at 12:55 PM Sergey Koposov @.***> wrote:

Thanks for looking into this. I guess the F star boundary is technically at ~7000K :) and the issue exists there as well, that's why I felt okay to report. I think I also saw something similar in the Paschen break for not so hot stars. I'll try to see if I can understand what's going there to suggest a fix (something that comes to mind -- not the drop the line but use the last point in the interpolation grid )

— Reply to this email directly, view it on GitHub https://github.com/ajwheeler/Korg.jl/issues/149#issuecomment-1322443751, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAFN2G2NN37P3PWHDSNFBJLWJOZSHANCNFSM6AAAAAASFQNLZA . You are receiving this because you commented.Message ID: <ajwheeler/Korg. @.***>

ajwheeler commented 1 year ago

I'm a bit torn here about which strategy to use. The issue with using the profile for the largest electron number density in the grid is that it is often too small by orders of magnitude. On the other hand, the naive treatment of the H lines with standard Voigt profiles doesn't seem to work very well. Here's one of the Balmer lines at several electron number densities: image The solid lines are from the Stehle & Hutcheon 1999 table, and the dashed lines are the standard treatment.

I'm going to think a bit more about this, and return later to make sure that I didn't make any silly mistakes with that plot.

edit: I fixed a mistake which made the agreement between profiles look worse, but general outlook is the same.

ajwheeler commented 1 year ago

edit: these are plotted with the wrong PHOENIX spectrum. Korg and PHOENIX agree until Korg goes wonky close to the break.

For some additional context, this is the current behavior (the phoenix atmosphere off by 100K, but that's the least of our worries): current

This is what I get if I use the profile for the largest electron number density available: const extrap

And if I use totally naive voigt profiles with the default broadening approximations. This actually has more H lines than the implementations above, but they appear more attenuated close to the break (I'm looking into this). download-1

mabruzzo commented 1 year ago

Fascinating! Is PHOENIX directly using Stehle & Hutcheon 1999?

Out of curiosity, do we have any idea what the other codes are doing in this regard?

As an aside, do you think this causes the differences in the continuum near the Balmer break when compared to turbospectrum (I seem to recall that they include H-lines as part of the continuum)? Or is that comparison at an irrelevant temperature?

image

ajwheeler commented 1 year ago

I have not yet tried to figure out what PHOENIX is doing. Turbospectrum uses HLINOP, which has the same basic setup as Korg (Stehle and Hucheon '99 for stark broadening and Barklem, Piskunov and O'Mara for self-broadening). That said, it does some other stuff too (e.g. Mihalas Hummer Dappen formalism), and it and Korg don't really agree on H lines.

As an aside, do you think this causes the differences in the continuum near the Balmer break when compared to turbospectrum (I seem to recall that they include H-lines as part of the continuum)? Or is that comparison at an irrelevant temperature?

I never figured out what is going on with the Balmer break in Turbospectrum. Honestly, I mostly stopped worrying about it when I saw that SME didn't reproduce it.

andycasey commented 1 year ago

I think this is being caused by upper energy levels being perturbed or dissolved in dense environments.

This would be important for the hydrogen partition function in dense and hot (>10,000 K) regions. The inner layers of the photosphere that Sergey highlighted exceed 10,000 K. We talked about including the probability occupation formalism, but I don't think it was ever implemented, right? I think there was some uncertainties about what exactly "The Right Thing To Do(tm)" was, but I don't remember what there was uncertainty about.

To test this, I used the MARCS model that Sergey originally referenced to synthesize spectra from 3400-3900 Angstroms (at 0.01 Angstrom steps) with SME, TurboSpectrum, and Korg.

image

Two experiments are shown:

  1. Continuum only: "no" transitions (with hydrogen_lines=false where applicable).
  2. Only hydrogen lines.

The "continuum only" test is a bit misleading because TurboSpectrum will load in the hydrogen line data even if you tell it not to. In any case we can see that TurboSpectrum seems to be the most realistic in both cases, and something funky is going on with SME in the lower panel. Ignore the SME spectra for now because it is a distraction. We'll come back to it later.

I found nothing informative by extracting the opacities from TurboSpectrum, so I dug into the code instead. I found that TurboSpectrum is computing the occupation probability for H. (It's also doing some hacky stuff nearby but that doesn't explain what is going on here).

Below is the original continuum spectrum from TurboSpectrum for the same photosphere, and in the dotted line is when I use TurboSpectrum to calculate the continuum spectrum for the same photosphere except here I have disabled the occupation probability calculations for hydrogen (by just forcing WCALC to return 1.0 always).

image

The dotted line matches the behaviour we see in Korg and SME, which is consistent with my suspicion that all of this is due to upper energy levels being perturbed or dissolved in dense environments.

Now, back to SME:

The spectra look very unphysical. The SME spectra in the 'hydrogen' experiment look similar to the strange behaviour we saw between 3660-3670 Angstrom with other model atmospheres. The hydrogen lines in the SME spectra also look slightly offset. I checked and SME has read in the line list correctly and knows the medium (air) that the wavelengths are given in, so it doesn't seem to be a vacuum/air wavelength issue. My guess is there is something wrong in SME's implementation of the probability occupation formalism, but that wouldn't explain the wavelength offsets... In short: I don't know what the hell is going on here, but whatever is going on in SME is unrelated to this issue.

ajwheeler commented 1 year ago

First of all, thank you for figuring this out. I'm embarrassed to have not realized it given that I knew that

It should have been obvious putting these things together it should have been obvious what causes the "roll off" before the Balmer break in turbospec.

This would be important for the hydrogen partition function in dense and hot (>10,000 K) regions. The inner layers of the photosphere that Sergey highlighted exceed 10,000 K

Unfortunately, for the upper energy levels involved near the Balmer break, the MHD formalism is important at lower temperatures: Screen Shot 2022-11-28 at 11 41 52 AM $\tau{5000} \approx 10^{-3}$ corresponds to $T\mathrm{eff} \approx 500 \mathrm{K}$ in the sun. It's impact on the partition function in this regime, so it doesn't affect chemical equilibrium, but the lines themselves are affected.

 I think there was some uncertainties about what exactly "The Right Thing To Do(tm)" was, but I don't remember what there was uncertainty about.

Basically, the problem is that the MHD formalism, while probably better than naive Saha-Boltzmann, sometimes over-attenuates the outer energy levels (See Figure 9 in the Korg paper). Anil Pradhan tells me this is a known thing among solar people, who compare it to a fully "non-chemical" (i.e. electrons/nuclei are treated "directly") equation of state. I prefer things which are more wrong in I ways I understand than things which are less wrong in subtle confusing ways, so I decided not to roll the MHD formalism into Korg. I'm very willing to reconsider, though.

andycasey commented 1 year ago

I neglected to write this comment earlier:

I think we should do what Turbospectrum does, even if it means the outer energy levels will be attenuated. The lines for the outer energy levels are already "wrong". This would just make them wrong in a different way, but would make the Balmer break more realistic.

ajwheeler commented 1 year ago

To give an update on this, I've started working on an implementation of H I bf opacities that incorporates level dissolution via MHD, but there are a couple things holding fix back. Most importantly, I'm not sure of the best way to calculate photoionization cross section redward of the jump. Turbospectrum (HBOP) extrapolates the tabulated cross sections and applies an ad-hoc "tapering" because they don't decline fast enough. It's possible that this is the best way forward, but I want to explore other options first.

Relatedly, I think it may be worth it to upgrade the hydrogenic photoionization cross-sections in general (probably to Nahar 2021). I'm going to talk to the atomic people here about the possibility, but it may take a couple of iterations to get to an understanding of the tradeoffs.

@segasai, is this something which is preventing you from moving forward? I'm going to continue to work on it in the background while other things take priority, but I will move it up to TODO list if it's a bottleneck for you.

segasai commented 1 year ago

@ajwheeler I was reluctant to say that earlier, but yes, with the current balmer/pashen break behaviour I don't think I could use synthesized korg spectra. A few cases where I was aiming to try korg were low res cases with broad wavelength coverage, so hydrogen is pretty visible there. (full disclosure -- I don't have a specific goal in mind, I was aiming to see how comparable the spectra are to phoenix and/or I can substitute one for another , and whether I can use korg to either get some individual element abundances or understand possible accuracy in things like alpha and what's recoverable from various LR spectra I have access to).

ajwheeler commented 1 year ago

I returned to this today, and noticed that the plots above are using the incorrect PHOENIX spectrum. This was a silly mistake on my part. The codes are in much closer agreement than I thought, except for Korg's "emission" directly red of the jump. This is at Teff=7500 (Korg/MARCS), 7600 (PHOENIX), logg=5.0, solar abundances: image

It's still the case that neither simple solution (using the closest available stark profile, or using naive Voigt profiles) fixes. The best solution still seems to be including bf opacity redward of the jump, possibly with a fudge to account for MHD being over-eager and/or the cross-sections being wrong.

I also wanted to check what real data looks like. This shows the Balmer jump for slightly different parameters according to Korg, PHOENIX, and MILES. Again, the parameters are slightly different for each code/library, so this is just a sanity check. Teff ~= 7500, logg ~=3.5, solar abundances. image

ajwheeler commented 1 year ago

I have a working implementation of H I bf that fixes the problem raised in this issue, but I am not yet satisfied with the correctness of the implementation. I wrote a prototype implementation of the H I bf absorption cross-section which dissolves orbitals according to the MHD formalism. Here's sigma as a function of wavelength. (The vertical lines mark where the cross-sections are no longer tabulated for the n=1 orbital, and where the tapering kicks in.) image

As you can see, the "taper" fudge used by turbospec makes a big difference. Interesting, the "extra absorption" in the uncorrected MHD case comes almost entirely from n=1 orbital, i.e. the adjusted Lyman break (this is noted in the turbospectrum source). In these calculations the taper is applied only to the Lyman break. I would prefer to get the cross-sections right, but the fact the details of this fudge matter the most below 2000 \AA is encouraging.

Here's some synthesize spectra (parameters the same as in the previous post): image Unfortunately, it looks like the tapering the Lyman bf absorption impacts things significantly far from the Lyman break. Turning the tapering off makes the flux blue of the balmer break accurate, but it's not clear whether or not it's helping or hurting in the red.

ajwheeler commented 1 year ago

I think #170 fully fixes this. It's merged now, so if you update Korg you shouldn't have this problem.

segasai commented 1 year ago

Thanks! I'll give it a try. I'm rerunning the full synthesis of the full marcs grid now.

ajwheeler commented 10 months ago

Closing this because I think it's solved, but please reopen if I am wrong.