Absolute parameters for AI Phoenicis using WASP photometry

Located at the Observatorio del Roque de los Muchachos, La Palma and at Sutherland Observatory, South Africa, the two WASP instruments are both formed from eight wide-field cameras each with a 2048 x 2048 pixel CCD. For AI Phe (1SWASP J010934.19-461556.0), over 170 000 photometric measurements were present in the WASP archive, taken between June 2006 and January 2014 by the instrument in South Africa. During this period, WASP-South has used two different types of lenses, the original 200-mm, f/1.8 lenses (Pollacco et al., 2006) and, from July 2012, 85-mm, f/1.2 lenses (Smith & WASP Consortium, 2014). The reduction procedure is identical for both lens types, and uses a dedicated pipeline (Pollacco et al., 2006) that has been optimised in each case. The data are then processed by a detrending algorithm, which was developed from the SysRem algorithm of Tamuz et al. (2005) and is described by Collier Cameron et al. (2006).

Use of the 85-mm lenses allows WASP-South to focus on brighter stars ($V\lesssim 9$). A larger reduction aperture is used in comparison to the 200-mm lenses (4 pixels from 3.5) and the limits on image and star rejection were also modified. This results in a shift of the photometry range to $6\lesssim V\lesssim 11$ for the 85-mm lenses from $9\lesssim V\lesssim 13$ for the 200-mm lenses. The 200-mm lenses use broad-band filters with a range of 400-700 nm, while the 85-mm use SDSS r’ filters. Observations from the two different types of lens have been analysed separately.

The WASP data can suffer from large amounts of scatter, due to clouds, etc. This section describes the methods used to remove unreliable observations in the two sets of data.

For the 200-mm data, a decision was made to only use data from cameras 225 and 226 during the analysis. These two cameras contribute more than 80% of the observations made with the 200-mm cameras. The data from the other cameras contained large amounts of scatter, perhaps related to the fact that the 200-mm lenses are very close to their saturation limit for AI Phe. They also only contributed a handful of observations to the eclipses, and would not help determine the radii of AI Phe.

As part of the WASP reduction pipeline, each photometric measurement gets assigned a weighting factor, $\sigma_{\rm xs}$, (denoted $\sigma_{t(i)}$ in Collier Cameron et al., 2006), to characterise the scatter present from external noise sources. This $\sigma_{\rm xs}$ value is set to zero if the pipeline deems the value to be missing or bad (Collier Cameron et al., 2006), so data with $\sigma_{\rm xs}=0$ were not included in any analysis. In some cases, the error in the measured flux was exceptionally high in comparison with the other measurements. It is likely that this is due to clouds being present at the time of exposure, as $\sigma_{\rm xs}$ is also higher for these data. If the error in a particular flux measurement was five or more times greater than the median error of all data from the same camera, then it was excluded from our analysis.

The flux measurement we have used for our analysis, $f$, is calculated from an aperture with a radius of 3.5 pixels or 4 pixels (for the 200-mm and 85-mm lenses, respectively) centred on AI Phe, with the detrending correction of Collier Cameron et al. (2006) applied, and $f_{\rm err}$ is the associated error in this value. The magnitude of each measurement was calculated from $f$ using the median flux of the data set as the zero point. The error in magnitude, $m_{\rm err}$, associated with each observation, was calculated using

One final check was used to look for sections of data that were significantly offset from the remaining data. The phase-folded data was split into 800 phase bins. The median magnitude and associated error in the median were then calculated for each bin. This formed a model from which we could estimate the expected magnitude of an observation given its phase. The entire data set was also split into blocks based on the night the observation was taken and the camera used. For each data block, the magnitude of each observation within that block was compared to its expected magnitude by matching it to a phase bin. If more than 80% of the observations from a particular data block were offset by more than ten times the errors associated with the appropriate phase-bins, then all data from that night/camera combination were excluded from the analysis.

Overall, from the initial 170 324 observations stored in the WASP archive, 126 780 remained for use during the analysis, 12 618 from the 200-mm lenses and 114 162 from the 85-mm lenses. Table 1 provides a summary of the number of points removed during this initial processing stage for each type of lens.

From an initial analysis of the data using the ephemeris given in Hrivnak & Milone (1984), we found that the primary eclipse in the 85-mm data was offset in phase by $0.00102\pm 0.00002$. This meant that the primary eclipse occurred more than 30 minutes later than predicted.

Table 2 details the times of primary minima, $t_{\rm pri}$, currently available for AI Phe. Two of the times have been previously published. The remaining five have been obtained by fitting data from WASP and All-Sky Automated Survey (ASAS, Pojmanski, 2002) using jktebop (Southworth et al., 2004). The 200-mm WASP data were split into three blocks – the range of Heliocentric Julian Date used in each block is given in Table 2. Only one block was used from the 85-mm data. Figure 1 shows the difference between observed $t_{\rm pri}$ and the calculated time based on the linear ephemeris of Hrivnak & Milone (1984). Based on this plot, a linear ephemeris cannot be used to describe the long-term periodicity of AI Phe. There are an insufficient number of minima to obtain a reliable quadratic ephemeris. Therefore, for this work, the following linear ephemeris has been fitted to only times of minima for the WASP data:

For the shorter timescale covered by the WASP data, a linear ephemeris is a suitable approximation. The new period is consistent with the value determined spectroscopically by Hełminiak et al. (2009), however, there is an 8.4-sigma difference between this new period and the value quoted in Hrivnak & Milone (1984). The zero-point of the ephemeris from Hełminiak et al. (2009) was not included in Fig. 1 and Table 2, because they use a different definition of this quantity. With the new ephemeris the phase-offset for the 85-mm data is reduced to 0.00005$\pm$0.00002.

In all the observations of AI Phe over the last 40 years or so, the secondary eclipse has not been observed in its entirety in one night, so it has not been possible to investigate the timings of the secondary eclipse. This investigation would have given an insight into the possible cause of the apparently quadratic nature of the ephemeris. If the deviations in calculated timings are common to both the primary and secondary minima, it may suggest a third body is involved, or if the deviations have opposite signs then it may suggest apsidal motion is involved.

Figure 2 shows the phase-folded lightcurves for both the 85-mm and 200-mm data using the ephemeris in equation 2, having been processed using the method described in Sect. 2.2.

An image of AI Phe from the Two Micron All Sky Survey (2MASS, Skrutskie et al., 2006) shows there is another star approximately 11″ to the east of AI Phe, as shown in Fig. 3. The close proximity of the star means it resides within the aperture used for the WASP photometry of AI Phe and so contributes to the measured flux. As such, it has been necessary to include third light as a parameter that is fitted during the lightcurve analysis.

Although seen in the 2MASS images, this survey does not present any photometric measurements for this companion star. Near-infrared photometry for this star was taken from the Deep Near-Infrared Survey (DENIS, Epchtein et al., 1997), with $I=13.50\pm 0.02$ mag, $J=12.95\pm 0.07$ mag and $K=12.60\pm 0.14$ mag. The $I-J$, $J-K$, and $I-K$ colours for the companion were compared to the Dartmouth stellar evolution models (Dotter et al., 2008), and by assuming a main-sequence star with $[{\rm Fe}/{\rm H}]=0.0$ and an age 2 Gyr, its mass was estimated as $0.91\pm 0.06M_{\sun}$. This was used to obtain expected absolute magnitudes from the Dartmouth model, which were in turn used to calculate the distance modulus in each band, resulting in an estimated distance of $590\pm 9$pc to the companion. Comparing this value with the distance to AI Phe of $162\pm 6\,{\rm pc}$ (Andersen et al., 1988), we conclude that the visual companion is not physically associated with the system.

This work uses the ebop lightcurve analysis code (Nelson & Davis, 1972; Popper & Etzel, 1981). The subroutine light from this code, which calculates the lightcurves at specified orbital phases for a given set of model parameters, was converted to double-precision floating point arithmetic and modifications were made to enable it to be called directly from the programming language, Python. The ebop lightcurve is then combined with the least-squares Levenberg-Marquardt, Python module, MPFIT (Markwardt, 2009) to find the best-fit parameters, i.e. the parameters that minimise the $\chi^{2}$ value between the model and data.

The following seven parameters were allowed to vary during the fitting process: surface brightness ratio at the centre of the stellar discs, $J$; sum of the radii, $r_{\rm sum}=r_{1}+r_{2}$; ratio of the radii, $k=r_{2}/r_{1}$; inclination, $i$; $e\cos{\omega}$, $e\sin{\omega}$; third-light, $l_{3}$. In addition to these fitted parameters, the fractional radii, $r_{1}$ and $r_{2}$ are automatically calculated from $r_{\rm sum}$ and $k$, while the eccentricity $e$ and longitude of periastron $\omega$ are calculated from $e\cos{\omega}$ and $e\sin{\omega}$.

The mass ratio, $q=M_{2}/M_{1}$, was fixed at 1.0418, taken from Hełminiak et al. (2009). This is slightly larger than the value 1.034 quoted by Andersen et al. (1988), but modifying the mass ratio by such a small amount had no effect on the best-fit parameters. The gravity darkening exponents also have little impact on the shape of the lightcurve in the case of AI Phe as the system is well-detached. Using values taken from the table by Claret & Bloemen (2011), the gravity darkening exponent for the primary, $y_{\rm p}$, and secondary, $y_{\rm s}$, were fixed at 0.26 and 0.50 respectively. The same values were used for both the 85-mm and 200-mm data.

The linear limb-darkening coefficients were also held fixed. Attempts were made to include these as free parameters, but we found that they are not usefully constrained by the data. Instead, the limb darkening coefficients were estimated by interpolation in the table of Claret & Bloemen (2011) by using values for effective temperature, surface gravity and metallicity for the two components from Andersen et al. (1988). For the 85-mm data we adopted the primary limb darkening coefficient, $u_{\rm p}=0.54\pm 0.03$ and $u_{\rm s}=0.67\pm 0.03$ for the secondary. For the 200-mm data we used $u_{\rm p}=0.52\pm 0.05$ and $u_{\rm s}=0.67\pm 0.05$ using the Kepler pass-band to approximate the response of the WASP broad-band filter. To account for the uncertainties in the limb darkening coefficients, each coefficient was varied by its error and a new model fitted. For each parameter, the average absolute difference between these models and the best fit model has been added in quadrature to the uncertainties from the best-fit parameters. Table 3 details the typical contribution to the uncertainty of each parameter, for both the 200-mm and 85-mm data. These were calculated from the mean contributions for each best-fit detailed in Sects. 3.6 and 3.7.

One of the greatest risks when using a least squares minimisation method is that the solution that is found is a local minimum, rather than the overall global minimum of the problem. If local minima are present it is important that these are considered when the uncertainties for the final parameters are calculated. Markov Chain Monte Carlo (MCMC), in the form of the Python module, emcee, (Goodman & Weare, 2010) was used to check for these local minima, within the parameter-space of the same seven parameters used in the model fitting ($J$, $r_{\rm sum}$, $k$, $i$, $e\cos{\omega}$, $e\sin{\omega}$ and $l_{3}$). emcee uses the affine-invariant ensemble sampling (stretch-move) algorithm developed by Goodman & Weare (2010), where a group of walkers explore the parameter space. This group of walkers can be split allowing the process to be run in parallel and the affine-invariant transformations mean the algorithm can cope with skewed probability distributions. The positions of the walkers within a particular sub-group are updated using the positions of walkers in the other subgroups.

The probability that a model produced by the parameters from the walkers, corresponds to the best-fit model is evaluated using the log likelihood function

where y is a vector of length $N$ containing the magnitudes generated for a model, $\Theta$ is a vector containing the varying parameters ($J$, $r_{\rm sum}$, $k$, $i$, $e\cos{\omega}$, $e\sin{\omega}$ and $l_{3}$), $m$ is the observed magnitude and $m_{\rm err}$ is the standard error on the magnitude. Priors were applied, but these are only used to prevent the parameter exploring areas that are unphysical, eg. $r_{\rm sum}$ or $J$ being less than zero.

The process ran using 150 walkers for 2 500 steps, of which the first 200 were discarded to allow for an adequate burn-in stage. For each of the walkers, a starting point for each parameter was chosen by choosing a number at random from a normal distribution (with a mean of zero and variance of 0.01) and adding it to the best-fit parameter. While this method would initially create a ball of walkers close to the solution found by the model, the burn-in stage allows the walker to spread out from this ball. Each parameter was subsequently checked, by plotting the walkers’ positions against step number, to ensure the burn-in stage was completed within these first 200 steps. To ensure the was no bias from the walkers’ starting positions, a test was run that allowed the walkers to start further from the model solution, up to three times the parameters’ uncertainties as determined from the covariance matrix of the best-fitting model. The burn-in stage for this test was much longer, but also failed to reveal any other local minima, so we concluded the original starting points were adequate.

Figure 4 shows example distributions for the 85-mm data. The contours plotted over the top of the distributions indicate the density of points with the darker regions showing higher densities towards the centre of each plot. The grey lines across each of the distributions show where the best-fit values lie. Generally the resulting distributions are symmetric, as shown by the histograms for each parameter. However, there are correlations between some of the parameters, as seen by the tilted ellipsoid shaped regions. Many of these correlations link the third-light parameter to the other parameters, more specifically there are strong correlations between $l_{3}$ and $J$, $k$ and $i$. The correlation with $k$ highlights the importance for the inclusion of the third-light parameter in the model, because without it, the resulting relative radii would be subject to a systematic error.

It is inappropriate to use the covariance matrix from the least-squares fit to estimate the standard errors on the model parameters. This is because the assumption of uncorrelated noise with a Gaussian distribution is not satisfied for WASP photometry. It is not clear that the distribution of points from the MCMC chain gives an accurate impression of the posterior probability distribution of the parameters, for the same reason. Therefore, it is not appropriate to use the likelihood calculated using equation 3 either. As such, the standard errors on the best-fit parameters have been calculated using a prayer-bead (residual permutation) as described by Southworth (2008) and based on an algorithm developed by Jenkins et al. (2002). The residuals between the data and the best-fit model are shifted by a number of steps to create a synthetic data set. A new model is then fitted to this synthetic data, and the process is repeated across the entirety of the original data set. The uncertainties are the standard deviation of the fitted parameters from all the synthetic data models. Ideally, the number of shifts would be $N-1$, where $N$ is the number of observations used in generating the best-fit model. However, due to the large number of points involved in the WASP lightcurves, the number of shifts used was restricted to 500 spread evenly across the data. To ensure the best-fit solution was not affected by the choice of initial parameters, the initial parameters were taken at random from positions within the MCMC analysis.

As mentioned in Sec. 2.1, WASP photometry is processed using the detrending algorithm described in Collier Cameron et al. (2006). The algorithm is used to remove four trends of systematic errors which are common to all stars in a particular field and is given by the equation

where $m_{i,j}$ and $\widetilde{m}_{i,j}$ are the observed and corrected magnitude, respectively, for star $j$ at time $i$. M is the total number of trends, $a_{i}$ are basis functions detailing the patterns of systematic errors and $c_{j}$ describes to what extent each basis function affects a particular star. The algorithm is applied separately to each unique combination of camera and season, and aids the process locating planetary transits by flattening lightcurves over the time scale of a typical transit length (2.5 hours). The work in this paper aims to determine the radii of AI Phe to the highest accuracy possible using the WASP photometry, and so the effects of the detrending process on the resulting parameters have been investigated.

Reversing the detrending entirely would re-introduce all the systematic errors the algorithm was designed to remove, making comparisons between the resulting parameters more difficult. Instead, effective detrending coefficients, $c{{}^{\prime}}$, have been calculated. These coefficients take into consideration the variability of an eclipsing binary, by including a lightcurve model, $L$, when they are calculated. For AI Phe, the situation can be described using

where $m_{i}$ and $\widetilde{m}_{i}$ are the observed and corrected magnitude (respectively) for AI Phe, and ${{}^{(k)}}a_{i}$ are the same detrending basis functions as before. The effective detrending coefficients for each of the basis functions were calculated using singular value decomposition (SVD). Effects from the altered detrending coefficients were removed from the observed data, then using ebop and MPFIT, best-fit parameters were obtained. A best-fit model generated using the method in Sec. 3.3, was used for $L$ initially. However, to ensure the choice of initial model did not bias the results, the values of $c{{}^{\prime}}$ were calculated for multiple models. This was done as an iterative process, which continued until all parameters change by less than 0.005% with each new model. Once the final set of best-fit parameters was determined, MCMC and prayer-bead analysis were used to calculate their associated uncertainties.

Table 4 contains the best-fit parameters for both the 85-mm and 200-mm, with the difference between the original data and detrended data included for comparison. With the exceptions of $J$ and $l_{3}$ for the 200-mm data and $e\cos\omega$ for the 85-mm, all the best-fit parameters for the detrended data lie within the uncertainties of the parameters from the original data. Therefore, in general, the detrending algorithm applied during the WASP pipeline does not affect the overall best-fit parameters. The large difference present for $J$ and $l_{3}$ for the 200-mm data is due to the quality of the data present in the primary eclipse. The eclipse is covered by the two cameras on one night, and there is an offset between the data from the two cameras. This maybe be due to the differences in the transmission profile of the filters used in different cameras. In the detrended case, the model shifts enough to favour the data from camera 226, which forms a deeper eclipse. This results in the larger value of $J$ and, because the surface brightness is strongly correlated with $l_{3}$, there is the corresponding decrease in $l_{3}$. Figures 6 and 6 show the best-fit model plotted against the detrended data for the 85-mm and 200-mm data, respectively. The fit residuals are also shown in these figures.

As the sets of parameters from the original and detrended data are generally consistent, but the uncertainties associated with the detrended parameters are smaller, it is the detrended parameters that will be used from now on.

In Table 4, between the 85-mm and 200-mm parameters, there are differences of $2.7\sigma$ and $4\sigma$ for $e$ and $\omega$, respectively. The values also differ from the values quoted by Hełminiak et al. (2009) and Andersen et al. (1988), $e$: 0.187(4) and 0.188(2) respectively; $\omega$: 110.1(9)° and 109.9(6)°, respectively. To investigate further, $e\cos\omega$ and $e\sin\omega$ were fixed at $-0.06424$ and $0.17561$, respectively. These values were calculated from the spectroscopic $e$ and $\omega$ of Hełminiak et al. (2009) and were chosen over Andersen et al. (1988), as the Hełminiak et al. values were obtained more recently. If the orbit has been varying (see Sec. 3.1), the most recent values should be nearest to those applicable to the time span covered by the WASP data.

However, fixing $e\cos\omega$ and $e\sin\omega$ produced models that were poor fits to the data. This was the case for both the 85-mm and the 200-mm data, but was more significant in for the 85-mm. For the 85-mm data, there was a phase offset of 0.001 between the observed data and resulting model of the secondary eclipse.

As an alternative to fixing $e\cos\omega$ and $e\sin\omega$, the values calculated from Hełminiak et al. (2009) with their standard errors were used as Gaussian priors during the model fitting. Table 5 contains the best-fit parameters for the detrended 85-mm and 200-mm data, with inclusion of the priors. Again, the uncertainties have been calculated through MCMC and prayer-bead analysis, and the error contribution from the uncertainties in the limb darkening value has also been included. For the 85-mm data, the inclusion of the priors has altered the best-fit parameters by less than their uncertainties. $e$ and $\omega$ have been altered more significantly for the 200-mm data, bringing their values closer to those of the 85-mm data. However, the values for both are still inconsistent with each other.

The exact cause of the differing $e$ and $\omega$ values remains unclear. Previous observations of AI Phe have yielded a range of values for $e$ and $\omega$. For example, Hrivnak & Milone (1984) found $e=0.1726\pm 0.0006$ and $\omega=(111.8\pm 0.1)\degr$ giving $e\cos\omega=-0.06410\pm 0.00007$ and $e\sin\omega=0.1603\pm 0.0007$, while the same UBVRI lightcurve analysed by Andersen et al. (1988) yielded mean values of $e\cos\omega=-0.064$ and $e\sin\omega=0.183$. No clear trend is present when all available value of $\omega$ were plotted against time, as might be expected if these differences are due to apsidal motion. The parameter $e\sin\omega$ is very sensitive to the shape of the secondary eclipse. Without observing the base of the secondary eclipse in one night, defining the exact shape of the eclipse and therefore determining the value of $e\sin\omega$ can be difficult. Values of $e$ and $\omega$ determined from spectroscopic orbits will not suffer these problems and should therefore be more accurate. It seems there is still work to be done in order to completely understand the behaviour of AI Phe’s orbit.

Despite the inconsistency of $e$ and $\omega$, the addition of the priors has had very little impact on $r_{1}$ and $r_{2}$. They have remained consistent with each other, with a small reduction in the uncertainties of the 200-mm data. The radii are determined by the contact points, which are well defined by the primary eclipse, and the ratio of the eclipses $k$, but only as $k^{0.25}$. Therefore, $r_{1}$ and $r_{2}$ are robustly measured despite problems with the secondary eclipse and small changes in the 200-mm eclipse depth. A number of the other parameters have also shown small reductions in their uncertainties, and therefore the best-fit parameters obtained with the priors have been used in further analysis.

Table 6 summarises the results from the lightcurve analysis, and includes the sets of best-fit parameters from the detrended 85-mm and 200-mm data, obtained with priors. For comparison, the results of Andersen et al. (1988) are also included in the table. As the two WASP data sets and the results of Andersen et al. (1988) are all independent, they have been combined to produce an overall weighted mean. Uncertainties from each of the parameters were used to calculate internal and external standard errors on these mean values (i.e., based on the error bars and scatter, respectively), with the larger of the two being quoted alongside the weighted means. For $r_{1}$ and $r_{2}$ the results were internally consistent. The mean surface brightness ratio and third-light are not included because the different filters were used to obtain the 200-mm and 85-mm data.

For the MCMC analysis, a vector of parameters $\@vec{d}$=($T_{1}$, $\rho_{1}$, $\rho_{2}$, $T_{\rm ratio}$, $M_{\rm sum}$, $q$, ${\rm[Fe/H]_{s}}$) can be used to define the observed quantities for AI Phe. These parameters have been chosen because each quantity is determined by an independent feature of the data used in the analysis, with little or no dependence on the other parameters. $T_{\rm ratio}$ is the ratio of the effective temperatures given by $T_{2}/T_{1}$, $T_{1}$ and $T_{2}$ are the effective temperatures of the two stars, $\rho_{1}$ and $\rho_{2}$ are the average stellar densities of the two stars, $M_{\rm sum}$ is the sum of their masses, $q$ is the mass ratio given by $M_{2}/M_{1}$ and [Fe/H]_s is the observed surface metal abundance.

The mass ratio and sum of the masses were chosen over directly using the individual masses $M_{1}$ and $M_{2}$, as the individual masses have stronger correlations. The densities of the two stars ($\rho_{1}=0.1935\pm 0.0044\,\rho_{\sun}$ and $\rho_{2}=0.0505\pm 0.0007\,\rho_{\sun}$) were calculated using

where $R$ is the radius for star $n=1,2$, $a$ is the semi-major axis of the orbit, $P$ is the orbital period, $G$ is Newton’s gravitational constant (Maxted et al., 2015b). $Q_{\rm n}$ is a function of the mass ratio, where $Q_{1}=q$ and $Q_{2}=1/q$. Equation 8 allows the density to be calculated directly from values of $r_{1}$ and $r_{2}$ derived from the lightcurve analysis using Kepler’s law, and independently from the mass estimates from the spectroscopic orbit. $T_{1}$ was taken to be $6310\pm 150$  K (Vandenberg & Hrivnak, 1985).

The ratio of the stars’ effective temperatures can be determined directly from the surface brightness ratio derived from the lightcurve analysis. The surface brightness ratio is related to the ratio of the eclipse depths in a totally-eclipsing binary, with very little dependence on the other parameters of the lightcurve. An approximate value of $T_{\rm eff}$ is needed for one of the stars in the binary, but this value has only a small effect on the derived $T_{\rm ratio}$. Our method uses Kurucz (1993) model atmospheres and the profiles for numerous passbands, (Bessell, 1990; Crawford & Barnes, 1970; Doi et al., 2010) and is similar to the method described by Maxted et al. (2015a). For each of the bands, Johnson $BVRI$, Strömgen $y$ and SDSS $r^{\prime}$, a relationship was established between effective temperature and surface brightness for $\log g=4.0$ and $\log g=3.6$. Interpolation between the values from the models for $\log g=4.0$ and $\log g=3.5$ was used in the case of $\log g=3.6$ as no model was available. Taking $T_{1}=6310\pm 150$  K we used these relationships to find a value for $T_{2}$ that gave the measured average surface brightness ratio for each band. The average surface brightness ratios were calculated from the central surface brightness ratios for the BVRI and y passbands from Andersen et al. (1988) and the 85-mm central surface brightness ratio was used for SDSS $r^{\prime}$ passband. Andersen et al. (1988) did not define their surface brightness ratios as average or central values. We have assumed the values are central surface brightness ratios because their lightcurve model uses the same methods as ebop and should therefore produce the same ratios as ebop. Following the example of Andersen et al. (1988), limb darkening coefficients for the BVRI were taken from Hrivnak & Milone (1984) and limb darkening coefficients for y from Wade & Rucinski (1985). The values of $T_{\rm ratio}$ derived using this method for the different passbands were found to be consistent with each other. The weighted mean and standard error are $T_{\rm ratio}=0.83\pm 0.01$, where the standard error estimate accounts for the uncertainties on both $J$ and $T_{1}$. There are surface brightness ratios available for the Johnson $U$ and Strömgen $uvb$ passbands in the Andersen et al. (1988) paper, however this method is less reliable in the bluer passbands as line-blanketing is more prevalent. These bands were therefore excluded from the determination of $T_{\rm ratio}$. Another consideration is the metallicity used by the models. The Kurucz (1993) model atmospheres consider solar metallicity, instead of the measured value of AI Phe ($-0.14\pm 0.1$, Andersen et al., 1988). However, Maxted et al. (2015a) found that changing the metallicity by $\pm 0.1$ dex changes the resulting $T_{2}$ by less than 10 K, or $\approx 0.002$ in $T_{\rm ratio}$. As such, the resulting error is within the uncertainties of $T_{\rm ratio}$.

The model parameters used to predict the observed data can be represented as $\@vec{m}=\left(\mbox{$\tau_{\rm sys}$},\mbox{$M_{1}$},\mbox{$M_{2}$},\mathrm{[Fe/H]}_{\mathrm{i}}\right)$, where $\tau_{\rm sys}$, $M_{1}$, $M_{2}$ and $\mathrm{[Fe/H]}_{\mathrm{i}}$ are the age, mass of star 1, mass of star 2, and initial metal abundance of the system, respectively. Due to diffusion and mixing processes occurring in the star during its evolution, the initial metal abundance, $\mathrm{[Fe/H]}_{\mathrm{i}}$ differs from the observed surface metal abundance, $\mathrm{[Fe/H]}_{\mathrm{s}}$.

The probability distribution function $p(\@vec{m}|\@vec{d})\propto{\cal L}(\@vec{d}|\@vec{m})p(\@vec{m})$ was determined using a MCMC method. ${\cal L}(\@vec{d}|\@vec{m})=\exp(-\chi^{2}/2)$ is used to estimate the likelihood of observing the data $\@vec{d}$ for a given model $\@vec{m}$, where

The Markov chain setup uses the same approach described by Maxted et al. (2015b), in that a Markov chain of points $\@vec{m}_{i}$ is created with the probability distribution $p(\@vec{m}|\@vec{d})$ using a jump probability distribution $f(\Delta\@vec{m})$. The generation of each trial point $\@vec{m}^{\prime}=\@vec{m}_{i}+\Delta\@vec{m}$ is dictated by probability distribution, with the point being rejected if any of the model parameters are outside of the ranges set by the priors. If ${\cal L}(\@vec{m}^{\prime}|\@vec{d})>{\cal L}(\@vec{m}_{i}|\@vec{d})$, the point is always included, and when ${\cal L}(\@vec{m}^{\prime}|\@vec{d})<{\cal L}(\@vec{m}_{i}|\@vec{d})$ the point is included with probability ${\cal L}(\@vec{m}^{\prime}|\@vec{d})/{\cal L}(\@vec{m}_{i}|\@vec{d})$ (Metropolis-Hastings algorithm). $\@vec{m}_{i+1}=\@vec{m}^{\prime}$ if the trial point is accepted, or $\@vec{m}_{i+1}=\@vec{m}_{i}$ otherwise (Metropolis et al., 1953; Hastings, 1970).

In comparison to Maxted et al. (2015b), a different approach is used to find a starting point. Their approach takes the point with the lowest $\chi^{2}$ when a sample of points are randomly generated across the model grid. In our case, because the parameter-space is so well constrained, the starting point needed to be guided towards the parameter-space. The measured masses and metallicity are fixed and used to generate an evolutionary track for each star in the system. From there, each point along the tracks is searched to find the point with the lowest $\chi^{2}$, using 2 000 steps in age. This is done separately for the main-sequence and post main-sequence stages, with the best-fit point from this process is used as an initial starting point.

For each parameter a step size is found such that $|\ln({\cal L}(\@vec{m}_{0}|\@vec{d})-\ln({\cal L}(\@vec{m}_{0}^{(j)}|\@vec{d})|\approx 0.5$, where $\@vec{m}_{0}^{(j)}$ denotes a set of model parameters that differs from $\@vec{m}_{0}$ only in the value of one parameter, $j$. From there, a burn-in stage of 10 000 steps is used to improve the initial set of parameters and to determine correlations between parameters by calculating a covariance matrix. The eigenvectors and eigenvalues of this matrix are used to determine a set of uncorrelated, transformed parameters, $\@vec{q}=(q_{1},q_{2},q_{3},q_{4})$, where each of the transformed parameters has unit variance (Tegmark et al., 2004). To estimate the probability distribution $p(\@vec{m}|\@vec{d})$, a second Markov chain of 1 000 000 steps is produced, using a Gaussian distribution for $f(\Delta\@vec{q})$ with unit standard deviation for each of the transformed parameters and the most likely model parameters as a starting point. With the tight constraints given by the priors, some of the model grids have a severely restricted parameter-space and so required a large number of steps to ensure the space was explored. All chains were visually inspected and checked via a running mean to ensure suitable mixing.

Using the Bayesian method described in Sects. 5.1 and 5.2 an age distribution has been produced for a number of model grids, with different mixing lengths $\alpha_{\rm ml}$ and helium abundance $Y$. In total, distributions were produced for six mixing lengths (1.22, 1.36, 1.50, 1.78, 2.04 and 2.32) with $\Delta Y$ held fixed at zero (see Fig. 7) and an additional ten distributions were produced for helium abundances, where $\Delta Y$ ranged from $-0.05$ to $0.05$ in increments of 0.01, with the mixing length held fixed at 1.78 (Fig. 8). Table 8 contains the parameters from the model with the lowest $\chi^{2}$ for each model grid, alongside the mean $\tau_{\rm mean}$ and standard deviation $\sigma_{\tau_{\rm mean}}$ of the age distributions produced by the 1 000 000-step chain. In some cases the age distributions produced bimodal distributions which do not fit a Gaussian profile, meaning there is a noticeable difference between the best-fit age $\tau_{\rm best}$ and $\tau_{\rm mean}$. The first section of the table contains the model grids where ${\alpha_{\rm ml}}$ is varied, and the second section contains model grids where the helium abundance, $Y_{\rm AI}$, is varied by $\Delta Y$ from the initial value $Y_{0}$.

The $\chi^{2}$ values in Table 8 show that some of the model grids fit in the parameter-space much better than others. In terms of the mixing lengths, values of 2.04 and 1.78 are favoured, with 1.78 producing the lowest $\chi^{2}$. Using the 16th, 50th and 84th percentiles, a mixing length of 1.78 gives an age of $4.35^{+0.23}_{-0.19}$ Gyr. As the mixing length in the model grid is increased, the best-fit model tends to increase the mass of the two stars. Comparing evolutionary tracks for star 2 (the sub-giant) at a fixed age, mass and metallicity but with increasing mixing length, results in larger densities, as the stars are more compact (Lebreton et al., 2014). However, with the tight observational constraints on $\rho_{2}$, the models attempt to improve the overall $\chi^{2}$ by increasing $M_{2}$ and by decreasing the metallicity. Increasing $M_{2}$ results in a star with a larger radius with the same density, but also reduces the age of the sub-giant and therefore the system. The effect on the age can be seen in the mean and best-fit ages as $\alpha_{\rm ml}$ is changed from 1.50 to 2.32. Decreasing the metallicity allows the system to evolve to the same evolutionary stage faster.

The $\alpha_{\rm ml}=1.22$ model grid does not follow the mass trend mentioned above, because the mass distribution is bi-modal, corresponding to two different evolutionary stages. The most dominant part of the distribution places the secondary very early in the contraction phase before the ascent to the red-giant branch, with a larger effective temperature. The less prominent part of the distribution, places the secondary further into the contraction phase. A similar explanation can be used to explain the sudden change in best-fit age between $\alpha_{\rm ml}=1.36$ and $\alpha_{\rm ml}=1.50$. For the $\alpha_{\rm ml}=1.36$ model, the secondary component sits firmly in the contraction phase, whereas for $\alpha_{\rm ml}\geqslant 1.50$, the secondary sits at the base of the red-giant branch.

Increasing the helium abundance used in the model grid has decreased the best-fit age of the system, with the models preferring to use smaller masses and cooler effective temperatures for both stars. The best-fit values for $\rho_{2}$ show very little variation with helium abundance, changing by less than 1%. Meanwhile, $\rho_{1}$ shows a much larger variation of 14%. The small uncertainty in $M_{2}$ tightly constraints the age of the sub-giant, allowing little variation in $\rho_{2}$. This is not the case with the main-sequence star. In order to find a model for the main-sequence star that fits the age determined by the sub-giant, parameters such as its effective temperature and density need to be varied significantly from their observed value for some values of the helium abundance. This is reflected in their much larger $\chi^{2}$ values.

Looking at the $\chi^{2}$ for the model grid with varying $\Delta Y$ in Table 8, the preferred models have a helium abundance that is closer to the initial value, effectively excluding values where $\Delta Y\leqslant-0.03$ and $\Delta Y\geqslant 0.04$. The preferred model is where $\Delta Y=0.0$, meaning $Y_{\rm AI}=0.261$. Using $\Delta\chi^{2}=1$ to define a 68.3% confidence interval on $\Delta Y$, based on the projection into the $\alpha_{\rm ml}$-$\Delta Y$ parameter space (Press et al., 1992), the models for $\Delta Y=-0.01$ and $\Delta Y=0.02$ fall just outside this interval. A $\Delta\chi^{2}=2.71$ defines a 90% confidence interval and covers the models $-0.02<\Delta Y<0.03$. We find $Y_{\rm AI}=0.26^{+0.02}_{-0.01}$ for a fixed mixing length of $1.78$, giving an age of $4.39\pm 0.32$ Gyr using $\tau_{\rm best}$ from the best-fit model. The uncertainty is estimated by directly adding $\sigma_{\tau_{\rm mean}}$ from $\Delta Y=0$ ($0.20\,$Gyr) and a systematic uncertainty of 0.12 Gyr (from not knowing $Y_{\rm AI}$). For comparison, Torres et al. (2010) found an age of 4.1 Gyr using experimental Victoria models (VandenBerg et al., 2006) and 5.0 Gyr from Yonsei-Yale models (Demarque et al., 2004), although no uncertainties are given. Spada et al. (2013) found the age of two components of AI Phe separately using an updated version of the Yale Rotational stellar Evolution Code (YREC). They find an age of $4.44\pm 0.08\,$Gyr for the hotter component, and an age of $4.54\pm 0.02\,$Gyr for the cooler component meaning our value is consistent with their ages.

Comparisons between mixing lengths are quite difficult because of the different approaches to used to calibrate the parameter, and how the mixing length is included in the models. Andersen et al. (1988) used a fixed mixing length of 1.50, which was calibrated by the producing a $1.0\,M_{\sun}$ star with $1.0\,R_{\sun}$ and the solar age of $4.7\,{\rm Gyr}$. They used interpolation between a solar model and a model with a $Z$ of 0.01 to find the age and helium abundance of AI Phe. Spada et al. (2013) used a similar solar calibration method but found the mixing-length was $\alpha=1.875$, and Torres et al. (2010) noted that at the time the Victoria models (VandenBerg et al., 2006) that had not been fully calibrated. Therefore, we note that our best-fit mixing length is the same as the solar value for our set of models. This is consistent with what was found by Andersen et al. (1988).

Figures 9 and 10 show the theoretical tracks for the two components of AI Phe which result in the best-fit age for the model grid with $\alpha_{\rm ml}=1.78$ and $\Delta Y=0.0$. In Fig. 10, we have plotted the observed densities calculated in this work and the effective temperatures from Milone et al. (1992).