Unexpected behavior of asymptotic period spacing

[earlbellinger]

Describe the bug Currently, MESA defaults to computing the asymptotic period spacing at $\nu_\max$ (see get_delta_Pg in star_utils). This is probably appropriate for red giant stars, where turbulent convection excites oscillations around $\nu\max$, but probably not appropriate for $g$-mode pulsators, whose modes may be excited in a different range by a different mechanism. This default behavior also means that MESA does not compute the asymptotic period spacing for many models, such as for a solar model (see example below), because the entire g mode cavity is below $\nu\max$.

The output also appears to be in general incorrect, as I will show below.

To Reproduce Include delta_Pg in history_columns.list.

Expected behavior $\Delta \Pi = \frac{2\pi^2}{\sqrt{2}} \left( \int{N^2 > 0} \dfrac{N}{r} \rm{d}r \right)^{-1}$ instead of $\Delta \Pi = \frac{2\pi^2}{\sqrt{2}} \left( \int{N^2 \geq (2\pi\nu_\max/10^6)^2} \dfrac{N}{r} \rm{d}r \right)^{-1}.$ There is also some additional behavior in MESA's computation of $\Delta\Pi$ which I do not really understand and may contribute to the incorrect figure as seen below.

Proposed solution Change the behavior of delta_Pg_mode_freq to the following:

1. delta_Pg_mode_freq $< 0$: use $\nu_\max$
2. delta_Pg_mode_freq $\geq 0$: use delta_Pg_mode_freq

Alternatively, we could change the default value of delta_Pg_mode_freq to 1d-99 (i.e., essentially 0, but not triggering the flag to use $\nu_\max$).

Example Inlist:

&star_job
    create_pre_main_sequence_model = .true.
/

&eos
/

&kap
  use_Type2_opacities = .true.
  Zbase = 0.02
/

&controls
    initial_mass = 1
    initial_z = 0.02
    xa_central_lower_limit_species(1) = 'h1'
    xa_central_lower_limit(1) = 1d-3
/

History columns file:

star_age
delta_Pg

History file output: (always zero)

                                      1                                        2 
                                star_age                                 delta_Pg 
                 1.0000000000000001E-005                  0.0000000000000000E+000 
                 7.4416000000000014E-005                  0.0000000000000000E+000 
                 2.5958682111999997E-004                  0.0000000000000000E+000 
                 7.2035107872931814E-004                  0.0000000000000000E+000 
                 1.8668799962237366E-003                  0.0000000000000000E+000 
                 4.7198108322034473E-003                  0.0000000000000000E+000 
...

If I now add delta_Pg_mode_freq = 1d-99 to controls, we get:

                                       1                                        2 
                                star_age                                 delta_Pg 
                 1.0000000000000001E-005                  1.2664938046875831E+007 
                 7.4416000000000014E-005                  1.2665604998279545E+007 
                 2.5958682111999997E-004                  1.2665603301258113E+007 
                 7.2035107872931814E-004                  1.2665607943531716E+007 
                 1.8668799962237366E-003                  1.2665607199246502E+007 
                 4.7198108322034473E-003                  1.2665605975245405E+007 
...

But if I then compute the periods using GYRE, calculate the period spacing, and compare with the integral I give above, I find that these values for delta_Pg are several orders of magnitude too high: Here the black dots come from GYRE, the red line is when I compute the integral over $N^2$ using the profile in the .GYRE file, and the red line is the output from MESA with delta_Pg_mode_freq = 1d-99. The latter is too high by a factor of 4350. Clearly something is amiss, but I don't know what.

I have also tried with a 5 $\rm{M}_\odot$ main-sequence star and found the same behavior, but the MESA output is off by a different factor. In this case I have also tried with the default value for delta_Pg_mode_freq and in both cases gotten the same wrong value: Now the output from MESA is too high by a factor of 1488. Since the factors are different it is unlikely to be a units issue (?). This issue may be related to these lines about the decay region as I don't have them in my calculation, but I am not sure.

[earlbellinger]

It seems this control has been discussed on the mailing list here: https://sourceforge.net/p/mesa/mailman/mesa-users/thread/B9CB631D-2B7A-498B-9525-983A46B3E805%40gmail.com/#msg32645250

I've pushed a fix in f3ee16d. I think the problem was that pow2(2*pi*1d-99/1d6) goes beserk instead of becoming ~0.

The result is improved but possibly not complete:

Solar model

Massive star

I'm a bit scattered as to why they disagree as the calculation should be basically exactly the same. In MESA we have:

omega2 = pow2(2*pi*nu_max/1d6)
if (nu_max <= 1d-99) omega2 = 0
el = 1
L2 = el*(el+1)

integral = 0
do k = 2, s% nz
   N2 = s% brunt_N2(k)
   r  = s% r(k)
   dr = s% rmid(k-1) - s% rmid(k)
   if (N2 > omega2) integral = integral + sqrt(N2)*dr/r
end do

delta_Pg = sqrt(2d0)*pi*pi/integral

and I've mimicked this in Python with

import pandas as pd
import numpy as np 
import os

prof = pd.read_table(os.path.join('work_dev', 'LOGS', 'profile5.data'), skiprows=5, sep='\s+')
N2 = prof.brunt_N2
N = np.sqrt([x if x>=0 else 0 for x in N2])
r = 10**mesaprof.logR

integral = 0
for ii in range(1, len(N)):
    integral += N[ii]*(r[ii-1] - r[ii])/r[ii]

delta_Pg = np.sqrt(2)*np.pi**2/integral

yet my Python estimate is clearly more accurate. Any ideas why?

[warrickball]

I think I have a diagnosis for all of this but won't have a chance to post until about ~7 hours from now.

[warrickball]

I sunk some time into this over the weekend and think I can see what is happening, though not why. I think we need to figure out why it's implemented the way it is and, if useful to keep, figure out precisely what we want to compute, and what controls we might need to signal that.

tl;dr

I agree that the default should reproduce the integral expression I think most people would expect (in particular, over N²>0). But for νmax ≈ 0 in stars with convective envelopes, the current implementation basically only integrates over the region with N²>0 near the surface, giving the absurd values. Before changing anything, I think we should try to figure out (or someone remember?) why 'delta_Pg' is implemented the way it is.

Details

To get everyone on the same page, here's a relevant snippet for what I think most people might expect code for ΔΠ to look like (in vals(1)):

         vals(1) = 0
         do k = 2, s% nz
            dr = s% rmid(k-1) - s% rmid(k)
            if (s% brunt_N2(k) > 0) then
               vals(1) = vals(1) + sqrt(s% brunt_N2(k))/s% r(k)*dr
            end if
         end do
         vals(1) = sqrt(2.0_dp)*pi**2/vals(1)

Here, I accumulate the integral in vals(1), then compute ΔΠ from π²√2/vals(1). Let's see what MESA actually does.

https://github.com/MESAHub/mesa/blob/9bae36e794ab4c87ce5c9d2754365fe415b5bc06/star/private/star_utils.f90#L791-L853

There are two important changes in the logic that I haven't used. First, it skips regions with ω² > S₁² (i.e., where p-modes propagate). Second, after passing (inward) over any g-mode region (N²>0), it starts accumulating a second integral with integrand √(-k), where k² = (1-N²/ω²)(1-S₁²/ω²)ω²/c² (which looks like the usual approximate wavenumber that we use to define the boundaries in propagation diagrams). When this second integral reaches some limit (hardcoded to 0.5), the whole loop exits.

I see two reasons why the g-mode calculation goes awry for a solar model. First, using ω²≈0 means that we don't skip the radiative near-surface layers (because ω² ≈ 0 ≤ S₁²), so that gets added to integral for ΔΠ. Then, the second integral starts accumulating and rapidly reaches the limiting value, so integral also stops accumulating and we exit with a small value for integral and thus a large ΔΠ. I get reasonable numbers by choosing delta_Pg_mode_freq to be anything greater than the Lamb frequency near the surface. e.g. 10 μHz worked.

The mystery to me, at the moment, is why the ΔΠ calculation is implemented like this. Like the previous thread Earl linked to, I presume it's something about red giants. I guess the second integral is trying to only select one g-mode region by avoiding connection to regions that are separated by too many decay-length-scales...? So I don't know if red giant asteroseismologists would want to preserve the current behaviour. I don't know offhand but will ask colleagues and collaborators about it.

[warrickball]

When I next have a chance, I'll make some red giant models, compute some dipole mixed modes and see how various estimates of ΔΠ compare with actual models.

If we do just integrate over regions with N²>0, should we neglect the near-surface regions with N²>0? They have to exist because the surface of the star must be radiative.

[warrickball]

I asked Dennis Stello reminded me of Bildsten et al. (2012), which defines the asymptotic period spacing as $$\Delta\Pi=\frac{\sqrt{2}\pi^2}{\int Nd\ln r}$$ over $\omega^2 < N^2, S_1^2$, rather than $0 < N^2$. The additional loop in the MESA implementation is probably related to the evanescent integral mentioned in the text, which is in turn possibly something about how well coupled the two g-mode cavities during the helium flash.

(The off-centre flash drives a convection zone, where $N^2 < 0$, which divides the usual g-mode cavity into two. See their Fig. 3.)

Aside from speculating on the observable properties of stars during the flash, I'm not sure how much weight this definition carries. I'd be interested to see a calculation of the actual period spacing from a mode frequency calculation compared to what one would estimate from different definitions of $\Delta\Pi$ but don't have time to do this myself.

I lean towards the default integral limits being the innermost region with $N^2 > 0$, which is easy enough to implement.

[evbauer]

Nice find! It seems plausible that the overly restrictive definition was motivated by this particular application, while our output should probably be the simpler more general definition (if I'm following this correctly).

[earlbellinger]

Thanks for tracking this down.

I'd be interested to see a calculation of the actual period spacing from a mode frequency calculation compared to what one would estimate from different definitions of $\Delta\Pi$ but don't have time to do this myself.

I plotted this in the original post (or did you have something else in mind?)

So perhaps I should add some documentation to my patch in https://github.com/MESAHub/mesa/commit/f3ee16daaa8707f71b4f3bbfb5f8b2b48d726636 and submit a pull request.

[warrickball]

I'd be interested to see a calculation of the actual period spacing from a mode frequency calculation compared to what one would estimate from different definitions of ΔΠ but don't have time to do this myself.

I plotted this in the original post (or did you have something else in mind?)

I mean in a red giant going through the flash, which was the application in Bildsten et al. (2012).

I think what'll happen in practice is that limiting the integral to the innermost g-mode cavity will have little effect on the original behaviour. So then it's just a case of changing the default to 0, and < 0 (not <= 0) to use nu_max, so that the old behaviour is still available, aside from trying to link multiple cavities depending on the "evanescent" integral, which doesn't appear to be documented anywhere.

So perhaps I should add some documentation to my patch in https://github.com/MESAHub/mesa/commit/f3ee16daaa8707f71b4f3bbfb5f8b2b48d726636 and submit a pull request.

I'm not sure your patch cover don't think your patch fully captures the issue. If we just set delta_Pg_mode_freq = 1d-99, then I think we end up integrating ΔΠ for the convectively-stable layer near the surface. IIRC I got sensible ΔΠ for a solar model only by choosing delta_Pg_mode_freq to be small but large enough for p-modes to propagate near the surface. I think we need to rewrite the integral because when we integrate where N²>0, I think we secretly neglect any convectively-stable near-surface layers. Loop outwards from the centre to the surface, start accumulating in the first g-mode cavity and return the result as soon as we leave it.

I'll try to revisit this in the next few days, now that the way forward is clearer.

[warrickball]

Here's a demonstration of why just patching the existing parameter isn't the whole solution. I've just computed a 1 Msun track with 5 definitions of ΔΠ.

1. MESA's current method with delta_Pg_mode_freq = 0 (so using νmax, "DPi_bildsten_numax").
2. MESA's current method with delta_Pg_mode_freq = 10 ("DPi_bildsten_10").
3. MESA's current method with delta_Pg_mode_freq = 1d-99 (so roughly 0, "DPi_bildsten_0").
4. An integral of N²dlnr over the inner g-mode cavity (first continuous region with N²>0, "DPi_outward") .
5. An integral of N²dlnr over the whole star ("DPi_everywhere").

Note the following.

• "DPi_bildsten_numax" doesn't appear because it's always zero.
• "DPi_bildsten_0" is orders of magnitude too large because the integral sees the near-surface radiative zone as a valid g-mode cavity, integrates the period spacing there, then gives up because it doesn't reach the central g-mode cavity before the evanescent integral crosses a threshold that says stop. This is why we can't just change the default value of delta_Pg_mode_freq.
• If we choose delta_Pg_mode_freq to be high enough that the near-surface radiative zone appears to be a p-mode cavity (in this case 10), the current implementation doesn't start integrating until it reaches the actual radiation zone. All is now well.
• If we integrate everywhere where N²>0, we get a small offset because our integral includes the near-surface layers.
• If we integrate by looping outwards until we find N²>0, start integrating, then exit when we encounter N²<0, we get what looks like the right number, but it seems to be unstable because sometimes on the pre-MS the first central cell has N²>0 but the neighbouring one doesn't, so the integral aborts on a small number and ΔΠ is huge. Specifically, this is my code:

       names(4) = 'DPi_outward'
       vals(4) = 0
       entered_g_mode_cavity = .false.
       do k = s% nz, 2, -1
          dr = s% rmid(k-1) - s% rmid(k)
          if (s% brunt_N2(k) > 0) then
             entered_g_mode_cavity = .true.
             vals(4) = vals(4) + sqrt(s% brunt_N2(k))/s% r(k)*dr
          else
             if (entered_g_mode_cavity) exit
          end if
       end do

  I think this integral over the innermost g-mode cavity is what we want (this is the definition in e.g. §7.4 of JCD's notes) but now there's this numerical behaviour to deal with. Not sure how best to work around this besides just accepting some strange behaviour, probably when there are small cores on the threshold of convective stability.

[fxt44]

... the first central cell has N²>0 but the neighbouring one doesn't, ...

maybe treat this case with an if statement? if s% brunt_N2(k) and s% brunt_N2(k-1) are both > 0

[warrickball]

I see this keeps coming up in the telecons so I'll chip in again here. IMO 'DPi_outward', for which I've copied my code above, is what we want, even though things are a bit weird on the pre-MS. I never tried Frank's suggestion.

The outstanding issue to me is what to do with the old behaviour. Do we add a flag for which formulation to use? 'traditional' vs 'bildsten2012' or something? Or keep the current behaviour for delta_Pg_mode_freq ≥ 0 and use different (negative) special values like -1 for the old formulation at $\nu_\mathrm{max}$ and -2 for the 'DPi_outward' I've defined above?