StellaR: a Package to Manage Stellar Evolution Tracks and Isochrones

We present the R package stellaR, which is designed to access and manipulate publicly available stellar evolutionary tracks and isochrones from the Pisa low-mass database. The procedures of the extraction of important stages in the evolution of a star from the database, of the isochrones construction from stellar tracks and of the interpolation among tracks are discussed and demonstrated.

lar clusters field. The theoretical effort of recovering the evolutionary history of the cluster benefits by the computation of extensive databases of stellar tracks and isochrones, such as Pietrinferni et al. (2006); Dotter et al. (2008); Bertelli et al. (2008). We recently computed a large data set of stellar tracks and isochrones, "The Pisa low-mass database" , with up to date physical and chemical inputs, and made available all the calculations to the astrophysical community at the Centre de Données astronomiques de Strasbourg (CDS) 1 , a data center dedicated to the collection and worldwide distribution of astronomical data.
In most databases, the management of the information and the extraction of the relevant evolutionary properties from libraries of tracks and/or isochrones is the responsibility of the end users. Due to its extensive capabilities of data manipulation and analysis, however, R is an ideal choice for these tasks. Nevertheless R is not yet well known in astrophysics; up to December 2012 only seven astronomical or astrophysical-oriented packages have been published on CRAN (see the taskview Chemometrics and Computational Physics).
The package stellaR  is an effort to make available to the astrophysical community a basic tool-set with the following capabilities: retrieve the required calculations from CDS; plot the information in a suitable form; construct by interpolation tracks or isochrones of compositions different to the ones available in the database; construct isochrones for age not included in the database; extract relevant evolutionary points from tracks or isochrones.

Get stellar evolutionary data
The Pisa low-mass database contains computations classified according to four parameters: the metallicity z of the star, its initial helium value y, the value of α-enhancement of the heavy elements mixture with respect to the reference mixture and the mixing length parameter α ml used to model external convection efficiency. The values of the parameters available in the database can be displayed using the function showComposition(): Upon specification of the aforementioned parameters, the stellaR package can import data from CDS (via anonymous ftp) from an active Internet connection. The CDS data are stored in ASCII format, and include a header with calculation metadata, such as the metallicity, the inital helium abundance, the mixing-length. The import is done via a read.table call, skipping the header of the files.
The following data objects can be downloaded from the database site: • Stellar track: a stellar evolutionary track computed starting from Pre-Main Sequence (PMS) and ending at the onset of helium flash (for masses M ≥ 0.55 M ) or at the exhaustion of central hydrogen (for 0.30 M ≤ M ≤ 0.50 M ). The functions getTrk() and getTrkSet() can be used to access such data; they respectively return objects of classes "trk" and "trkset".
• Stellar ZAHB: Zero-Age Horizontal-Branch models. The functions getZahb() can be used to access such data; it returns an object of class "zahb".
• HB models: computed from ZAHB to the onset of thermal pulses. The functions getHb() and getHbgrid() can be used to access such data; they respectively return objects of classes "hb" and "hbgrid".
• Stellar isochrones: computed in the age range [8 -15] Gyr. The functions getIso() and getIsoSet() can be used to access such data; they respectively return objects of classes "iso" and "isoset".
Readers interested in details about the computation procedure are referred to . The data gathered from CDS are organized into objects of appropriate classes. The package includes print and plot S3 methods for the classes "trk", "trkset", "hb", "zahb", "hbgrid", "iso", and "isoset".
As an example, we illustrate the recovering of the stellar track for a model of mass M = 0.80 M , metallicity z = 0.001, initial helium abundance y = 0.25, mixing-length α ml = [1] "mass" "z" "y" "ml" "alpha.enh" "data" > class(track) [1] "trk" "stellar" The function getTrk() returns an object of class "trk", which is a list containing the track metadata, i.e. as the star mass, the metallicity, the initial He abundance, the mixing-length and the αenhancement, and the computed data in the data frame data. Track data contains the values of 15 variables: > names(track$data) [1] "mod" "time" "logL" "logTe" "mass" "Hc" "logTc" "logRHOc" "MHEc" "Lpp" "LCNO" "L3a" "Lg" "radius" "logg" The included variables are: mod the progressive model number; time the logarithm of the stellar age (in yr); logL the logarithm of the surface luminosity (in unit of solar luminosity); logTe the logarithm of the effective temperature (in K); mass the stellar mass (in unit of solar mass); Hc the central hydrogen abundance ( [1] "mass" "massRGB" "z" "y" "ml" "alpha.enh" "data" > class(hbtk) [1] "hb" "stellar" which returns an object of class "hb", which differs from an object of class "trk" only for the presence of the variable massRGB, i.e. the Red-Giant Branch (RGB) progenitor mass. Usually a set of tracks with different mass and/or metallicity are needed for computations. The package stellaR provides the function getTrkSet(), which can download a set of tracks with different values of mass, metallicity, helium and mixinglength. As an example the whole set of masses (from 0.30 to 1.10 M , step of 0.05 M ), for metallicity z = 0.001, initial helium abundance y = 0.25, mixinglength α ml = 1.90, α-enhancement [α/Fe] = 0.0 can be downloaded as follows:  The function getTrkSet() returns an object of class "trkset", a list containing objects of class "trk". The track set can be displayed in the usual (log T eff , log L/L ) plane by a call of the function plot(): > plot(trks, lty = 1:2) The output of the function is shown in Figure 1. The plot is produced by a call to the function plotAstro() which allows the user to customize several aspects of the plot, such as the axes label, the number of minor ticks between two major ticks, the limits of the axes, the color and type of the lines (as in the example), the type of the plot (lines, points, both, . . . ).
The use of plot() function is further demonstrated in Figure 2 Apart from the plots discussed before, it is easy to display other relation between the computed variables. In the following example we get the data for two masses, namely 0.50 and 1.00 M and plot the trend of the radius (in units of solar radius) versus the logarithm of the age for the first 100 models. The resulting plot is displayed in Fig. 3.

Tools for interpolating among data structures
Even if the database provides a fine grid of models it is possible that the specific model needed for a comparison with observational data has not been calculated. To address this issue, the package stellaR includes tools for construction of new sets of isochrones. The simplest case is whenever one desires an isochrone for ages not included in the database, but for a combination of metallicity, He abundance, mixing-length and α-enhanchment existing in the database. The function makeIso() is able to compute the required isochrones by mean of interpolation on a tracks set. The interpolation can be performed for age in the range [7 -15] Gyr. The user has the choice to explicitly provide to the function a set of tracks, previously downloaded from CDS, or to specify the required composition of the tracks to be downloaded for the interpolation. To show the usage of the function we use the object trks downloaded before to obtain an isochrone of age 9.7 Gyr: > iso.ip <-makeIso(age = 9.7, tr = trks) > iso.ip Stellar isochrone Age = 9.7 Gyr Z = 0.001 , Y = 0.25 Mixing length = 1.9 [alpha/Fe] = 0 the call produces a result identical to makeIso(age = 9.7,z = 0.001,y = 0.25,ml = 1.9,afe = 0); in this last case the data are taken from CDS before the interpolation procedure. The interpolation technique is based upon the fact that all the tracks contain the same number of points by construction, and that a given point corresponds to the same evolutionary phase on all the tracks. We define S(M) as the set of tracks to be used for interpolation, parametrized by the value of the mass M. Let t i (M) be the evolutionary time for the ith point on the track of mass M, and A be the age of the required isochrone. Let k be the point on the track of lower mass of the set S(M) for which t k (M) ≥ A. For each point j ≥ k on the tracks in S(M), a linear interpolation in age of the values of mass, logarithm of the effective temperature and logarithm of the luminosity is performed among tracks. These points define the required isochrone. A potential problem of this simple procedure will occur whenever massive stars develop a convective core during the Main Sequence (MS). In this case, as shown for example in Mowlavi et al. (2012), the monotonic trend of the evolutionary time -that decreases with increasing stellar mass at the end of the MS -inverts at the middle of the MS. However the problem will be encountered only for early-time isochrones, for which the mass at the isochrone Turn-Off will be in the interval during which the convective core develops. The procedure outlined in this Section is adequate for construction of isochrones throughout the range allowed by the function.

Tracks interpolation
The package stellaR provides also a tool for performing a 3D interpolation on the database to construct a set of tracks for values of metallicity, helium abundance and mixing-length not included in the computations available at CDS. The function interpTrk() can be used for this procedure. A call to this function causes the download from CDS of the sets of tracks needed for the interpolation.
The new set of tracks is computed by mean of a linear interpolation. The metallicity is logtransformed before the interpolation procedure. Let T z,y,α ml (M) i be ith point on the dataset containing the evolutionary time, the effective temperature and the logarithm of surface luminosity for the track of mass M, and given composition. The interpolation algo-rithm proceeds as follows: where the symbol * means that interpolation occurred in the substituted variable. The selection of the tracks set which enter in the interpolation is based upon the identification of the vertexes of the cell of the (z, y, α ml ) space containing the point identified by the required parameters. Then, for all the 17 masses at the vertexes, the linear interpolation described above is performed. In the worst case scenario, whenever no one of the supplied parameters values exist in the database, the interpolation requires 2 3 = 8 sets of 17 masses. The algorithm is however able to reduce the dimensionality of the process if some of the variables values exist in the database.
As a demonstration, let us compute a set of tracks with mixing-length value α ml = 1.74, z = 0.002, y = 0.25, [α/Fe] = 0.0: > ip.trk <-interpTrk(z = 0.002, y = 0.25, ml = 1.74, afe = 0) Since the values of z and y exist in the database, only an interpolation on the mixing-length value is performed by the code. The set of tracks can be used for isochrones construction, like a standard set of tracks: > ip.iso <-makeIso(age = 12, tr = ip.trk)

Keypoints extraction
Important stages in the evolution of a star are defined as "keypoints", e.g. hydrogen core exhuastion, Turn-Off luminosity, RGB tip luminosity. To simplify their extraction the package stellaR provides the function keypoints(), which operates on an object of class "trk" or "iso". The function extracts from the data stored in objects of class "trk" the rows of the data frame relative to the following evolutionary stages: When the function is called on an object of class "iso" it returns a data frame containing only TO and BTO phases. In both cases the function inserts in the returned data frame the columns relative to mass (or age for isochrones), metallicity, initial He value, mixing-length, α-enhancement, and evolutionary phase identifier.
As a demonstration we extract the TO and BTO points from the object isoc generated in a previous example: The points can be easily superimposed to the isochrones in the theoretical plane. The top panel of Figure 5 is obtained with the following commands: > plot(isoc) > points(logL~logTe, data = kp, pch = id, col = id) > legend("topleft", c("TO", "BTO"), pch = 1:2, col = 1:2) As a last example we extract a set of tracks for masses in the range [0.40 -1.10] M and three metallicity z = 0.0001, 0.001, 0.01 and we display (bottom panel in Figure 5) the time of exhaustion of central hydrogen as a function of the mass of the star.

Summary
This paper demonstrated how the package stellaR can be useful to the astrophysical community as a tool to simplify the access to stellar tracks and isochrones calculations available on-line. A set of tools to access, manipulate and plot data are included in the package and their usage was shown. The interpolation functions included in the package can be used to safely produce tracks or isochrones at compositions not included in the database, without the need that users develop software on their own. A planned extension of the package is the modification of the algorithm of isochrone construction to make feasible the calculation of isochrone of young ages. This step can be useful in view of a possible extension of the Pisa database to higher masses, or to manipulation of data stored in other databases. In fact, while the package is currently developed for accessing data from Pisa low-mass database, other public databases can be in principle accessed in the same way. This step is however complicated by the fact that no standard for stellar model output exists in the astrophysical community, requiring the personalization of the interface functions for each set of data.