ajwheeler / Korg.jl

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

Comparison to PHOENIX in the CaT region #139

Open segasai opened 1 year ago

segasai commented 1 year ago

Hi,

I've synthesized a bunch of spectra in the 3000A to 10000A wavelength range using Korg on a logg/teff/feh/alpha grid that I can compare with PHOENIX. One of the features I observe is major difference around the CaT (in the wings particularly). compar here's a T=4500K, logg=5, [Fe/H]=0 [a/Fe]=0 case I understand that you are not responsible for PHOENIX, but I heard a view from an expert )) that we see here is more likely korg's deficiency rather than PHOENIX due to some kind of strong line treatment. Do you have any advice/opinion on that ? Or maybe I'm not using the right options ?

Thanks

PS I just used the standard Korg synthesize call with MARCS models + all the lines from VALD + individual alpha element abundances specified according to [Fe/H], [a/Fe].

ajwheeler commented 1 year ago

that we see here is more likely korg's deficiency rather than PHOENIX due to some kind of strong line treatment. Do you have any advice/opinion on that ?

My expectation is that this is correct, but I'm going to have to do some reading about PHOENIX's treatment of strong lines. Thanks for bringing this to my attention.

segasai commented 1 year ago

The PHOENIX paper says this on strong lines (sec 2.2 https://www.aanda.org/articles/aa/full_html/2013/05/aa19058-12/aa19058-12.html ) "The next step was to activate a mode in PHOENIX that triggers the use of special line profiles, which we used for the Ca lines, for example. In this mode no convergence is checked, so we carried out five more iterations by default in order to stratify the atmosphere again properly. Finally, the high resolution spectrum was synthesised from this intermediate model." Pointed by @alexji

ajwheeler commented 1 year ago

Thanks.

I've been skimming the various PHOENIX papers, and I'm still unsure what they mean by "special" profiles. Several papers mention that they use Doppler profiles for weak lines, and Voigt profiles for sufficiently strong lines, but it doesn't seem like that's what's being described here.

It's worth nothing (this is largely a note for my future self) that in section 2.3 of that paper, they say that they treat Ca II in NLTE. That is likely to be more accurate than Korg, with the big caveat that NLTE calculations involve even more uncertain physical quantities than LTE ones, so systematics may be very important. It's almost certainly not what's going on here, though, since NLTE will affect the line core much more strongly than the wings.

This may also come down to a difference between the MARCS and PHOENIX atmospheres, though the continua agree, so I don't think that's likely either. image

andycasey commented 1 year ago

The description of how strong lines is treated is a little short. One other potential point of difference here is that I think the PHOENIX grid is entirely spherical geometry, and the MARCS models in this parameter range are all plane-parallel. But I doubt that would explain these differences.

segasai commented 1 year ago

I googled for 'special line profile' and could only find this reference https://iopscience.iop.org/article/10.3847/1538-3881/ab8d33 which also talks about special line profiles, but without going into details. I'm going to ping @alexji again since he mentioned previously that he has seen aome strong line flag in MOOG. (i have no clue myself).

andycasey commented 1 year ago

The strong line flag in MOOG is used to tell MOOG that it should calculate opacities for a given line at all wavelengths. Otherwise MOOG calculates line contributions +/- some user-defined window in wavelength (usually ~2 Angstroms). Korg will calculate opacities until the line contribution relative to the continuum is small (as you saw with Balmer lines), which should work well for weak and strong lines. If it didn’t work well for CaT then we’d see the sharp cutoff in wavelength like you saw elsewhere. I think the “special line treatment” here (whatever it is) might be different. I’ll read the paper you linked.

ajwheeler commented 1 year ago

Agreed that spherical geometry is unlikely to be the culprit. This paper lists H lines as the only ones with "special" profiles, but they be leaving something out.

andycasey commented 1 year ago

@segasai do you have the PHOENIX code, or just the pre-computed spectra using PHOENIX?

segasai commented 1 year ago

I don't have PHOENIX code. Just the v2.0 synthesized spectra.

ajwheeler commented 1 year ago

They don't make their atmospheres (T, P, etc at a function of tau or depth) available, do they? I couldn't find anything online, but I may have missed something.

segasai commented 1 year ago

Here ftp://phoenix.astro.physik.uni-goettingen.de/AtmosFITS/ they have tau,temp,pgas,rho for 64 layers in their atmospheres and even element abundances

ajwheeler commented 1 year ago

Andy tells me that these don't have electron partial pressures, which Korg relies on, so I can't compare Korg+PHOENIX atmosphere vs Korg+MARCS. I'm hoping to remedy that soon, so we will be able to do the comparison eventually, just not yet.

andycasey commented 1 year ago

Out of interest, here is a comparison of the photospheric quantities. The models used were:

So both have Teff = 4500, logg = 5.0, [M/H] = 0, [alpha/M] = 0, microturbulence = 1 km/s. By selecting PHOENIX models with [M/H] = 0 and [alpha/M] I could only find the "lte*" photospheres.

I assume \tau in PHOENIX files are Rosseland optical depths. Here's the comparison:

phoenix

The PHOENIX models go much shallower than the MARCS models, but the Ca II triplet lines form much deeper, so that's an unlikely explanation for the difference.

The PHOENIX documentation page (https://phoenix.astro.physik.uni-goettingen.de/?page_id=108) says that they store electron partial pressure, but the FITS file I found on the FTP server definitely doesn't:

image

(I think there's also a unit bug in the PHOENIX headers for stellar masses, because it implies the PHOENIX model I used is 600 M_\odot).

alexji commented 1 year ago

I think the easiest/best thing to do would be to contact one of the PHOENIX authors if they would be willing to send a few model atmospheres with the required quantities (and ask about the Ca treatment as well).

For the MOOG stuff: there was an ancient (no longer used) flag to do magic things with damping coefficients for specific strong lines. I suspect something similar happens in PHOENIX to better reproduce observed lines. As Adam said I think they use NLTE for some specific things too, but I think the damping treatment is more important here as it's a wing mismatch.

segasai commented 1 year ago

I gave a wrong link. this is the link to the v2.0 atmospheres ftp://phoenix.astro.physik.uni-goettingen.de/v2.0/AtmosFITS/ and these have the pe column

andycasey commented 1 year ago

With v2.0 atmosphere. Seems to be better agreement if \tau is at 5000 Angstroms. I'll try synthesize something now.

image
segasai commented 1 year ago

I've just briefly looked at PHOENIX vs actual data. I looked at stars from APOGEE with S5 1700D data. This is a star with logg=4.566099, Teff=4992.816, M/H=-0.08617 according to APOGEE DR17. and the 2DF data on that star was fitted by my PHOENIX interpolated models with teff=5250, logg=5, feh=-0.3 And here is the fit s5star

Red is PHOENIX, black is data. The cores are not well fit, but wings look sensible, so I'd say that supports the view that PHOENIX is accurate here.

andycasey commented 1 year ago

Unfortunately we still need (physical) depth in these photosphere models in order to synthesise a spectrum with Korg.

I'm also inclined to believe that PHOENIX is more accurate here, but we will need to understand why.

Out of interest, if you synthesize a Korg spectrum at the two possible sets of parameters (4993/4.57/-0.09 vs 5250/5.0/-0.3) and convolve them to match the data, how do those compare?

segasai commented 1 year ago

To answer the last question, I finally managed to substitute PHOENIX by Korg grid for my spectral fitter (rvspecfit), so I can now side to side compare interpolated Korg/PHOENIX. as well as refit say S5 data with Korg grid. For the example star show above here're the plots for the two sets of parameters (4993/4.57/-0.09 vs 5250/5.0/-0.3) The spectra are convolved to 1700D resolution. image

segasai commented 1 year ago

FWIW

here's a side-by-side comparison of fits by interpolated PHOENIX vs vs interpolated Korg for the same star as shown above: We can see the same CaT problems as discussed above. And also Alex pointed at the 8806 Mg line. (Top is phoenix, bottom is Korg) fig_phoe fig_korg

ajwheeler commented 1 year ago

Korg doesn't need depth when doing calculations in plane-parallel, which a good approximation at logg=5. I've synthesized spectra using the two atmospheres compared by Andy, and the difference looks very close to Sergey's original plot. image

segasai commented 1 year ago

I've looked at the literature and there is a substantial number of papers on CaT specifically. They talk quite a bit about collisional broadening with neutral atoms and van der Waals: https://ui.adsabs.harvard.edu/abs/2000A%26A...353..666C/abstract https://ui.adsabs.harvard.edu/abs/1990A%26A...231..125S/abstract https://ui.adsabs.harvard.edu/abs/1988MNRAS.231..115S/abstract https://ui.adsabs.harvard.edu/abs/1992A%26A...254..258J/abstract And in these papers they do generally get reasonably looking CaT

I'm wondering if there could be some issue with this code https://github.com/ajwheeler/Korg.jl/blob/3716575e67e6cd67ac55f7de3e8f86914bda4de3/src/linelist.jl#L291 as it looks quite dense : )

alexji commented 1 year ago

I took a quick look at the parser, it should be getting the VDW coefficient from the line list so not using the default Unsold approx.

Sergey, I think we should check your VALD linelist to make sure it was queried correctly; it looks like you may not have the ABO theory damping coefficient. For the 8662 line, I see your VDW values are -7.675ish, while the ABO value should be 291.275 (Adam has correctly implemented parsing this as far as I can tell; details here https://www.astro.uu.se/~barklem/howto.html). I don't know how big of an effect this could be or whether it is sufficient to explain the difference, but at least we should fix that.

ajwheeler commented 1 year ago

Here's a zoomed-in version of the plot above showing that the atmospheric structure makes a huge difference, but doesn't fix the discrepancy--it makes it just as large in the opposite direction. The Ca II line in the Korg + PHOENIX spectrum is wider than the PHOENIX spectrum instead of narrower. I agree that the broadening parameters are likely to be important, which is why I'm very curious about the potential "special treatment" in PHOENIX.

image

I took a quick look at the parser, it should be getting the VDW coefficient from the line list so not using the default Unsold approx.

This is correct. The ABO calculation code is here if you are curious. I wouldn't bet my life on its correctness, but I think it's right. The linelist I've been using (including for the plot here) has the correct ABO parameter, so there is still unexplained disagreement.

alexji commented 1 year ago

Based on Peter's e-mail, can you maybe try running the 8662 (air) line, with damping values -7.760 instead of 291.275 (matching what VALD says is the Kurucz damping data instead of ABO)?

alexji commented 1 year ago

For the ABO damping, here's where it is compared to turbospectrum. https://github.com/bertrandplez/Turbospectrum2019/blob/master/source-v19.1/depth.f#L193

The exponent of (vbar/v0) is maybe different (^-alpha vs ^(1-alpha)), though I might be missing a term elsewhere that cancels it out. https://github.com/ajwheeler/Korg.jl/blob/main/src/line_absorption.jl#L241 https://github.com/bertrandplez/Turbospectrum2019/blob/master/source-v19.1/depth.f#L200

ajwheeler commented 1 year ago

I don't see anywhere the the -alpha becomes (1-alpha) in Turbospectrum, but I think it's correct. (see https://www.astro.uu.se/~barklem/howto.html). The fortran implementation on that page also uses (1-alpha):

GV=(4./PI)**(XALPHA*0.5)*GAMMAF*1.E4*SIGMA
VBAR=SQRT(21173.*T*(1./1.008+1./(AWGT/UU)))
GV=GV*((VBAR/1.E4)**(1.-XALPHA))

ABO vs no ABO does not make a big difference. Here's the above plot showing both Korg spectra computed with -7.760 and 291.275 for the Ca II line. ABO vs no

ajwheeler commented 1 year ago

To get a better perspective on the tension, I've plotted here one of the CaT lines for 6 different Teff/logg pairs. I've plotted for each PHOENIX, Korg + MARCS, Korg + PHOENIX atmosphere, and Korg + PHOENIX atmosphere with my best guess at the PHOENIX line parameters. I'm also matching the PHOENIX microturbulence and solar abundances. phoenix vs korg

I'm pulling the Kurucz line parameters from here, which may be wrong, but I don't know where to find the specific version Peter said they used. If I understand his email correctly, the "special treatment" of the CaT lines is just the ABO parameters from Barklem, Piskunov & O'Mara (2000), which are the same as what is currently on VALD (at least for the Ca triplet--I just checked.)

I think it's clear that the choice of model atmosphere is very important. The specific line parameters are too, and if anyone knows where I can get my hands on the 2005 Kurucz data, please let me know. Other possible sources of disagreement are the continuum opacities and chemical equilibrium solvers. I didn't see details about either of those in Husser+ 2013, but I haven't done a lot of digging. I think it's unlikely that this comes down to the radiative transfer scheme, line broadening implementation (as opposed to input data), or NLTE effects.