Brackett line profiles, normalization, and level dissolution

[ajwheeler]

In #200, I added support for the Brackett series based heavily on Barklem's HLINOP.f. My Korg.brackett_line_profile and STARK1 in HLINOP.f are doing the same calculations (though I have stripped out a lot as I only need to calculate profiles for Brackett lines). I have confirmed by running the fortran directly that I am producing the same values, and Korg produces spectra with Brackett lines nearly identical to Turbospecrum (see below), which uses HLINOP.

Here I show synthesized spectra in the APOGEE region (using the APOGEE linelist and model atmospheres) for a few sets of stellar parameters. All have [M/H] = -1.5. The last column shows the probability that the upper levels of the two zoomed-in lines are dissolved. (Tau=1 is marked with a vertical line.)

(Click to view bigger.)

I think this shows clearly that in the high logg-cases, Korg and turbospectrum show less absorption because they are dissolving the upper levels of the transitions. I was under the impression that synspec did the same, so it's worth tracking down the difference.

The Korg and Turbospectrum Brackett lines implementations are quite different from the synspec Brackett lines. Synspec's implementation of the profiles is simpler, having only a Doppler core and an asymptotic Holtsmark wing. The Korg/Turbospec implementation dynamically chooses how to model the core (Stark often dominates over Doppler), has a more sophisticated Holtsmark profile, and include the effect of impact broadening by electrons.

Click to view the synspec stark profile source code

```fortran C C ******************************************************************** C FUNCTION STARKA(BETA,A,DIV,FAC) C C Approximate expressions for the hydrogen Stark profile C C Input: BETA - delta lambda in beta units, C BETAD - Doppler width in beta units C A - auxiliary parameter C A=1.5*LOG(BETAD)-1.671 C DIV - only for A > 1; division point between Doppler C and asymptotic Stark wing, expressed in units C of betad. C DIV = solution of equation C exp(-(beta/betad)**2)/betad/sqrt(pi)= C = 1.5*FAC*beta**-5/2 C (ie. the point where Doppler profile is equal to C the asymptotic Holtsmark) C In order to save computer time, the division point C DIV is calculated in advance by routine DIVSTR. C FAC - factor by which the Holtsmark profile is to be C multiplied to get total Stark Profile C FAC should be taken to 2 for hydrogen, (and =1 C for He II) C INCLUDE 'PARAMS.FOR' c PARAMETER (F0=-0.5758228,F1=0.4796232,F2=0.07209481,AL=1.26) c PARAMETER (SD=0.5641895,SLO=-2.5,TRHA=1.5,BL1=1.14,BL2=11.4) PARAMETER (F0=-0.5758228,F1=0.4796232,F2=0.07209481/2.,AL=1.26) PARAMETER (SD=0.5641895,SLO=-2.5,TRHA=1.5,BL1=1.52,BL2=8.325) PARAMETER (SAC=0.07966/2.) XD=BETA/BETAD C C for a > 1 Doppler core + asymptotic Holtzmark wing with division C point DIV C IF(A.GT.AL) THEN IF(XD.LE.DIV) THEN STARKA=SD*EXP(-XD*XD)/BETAD ELSE STARKA=TRHA*FAC*EXP(SLO*LOG(BETA)) END IF ELSE C C empirical formula for a < 1 C IF(BETA.LE.BL1) THEN STARKA=SAC*FAC ELSE IF(BETA.LT.BL2) THEN XL=LOG(BETA) FL=(F0*XL+F1)*XL STARKA=F2*FAC*EXP(FL) ELSE ```

In summary, Synspec has simplified Stark broadening, and it's doing somthing differently with level dissolution. So, why think there's anything further to discuss? When I look directly at the profiles produced by Korg.brackett_line_profile or STARK1, they do not integrate to unity. They usually integrate to something between 1.8 and 2, with the lines with smaller upper principal quantum numbers integrating to almost exactly 2. This is because they are summing two profiles as an approximation of convolution: one for the impact broadening by near electrons, and one for the quasistatic Holtsmark broadening by ions and distant elections. (See these lines in Korg and this line in STARK1.) When the line core is set my the Doppler profile, this isn't a problem, but I'm concerned that this isn't often the case for Brackett lines. This plot shows the widths of the Doppler, Lorentz, and Stark profiles as a function of optical depth for a wide swath of MARCS atmospheres (5000 K < Teff < 8000 K, 0 < logg < 1, -2.5 < [M/H] < 1) for the n=4->11 transition.

Korg and Turbospectrum have slightly different approaches to approximate convolution, but both (if I am reading HLINOP correctly) will treat the core as pure-Stark when the Stark profile has the largest half-width.

[ajwheeler]

Here's a related experiment. I'm looking here at the n=4->11 line only for stars with lower logg and less level dissolution.

1. The first row is the same as in the above.
2. In the second, I'm dividing the stark profile by 2, but only in the core. (Doing it everywhere makes the wings too weak.) This gives lines closer to the synspec ones, but it seems like I might need to tweak my definition of the line "core".
3. In the third, I'm enforcing that the core always be doppler-shaped. I think this showing that for this line, the core isn't dominated by doppler broadening in much of the relevant parameter space.

(click to view bigger)

[ajwheeler]

Part of what's going on here is that Brackett lines form much deeper in the atmosphere than, e.g. the Balmer series. This means that the assumption that Doppler dominates the profile core breaks down. In the plot of Doppler, Lorentz, and Stark half-widths above, I show the half-width of each component of the n=4->11 profile (Br η) as a function of optical depth. But, what optical depth range actually matters?

Here, I show depth of formation as a function of wavelength for Br η, as well as H alpha, in contrast. Br η for below the (Roseland) tau=1 layer!

At this depth, the Stark profile is much broader than the Doppler profile, so the convolution looks almost like pure-Stark (even in the core). You can verify that stark broadening dominates by looking at the half-width plot above. Here, I'm plotting the broadening profiles (with stark broken into impact and quasistatic) at the depth of formation for Br η. (The only difference between the plots is linear vs log scale. Summing the two stark components gives a too-strong profile out to about 50 Å from the line core.

Let's examine the normalization as a function of depth. In this plot, we can see that each component (impact, quasistatic) of the stark profile is indeed roughly normalized. The sum is thus fairly close to 2, and Korg's calculated profile jumps to 2 at the depth when it switches from doppler to stark for the line core ($\tau = 2 \times 10^{-2}$ or so.)

I think the most promising way forward may be to do actual numerical convolution of the line profiles. Experimentation with this line reveals that we can get acceptable accuracy if the convolution is done at a 0.1 Å resolution. This takes about $10^{-4}$ s per profile, which is about 0.25s per synthesis, assuming 50 layers and 50 Brackett lines. That seems accepatable, though implementation isn't trivial.

[ajwheeler]

This was fixed by #218.

Here's an updated comparison plot with my fixed profiles and MHD off:

The synspec profiles are pretty close to mine, which makes sense, because they are normalized if not totally correct (I claim).