Stellar model grids

Background and EEPs

Stellar model grids are typically constructed as a set of evolutionary tracks, where models of stellar evolution are run on grids of initial mass and metallicity, often with some other physical parameter varied as well (e.g., rotation, helium fraction, \(\alpha\)-abundance, etc.). Each of these evolutionary tracks predicts various physical properties (temperature, luminosity, etc.) of a star with given initial mass and metallicity, as a function of age.

It is also often of interest to re-organize these evolution track grids into “isochrones”—sets of stars at a range of masses, all with the same age. As described in this reference, in order to construct these isochrones, the time axis of each evolution track gets mapped into a new coordinate, called “equivalent evolutionary phase,” or EEP. The principle of the EEPs is to first identify physically significant stages in stellar evolution, and then subdivide each of these stages into a number of equal steps. This adaptive sampling enables accurate interpolation between evolution tracks even at ages when stars are evolving quickly, in the post-main sequence phases.

Previous versions of isochrones relied directly on these precomputed isochrone grids and interpolated between grid points in (mass, age, feh) space. This returned inaccurate results for post-MS stages of stellar evolution, and thus was not reliable for modeling evolved stars. However, beginning with v2.0, isochrones now implements all interpolation using EEPs. In addition, it provides direct access to the evolution track grids, in addition to precomputed isochrone grids. Note that version 2.0 includes only the MIST models; future updates will include more (e.g. PARSEC, YAPSI).

Model Grid Objects and Interpolation

Isochrones provides a simple and direct interface to full grids of stellar models. Upon first access, the grids are downloaded in original form, reorganized, and written to disk in binary format in order to load quickly with subsequent access. The grids are loaded as pandas dataframes with multi-level indexing that reflects the structure of the grids: evolution track grids are indexed by metallicity, initial mass, and EEP; and isochone grids by metallicity, age, and EEP.

from isochrones.mist import MISTEvolutionTrackGrid, MISTIsochroneGrid

track_grid = MISTEvolutionTrackGrid()
track_grid.df.head()  # just show first few rows

			nu_max	logg	eep	initial_mass	radius	logTeff	mass	density	Mbol	phase	feh	Teff	logL	delta_nu	interpolated	star_age	age	dt_deep
initial_feh	initial_mass	EEP
-4.0	0.1	1	143.524548	3.033277	1.0	0.1	1.593804	3.620834	0.1	0.034823	5.132871	-1.0	-3.978406	4176.707371	-0.157148	21.776686	False	13343.289397	4.125263	0.026168
		2	145.419039	3.038935	2.0	0.1	1.583455	3.620769	0.1	0.035510	5.147664	-1.0	-3.978406	4176.085183	-0.163066	21.993078	False	14171.978264	4.151430	0.026121
		3	147.409881	3.044805	3.0	0.1	1.572790	3.620702	0.1	0.036237	5.163015	-1.0	-3.978406	4175.435381	-0.169206	22.219791	False	15048.910447	4.177505	0.026016
		4	149.499346	3.050886	4.0	0.1	1.561817	3.620631	0.1	0.037006	5.178922	-1.0	-3.978406	4174.757681	-0.175569	22.457004	False	15975.827275	4.203463	0.025996
		5	151.703570	3.057203	5.0	0.1	1.550499	3.620558	0.1	0.037823	5.195452	-1.0	-3.978406	4174.049081	-0.182181	22.706349	False	16962.744747	4.229496	0.025996

iso_grid = MISTIsochroneGrid()
iso_grid.df.head()  # just show first few rows

			eep	age	feh	mass	initial_mass	radius	density	logTeff	Teff	logg	logL	Mbol	delta_nu	nu_max	phase	dm_deep
log10_isochrone_age_yr	feh	EEP
5.0	-4.0	35	35	5.0	-3.978406	0.100000	0.100000	1.106082	0.104184	3.617011	4140.105252	3.350571	-0.489734	5.964335	37.987066	299.346079	-1.0	0.002885
		36	36	5.0	-3.978406	0.102885	0.102885	1.122675	0.102507	3.618039	4149.909661	3.347798	-0.472691	5.921728	37.739176	298.570836	-1.0	0.003573
		37	37	5.0	-3.978406	0.107147	0.107147	1.147702	0.099921	3.619556	4164.436984	3.343658	-0.447471	5.858678	37.345115	297.180748	-1.0	0.004247
		38	38	5.0	-3.978406	0.111379	0.111379	1.173015	0.097287	3.621062	4178.903372	3.339612	-0.422498	5.796244	36.923615	295.526946	-1.0	0.004217
		39	39	5.0	-3.978406	0.115581	0.115581	1.198615	0.094627	3.622555	4193.289262	3.335660	-0.397776	5.734440	36.473151	293.589960	-1.0	0.004189

This generally contains only a subset of the original columns provided by the underlying grid, with standardized names. There are also additional computed columns, such as stellar radius and density. The full, original grids, can be found with the .df_orig attribute if desired:

iso_grid.df_orig.head()  # just show first few rows

			EEP	log10_isochrone_age_yr	initial_mass	star_mass	star_mdot	he_core_mass	c_core_mass	o_core_mass	log_L	log_L_div_Ledd	...	nu_max	acoustic_cutoff	max_conv_vel_div_csound	max_gradT_div_grada	gradT_excess_alpha	min_Pgas_div_P	max_L_rad_div_Ledd	e_thermal	phase	feh
log10_isochrone_age_yr	feh	EEP
5.0	-4.0	35	35	5.0	0.100000	0.100000	-1.455094e-13	0.0	0.0	0.0	-0.489734	-4.167035	...	299.346079	2233.536029	0.127243	1.095544	0.0	0.999989	0.000016	3.002314e+46	-1.0	-4.0
		36	36	5.0	0.102885	0.102885	-1.562027e-13	0.0	0.0	0.0	-0.472691	-4.129623	...	298.570836	2228.014832	0.128938	1.101114	0.0	0.999989	0.000017	3.106838e+46	-1.0	-4.0
		37	37	5.0	0.107147	0.107147	-1.707298e-13	0.0	0.0	0.0	-0.447471	-4.073591	...	297.180748	2218.440338	0.130528	1.109114	0.0	0.999988	0.000019	3.264230e+46	-1.0	-4.0
		38	38	5.0	0.111379	0.111379	-1.836256e-13	0.0	0.0	0.0	-0.422498	-4.017245	...	295.526946	2207.403678	0.132657	1.116760	0.0	0.999987	0.000021	3.424460e+46	-1.0	-4.0
		39	39	5.0	0.115581	0.115581	-1.949639e-13	0.0	0.0	0.0	-0.397776	-3.960633	...	293.589960	2194.776391	0.134294	1.124050	0.0	0.999986	0.000022	3.587613e+46	-1.0	-4.0

5 rows × 80 columns

iso_grid.df_orig.columns

Index(['EEP', 'log10_isochrone_age_yr', 'initial_mass', 'star_mass',
       'star_mdot', 'he_core_mass', 'c_core_mass', 'o_core_mass', 'log_L',
       'log_L_div_Ledd', 'log_LH', 'log_LHe', 'log_LZ', 'log_Teff',
       'log_abs_Lgrav', 'log_R', 'log_g', 'log_surf_z', 'surf_avg_omega',
       'surf_avg_v_rot', 'surf_num_c12_div_num_o16', 'v_wind_Km_per_s',
       'surf_avg_omega_crit', 'surf_avg_omega_div_omega_crit',
       'surf_avg_v_crit', 'surf_avg_v_div_v_crit', 'surf_avg_Lrad_div_Ledd',
       'v_div_csound_surf', 'surface_h1', 'surface_he3', 'surface_he4',
       'surface_li7', 'surface_be9', 'surface_b11', 'surface_c12',
       'surface_c13', 'surface_n14', 'surface_o16', 'surface_f19',
       'surface_ne20', 'surface_na23', 'surface_mg24', 'surface_si28',
       'surface_s32', 'surface_ca40', 'surface_ti48', 'surface_fe56',
       'log_center_T', 'log_center_Rho', 'center_degeneracy', 'center_omega',
       'center_gamma', 'mass_conv_core', 'center_h1', 'center_he4',
       'center_c12', 'center_n14', 'center_o16', 'center_ne20', 'center_mg24',
       'center_si28', 'pp', 'cno', 'tri_alfa', 'burn_c', 'burn_n', 'burn_o',
       'c12_c12', 'delta_nu', 'delta_Pg', 'nu_max', 'acoustic_cutoff',
       'max_conv_vel_div_csound', 'max_gradT_div_grada', 'gradT_excess_alpha',
       'min_Pgas_div_P', 'max_L_rad_div_Ledd', 'e_thermal', 'phase', 'feh'],
      dtype='object')

Any property (or properties) of these grids can be interpolated to any value of the index parameters via the .interp method:

track_grid.interp([-0.12, 1.01, 353.1], ['mass', 'radius', 'logg', 'Teff'])

array([1.00983180e+00, 1.04691913e+00, 4.40266419e+00, 6.03383320e+03])

Similarly, the .interp_orig method interpolates any of the original columns by name:

track_grid.interp_orig([-0.12, 1.01, 353.1], ['v_wind_Km_per_s'])

array([2.87408918e-05])

Note that these interpolations are fast—30-40x faster than the equivalent interpolation in scipy, for evaluating at a single point:

from scipy.interpolate import RegularGridInterpolator

grid = track_grid.interp.grid[:, :, :, 4]  # subgrid corresponding to radius
interp = RegularGridInterpolator(track_grid.interp.index_columns, grid)
assert track_grid.interp([-0.12, 1.01, 353.1], ['radius']) == interp([-0.12, 1.01, 353.1])

%timeit interp([-0.12, 1.01, 353.1])
%timeit track_grid.interp([-0.12, 1.01, 353.1], ['radius'])

The slowest run took 14.68 times longer than the fastest. This could mean that an intermediate result is being cached.
100 loops, best of 3: 1.13 ms per loop
The slowest run took 5.04 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 12.5 µs per loop

In order to select a subset of these grids, you can use pandas multi-index magic:

iso_grid.df.xs((9.0, 0.0), level=(0, 1)).head()  # just show first few rows

	eep	age	feh	mass	initial_mass	radius	density	logTeff	Teff	logg	logL	Mbol	delta_nu	nu_max	phase	dm_deep
EEP
193	193	9.0	0.042799	0.100000	0.100000	0.126216	70.115891	3.460248	2885.680821	5.235913	-3.002129	12.245322	1045.120425	27533.135930	-1.0	0.003449
194	194	9.0	0.042805	0.103449	0.103449	0.129201	67.622476	3.462649	2901.679870	5.227759	-2.972222	12.170556	1025.261668	27025.521973	-1.0	0.004051
195	195	9.0	0.042814	0.108103	0.108103	0.133343	64.282466	3.465890	2923.411541	5.216743	-2.931855	12.069637	998.429682	26339.315728	-1.0	0.004742
196	196	9.0	0.042824	0.112932	0.112932	0.137791	60.858291	3.469256	2946.160133	5.205249	-2.889887	11.964717	970.463953	25622.963627	-1.0	0.004720
197	197	9.0	0.042834	0.117543	0.117544	0.142177	57.660112	3.472471	2968.048014	5.194272	-2.849811	11.864529	943.753111	24938.691941	-1.0	0.004735

Example visualization

Just for fun, let’s plot a few isochrones:

import hvplot.pandas

# Select two isochrones from the grid
iso_df1 = iso_grid.df.xs((9.0, 0.0), level=(0, 1))
iso_df2 = iso_grid.df.xs((9.5, 0.0), level=(0, 1))

options = dict(invert_xaxis=True, legend_position='bottom_left')

# Isn't hvplot/holoviews great?
plot1 = iso_df1.hvplot.line('logTeff', 'logL', label='Log(age) = 9.0')
plot2 = iso_df2.hvplot.line('logTeff', 'logL', label='Log(age) = 9.5')
(plot1 * plot2).options(**options)