Estimating \(\rm \nu_{max}\) hacks

For first-time users, we’ll assume you do not know what \(\rm \nu_{max}\) is for a given star [and that’s totally ok]

In these scenarios, we have a convenient built-in method that uses the available information and/or data on a target to estimate a value for numax (\(\rm \nu_{max}\)). This is really helpful (and *strongly encouraged*) for our non-expert users. The primary reason for this is that the main module (i.e. the global fit) derives all parameters but models the stellar background and solar-like oscillations separately in two steps.

Basically we use the estimated \(\rm \nu_{max}\) to mask out the region in the power spectrum believed to be exhibiting the solar-like oscillations so that it will not influence the stellar background estimates. The reverse is also true, hence why we do this in two steps. Consequently, we correct for the background contribution in the power spectrum before measuring our global properties \(\rm \nu_{max}\) and \(\Delta\nu\).

---

\(\rm \nu_{max}\) hack summary

There are three main ways to brute force \(\rm \nu_{max}\) while still running the pysyd.target.Target.estimate_parameters method:

• **Option 1:** enter the desired trial number by using the --ask flag

• **Option 2:** provide your own value for \(\rm \nu_{max}\) by using the same flag

• **Option 3:** use an upper frequency limit to cut out sharp (likely not astrophysical) features in the power spectrum

---

from pysyd.utils import Parameters
from pysyd import plots
from pysyd.target import Target

Load in default settings for KIC 1435467

Since pySYD is optimized for command-line use as well as processing multiple stars, a lot of the options such as showing figures and printing verbose output are disabled. However for demonstration purposes, we’ll change two of the defaults.

params = Parameters()
params.add_targets(stars=['1435467'])
star = Target('1435467', params)
star.params['verbose'], star.params['upper_ex'] = True, None
if star.load_data():
    star.estimate_parameters()

-----------------------------------------------------------
Target: 1435467
-----------------------------------------------------------
# LIGHT CURVE: 37919 lines of data read
# Time series cadence: 59 seconds
# POWER SPECTRUM: 99518 lines of data read
# PS oversampled by a factor of 5
# PS resolution: 0.426868 muHz
-----------------------------------------------------------
PS binned to 219 datapoints

Numax estimates
---------------
Estimate 1: 7859.47 +/- 8.34
S/N: 12.19
Estimate 2: 7876.29 +/- 27.05
S/N: 6.73
Estimate 3: 1457.85 +/- 90.46
S/N: 9.25
Selecting model 1

By default, it selects the trial that has the highest signal-to-noise (S/N) detection – which here was the first estimate with S/N ~12. However, the value for numax is really high (\(\rm \nu_{max} \sim 7860 \mu Hz\)) and is likely latching on to an artefact or something else that is not astrophysical.

To be sure though, let’s take a look at the plots and results to see what’s going on.

Plot estimates

star.params['show'] = True
plots.plot_estimates(star)

As suspected, it selected and highlighted the first “trial” due to some high frequency artefact. It was still behaving how it was expected to, but that’s not the type of feature we are looking for. In the upper righthand corner, we can see power excess in the ~1000-2000 muHz frequency region.

We have a few options of how we can go about this. The first is to provide an upper frequency limit that will restrict the power spectrum to below that large spike. Although, if you take a look at the lower righthand corner, the third trial does end up estimating a good value for numax. Therefore we can brute force this trial with the --ask option, which is False by default.

---

Option 1: enter the trial number

We will re-load the star and this time change the --ask command to True. This option will literally ask you which trial you prefer. This will also display the plot regardless of your preset options, that way you can make the most informed decision for which to select. As a result, we do not recommend using this for many stars.

star.params['ask'] = True
star.estimate_parameters()

-----------------------------------------------------------
PS binned to 219 datapoints

Numax estimates
---------------
Estimate 1: 7859.47 +/- 8.34
S/N: 12.19
Estimate 2: 7876.29 +/- 27.05
S/N: 6.73
Estimate 3: 1457.85 +/- 90.46
S/N: 9.25

Which estimate would you like to use? 3
Selecting model 3

Now we can feel more confident about the value that we provide to the main method.

Before we go to our other main alternative (masking the higher frequencies), there’s one other thing I’d like to point out. Let’s assume that all three of these trials had bogus values but we are pretty confident that the numax was around ~1450. Turns out that we an use the same flag but instead, provide our own value.

---

Option 2: provide your own estimate

star.params['ask'] = True
star.estimate_parameters()

-----------------------------------------------------------
PS binned to 219 datapoints

Numax estimates
---------------
Estimate 1: 7859.47 +/- 8.34
S/N: 12.19
Estimate 2: 7876.29 +/- 27.05
S/N: 6.73
Estimate 3: 1457.85 +/- 90.46
S/N: 9.25

Which estimate would you like to use? 5
ERROR: please select an integer between 1 and 3
       (or 0 to provide your own value for numax)


Which estimate would you like to use? 0



What is your value for numax? 1400.00
Using numax of 1400.00 muHz as an initial guess

So we intentionally provided an integer value for a trial that does not exist so you could see the alternate option(s) for input data. By entering 0, we are able to provide a float number that we will set numax to.

The plot will still show the three trials but you might’ve noticed that it didn’t highlight any trials since we are using our own value.

Finally, we can catch this problem earlier on by providing an upper frequency limit for the power spectrum that is used in this routine.

---

Option 3: Provide upper limit

Another option is to set the upper frequency limit (--upper_ex) to something below the high-frequency artefacts, say \(\rm 6000 \mu Hz\).

Note: you’ll notice there are a lot of lower/upper bounds but their naming will start to make sense as you use the software more. For example, this first routine that estimates numax used to be called find_excess() (since that’s technically what it does!) and hence, the “ex” for the flag. The same is true for the BackGround-fitting routine (--lower_bg/--upper_bg)

star.params['ask'], star.params['upper_ex'] = False, 6000.0
star.estimate_parameters()

-----------------------------------------------------------
PS binned to 189 datapoints

Numax estimates
---------------
Estimate 1: 1430.02 +/- 72.61
S/N: 2.43
Estimate 2: 1479.46 +/- 60.64
S/N: 4.87
Estimate 3: 1447.42 +/- 93.31
S/N: 13.72
Selecting model 3

Remember that the figure was only displaying because we enacted the --ask option, so let’s double check the figure to be sure.

plots.plot_estimates(star)

If you look in the upper righthand corner of the figure, you’ll notice that the power spectrum only goes up to \(6000 \mu \mathrm{Hz}\) this time. That enabled three nice, equally-usable trial runs but you can see that the values for \(\rm \nu_{max}\) are all consistent with one another to within \(\rm \pm 1\sigma\). Therefore any of these guesses would suffice!

---