cooling age from the photometric fit

[sparx00]

Hello, first of all, this is an awesome repo!

I have been playing around the photometric fit of PSO J1801+6254 that is showcased in the repo. In ftr.best_fit_params there is an entry 'age' - does this correspond to the cooling age of the WD? If so, is there any method of obtaining its uncertainty or percentiles, as is possible for logg, Teff, or mass? Thanks.

[cylammarco]

Hi @sparx00, thank you for your interest in this repo.

The fitter defaults to use scipy's minimize, it is fast and good for simple cases. However, it does not come with an error estimation in its fit. Hence we get nan in the error entries.

>>> ftr.best_fit_params
{'H': {'chi2': 16.331071596034946, 'chi2dof': -1, 'Teff': 3945.5931387718992, 'Teff_err': nan, 'logg': 7.883606976283595, 'logg_err': nan, 'g_ps1': 16.69712684278966, 'distance': 71.231, 'dist_mod': 4.263345206871898, 'r_ps1': 15.704034694889183, 'i_ps1': 15.283477107938552, 'z_ps1': 15.095069184688736, 'y_ps1': 15.027153529904982, 'G3': 15.715138717830666, 'G3_BP': 16.412011674291037, 'G3_RP': 14.912256142772883, 'J_mko': 14.184129344073243, 'H_mko': 14.349943690969402, 'K_mko': 14.462301268674423, 'Av_g_ps1': 0.0, 'Av_r_ps1': 0.0, 'Av_i_ps1': 0.0, 'Av_z_ps1': 0.0, 'Av_y_ps1': 0.0, 'Av_G3': 0.0, 'Av_G3_BP': 0.0, 'Av_G3_RP': 0.0, 'Av_J_mko': 0.0, 'Av_H_mko': 0.0, 'Av_K_mko': 0.0, 'mass': 0.506789576678794, 'Mbol': 15.751993832322924, 'age': 8412712336.806402}, 'He': {'chi2': 75.8052680625422, 'chi2dof': -1, 'Teff': 4076.465537192801, 'Teff_err': nan, 'logg': 8.010280533194848, 'logg_err': nan, 'g_ps1': 16.651799655189226, 'distance': 71.231, 'dist_mod': 4.263345206871898, 'r_ps1': 15.865175902642834, 'i_ps1': 15.475428665678795, 'z_ps1': 15.298474684843848, 'y_ps1': 15.219156145386872, 'G3': 15.850655379868732, 'G3_BP': 16.450733760160897, 'G3_RP': 15.104976599249024, 'J_mko': 14.25711164011887, 'H_mko': 14.00191185593229, 'K_mko': 14.065524105238422, 'Av_g_ps1': 0.0, 'Av_r_ps1': 0.0, 'Av_i_ps1': 0.0, 'Av_z_ps1': 0.0, 'Av_y_ps1': 0.0, 'Av_G3': 0.0, 'Av_G3_BP': 0.0, 'Av_G3_RP': 0.0, 'Av_J_mko': 0.0, 'Av_H_mko': 0.0, 'Av_K_mko': 0.0, 'mass': 0.5744639170135237, 'Mbol': 15.793462337016509, 'age': 7699932412.531036}}

If we choose the least squares method by providing the keyword argument, e.g. fitter.fit(..., ..., ..., method="least_squares", ...). We get the uncertainties estimated.

>>> ftr.best_fit_params
{'H': {'chi2': 16.52680031325726, 'chi2dof': 9, 'Teff': 3934.219671470367, 'Teff_err': 28.648180206953782, 'logg': 7.881944093560716, 'logg_err': 0.02940730383374193, 'g_ps1': 16.70888624682886, 'distance': 71.231, 'dist_mod': 4.263345206871898, 'r_ps1': 15.711652542307078, 'i_ps1': 15.289897484103534, 'z_ps1': 15.100896653873047, 'y_ps1': 15.033682278021793, 'G3': 15.722769962685309, 'G3_BP': 16.421765886432173, 'G3_RP': 14.918583927468394, 'J_mko': 14.195885928621847, 'H_mko': 14.369728253162275, 'K_mko': 14.485824544732441, 'Av_g_ps1': 0.0, 'Av_r_ps1': 0.0, 'Av_i_ps1': 0.0, 'Av_z_ps1': 0.0, 'Av_y_ps1': 0.0, 'Av_G3': 0.0, 'Av_G3_BP': 0.0, 'Av_G3_RP': 0.0, 'Av_J_mko': 0.0, 'Av_H_mko': 0.0, 'Av_K_mko': 0.0, 'mass': 0.5058299512803053, 'Mbol': 15.762438457025754, 'age': 8427793337.953575}, 'He': {'chi2': 76.35827997891005, 'chi2dof': 9, 'Teff': 4105.24361245898, 'Teff_err': 12.182092418109297, 'logg': 8.046206645238266, 'logg_err': 0.009590963590389428, 'g_ps1': 16.653178582729726, 'distance': 71.231, 'dist_mod': 4.263345206871898, 'r_ps1': 15.877209171932067, 'i_ps1': 15.493990720076352, 'z_ps1': 15.32115313413484, 'y_ps1': 15.24413213187482, 'G3': 15.864169045133979, 'G3_BP': 16.456298830618113, 'G3_RP': 15.124501125674778, 'J_mko': 14.286884579041375, 'H_mko': 14.03113767716474, 'K_mko': 14.08401755028258, 'Av_g_ps1': 0.0, 'Av_r_ps1': 0.0, 'Av_i_ps1': 0.0, 'Av_z_ps1': 0.0, 'Av_y_ps1': 0.0, 'Av_G3': 0.0, 'Av_G3_BP': 0.0, 'Av_G3_RP': 0.0, 'Av_J_mko': 0.0, 'Av_H_mko': 0.0, 'Av_K_mko': 0.0, 'mass': 0.5971444628635761, 'Mbol': 15.81143239002258, 'age': 7860897887.069938}}

Hope this helps!

[sparx00]

Hi @cylammarco, thanks for the response. I can see that in the output you posted there is an estimate for Teff_err and logg_err, but not for the WD age. I have also tried it on my machine with the same results. If I swap out 'logg' for the 'mass' in the 'independent' at the fit input I can get an error estimate for the mass using emcee (although I cannot make this work for the least squares method). Is there any way to get an error estimate on the age, or are the error estimates limited to Mbol, Teff, logg and mass? Thank you.

[cylammarco]

It is not currently possible to get an uncertainty in the cooling age. I'll think about if this can be implemented (easily).

I'll have a look at the case when "logg" is swapped to "mass".

[cylammarco]

(Regarding age, please see also the next thread.)

[sparx00]

Thanks. Is the 'default' model used for the ftr.fit() from Bedard et al. 2020?

[cylammarco]

The age is indeed the cooling age (as the age is 0 when the DA temperature is > 110000).

[cylammarco]

@sparx00 sorry I deleted your message, I meant to delete my reply on the wrong thread... I believe you asked if it is possible to get the cooling age from a different model. Currently, it is not possible to directly obtain that from the .fit(). However, it is possible to construct a cooling interpolator with a CoolingModelReader by to get a cooling age by providing the luminosity (at 1 sigma limit) and mass (at 1 sigma limit) by doing:

  from WDPhotTools.cooling_model_reader import CoolingModelReader
  
  cr = CoolingModelReader()
  
  # For example, choosing the LPCODE ONe DA model
  cr.get_cooling_model('lpcode_one_da_19')
  
  # Interpolate
  cr.compute_cooling_age_interpolator()
  
  # get the cooling age by provide (log10(L), M solar mass)
  # e.g. 10^30 and 1.0 solar mass
  cr.cooling_interpolator(30, 1.0)
  
  >>> array(3.31069092e+09)

I will look into making this a simpler procedure in the coming weeks.

[sparx00]

Hello. Thank you, this helps.

[sparx00]

@sparx00 sorry I deleted your message, I meant to delete my reply on the wrong thread... I believe you asked if it is possible to get the cooling age from a different model. Currently, it is not possible to directly obtain that from the .fit(). However, it is possible to construct a cooling interpolator with a CoolingModelReader by to get a cooling age by providing the luminosity (at 1 sigma limit) and mass (at 1 sigma limit) by doing:

  from WDPhotTools.cooling_model_reader import CoolingModelReader
  
  cr = CoolingModelReader()
  
  # For example, choosing the LPCODE ONe DA model
  cr.get_cooling_model('lpcode_one_da_19')
  
  # Interpolate
  cr.compute_cooling_age_interpolator()
  
  # get the cooling age by provide (log10(L), M solar mass)
  # e.g. 10^30 and 1.0 solar mass
  cr.cooling_interpolator(30, 1.0)
  
  >>> array(3.31069092e+09)

I will look into making this a simpler procedure in the coming weeks.

Hi, @cylammarco. There might be a bug somewhere in the interpolator. I have just noticed that the code example with (L=10^30 W; M_wd = 1.0 M_sol) gives unphysical results. That WD seems too bright. For instance, Sun has a luminosity around 4x10^26 W and WDs are typically a few orders fainter, especially the heavier ones.

Also, if you swap the cooling models used by the interpolator, the result stays the same. E.g., changing 'lpcode_one_da_19' to 'montreal_co_da_20' gives exactly the same output.

[cylammarco]

Hi, I think there may be some cgs units involved. I'll dig into it and see if I can find any documentation detailing the units for the various models.

[sparx00]

Hi, I did not take that into the consideration. If the luminosity is in ergs/s then it does not seem that strange to me.