THE CLUSTER AGES EXPERIMENT (CASE). VII. ANALYSIS OF TWO ECLIPSING BINARIES IN THE GLOBULAR CLUSTER NGC 6362^*

The Cluster AgeS Experiment is a long term project conducted mainly at Las Campanas Observatory with a principal aim of the detection and detailed study of detached eclipsing binaries (DEBs) in nearby globular clusters (GCs). The main goal is to measure masses, luminosities and radii of DEB components with a precision of better than 1% in order to determine GC ages and distances, and to test stellar evolution models (Kaluzny et al. 2005). The methods and assumptions we employ follow the ideas of Paczyński (1997) and Thompson et al. (2001). Thus far, we have presented results for seven DEBs from four clusters: ω Cen, 47 Tuc, M4 and M55. These were the first and, to the best of our knowledge, the only direct measurements of the fundamental parameters of main-sequence and subgiant stars in GCs. More details on the project and references to earlier papers can be found in Kaluzny et al. (2014a).

The present paper is devoted to an analysis of two DEBs, V40 and V41, belonging to the GC NGC 6362. Radial velocity and photometric observations are used to derive masses, luminosities, and radii of the component stars. In Section 2 we describe the photometric observations and determine orbital ephemerides. Section 3 presents spectroscopic observations and radial-velocity measurements. The combined photometric and spectroscopic solutions for the orbital elements and the resulting component parameters are obtained in Section 4. In Section 5, we compare the derived parameters to Dartmouth stellar evolution models, with an emphasis on estimating the age of the binaries. Finally, in Section 6 we summarize our findings.

Photometric observations analyzed in this paper were obtained at Las Campanas Observatory between 1995 and 2011 using the 2.5 m du Pont telescope equipped with the TEK5 and SITE2K CCD cameras, and the 1.0 m Swope telescope equipped with the SITE3 CCD camera. The scale was 0.26 arcsec pixel⁻¹ and 0.44 arcsec pixel⁻¹ for the du Pont and Swope telescopes, respectively. For each telescope the same set of BV filters was used for the entire duration of the project. The first likely eclipse events of V40 and V41 were observed in 1991 June. For the next eight years observations were conducted in a survey mode with the aim of detecting more eclipses and establishing ephemerides for both variables. The periodic variability of V40 and V41 was firmly established by Mazur et al. (1999). Starting with the 2001 season, we attempted to cover predicted eclipses of V40 or V41. Several eclipses were observed with both telescopes. The light curve analysis in Section 4 is based exclusively on the du Pont data. The less accurate Swope photometry was used solely to improve the ephemerides (the difference in accuracy between the two telescopes resulted from the difference in the light-gathering power and the poorer sampling of the point-spread function (PSF) of the Swope telescope).

The core and half-light radii of NGC 6362 are r_c = 70'' and r_h = 123'', respectively (Harris 1996, 2010 edition). V40 is located in the core region at a distance of d = 59'' from the cluster center. V41 resides at the half-light radius with a distance of d = 121''. The du Pont images show no evidence of unresolved visual companions to either binary. In the case of V40 this is confirmed by the examination of Hubble Space Telescope (HST)/Advanced Camera for Surveys (ACS) images obtained within the program GO-10775 (PI Ata Sarajedini). Both systems are proper-motion members of NGC 6362 (Zloczewski et al. 2012; Kaluzny et al. 2014b). Coordinates and finding charts for V40 and V41 can be found in Mazur et al. (1999).

The light curves presented here were measured with profile fitting photometry extracted with the Daophot/Allstar codes (Stetson 1987, 1990), and transformed to the standard Johnson system based on observations of Landolt (1992) standards. A more detailed description of the data and the reduction procedures can be found in a paper presenting the photometry of 69 variables from the field of NGC 6362 (Kaluzny et al. 2014b). We note that, in general, more accurate photometry of crowded fields can be obtained with an image subtraction technique. However, profile fitting photometry is less vulnerable to systematic errors, and as such it is a better option whenever the star is free from blending.

2.1. Ephemerides

Several eclipse events of V40 were covered densely enough to determine individual times of minima. They were found using the KWALEX code which is based on an improved version of the algorithm presented by Kwee & van Woerden (1956) (see appendix). The linear ephemerides:

$$\begin{eqnarray}{\mathrm{HJD}}^{\mathrm{MinI}} & = & 2453915.50966(25)+5.2961725(14)\times E\\ {\mathrm{HJD}}^{\mathrm{MinII}} & = & 2453891.82853(14)+5.29617486(99)\times E\end{eqnarray} \tag{ 1 }$$

proved to be an adequate fit to the moments of minima listed in Table 1. The orbit of V40 is eccentric. Since the secondary (shallower) minimum was observed more times than the primary one, we decided to choose it as a base of the phase count. Thus, we set ϕ^MinII = 0.5, which necessarily results in ${\phi }^{\mathrm{MinI}}\ne 0.$

Table 1.  Times of Minima for V40

E	HJD-2450000	σ	O−C
−480.5	1349.66400	0.00129	0.00060
−274.5	2440.67681	0.00083	−0.00019
0.5	3891.82809	0.00027	0.00045
3.5	3907.71733	0.00024	−0.00027
16.5	3976.56802	0.00139	−0.00069
77.5	4299.63424	0.00026	−0.00024
138.5	4622.70052	0.00025	0.00015
276.5	5353.57276	0.00031	0.00003
0.	3915.50987	0.00028	−0.00022
1.	3920.80387	0.00096	0.00196
139.0	4651.67750	0.00035	0.00014
335.0	5689.72753	0.00039	−0.00007

Download table as:  ASCIITypeset image

Altogether, 13 eclipses of V41 were observed between 2007 August and 2011 September. However, ingress and egress were only covered for a few of them. This prevented us from a classical period study based on the determination of individual moments of minima and the O−C technique. To find the ephemeris, a combined V light curve including data from all seasons and both telescopes was used. This was fitted with the help of the JKTEBOP v34 code⁷ (Southworth 2013, and references therein), allowing only for variations of P and moment of the primary eclipse T₀. Mode 8 of the code was used, which enabled a robust determination of the errors of the two parameters with the help of Monte Carlo simulations. The following linear ephemeris

$$\begin{eqnarray}&&{\mathrm{HJD}}^{\mathrm{MinI}}=2454263.65242(25)+17.8888441(40)\times E\end{eqnarray} \tag{ 2 }$$

was obtained. The light curves phased with ephemerides given by Equations (1) and (2), and used for the subsequent analysis, are shown in Figure 1.

Zoom In Zoom Out Reset image size

The spectra were taken with the MIKE echelle spectrograph (Bernstein et al. 2003) on the Magellan Baade and Clay 6.5 m telescopes, using a 07 slit which provided a resolution R ≈ 40,000. A typical observation consisted of two 1800 s exposures of the target, flanking an exposure of a thorium–argon hollow-cathode lamp. The raw spectra were reduced with the pipeline software written by Dan Kelson, following the approach outlined in Kelson (2003). For further processing, the IRAF package ECHELLE was used. The velocities were measured using software based on the TODCOR algorithm (Zucker & Mazeh 1994), kindly made available by Guillermo Torres. Synthetic echelle-resolution spectra from the library of Coelho et al. (2005) were used as velocity templates. These templates were interpolated to the values of $\mathrm{log}\;g$ and T_eff derived from the photometric solution (see Section 4) and Gaussian-smoothed to match the resolution of the observed spectra. For the interpolation in chemical composition we assumed [Fe/H] = −1.07 and [α/Fe] = 0.4, i.e., values adopted for the cluster as a whole (Carretta et al. 2009; Dotter et al. 2010). These values are not well determined (see the discussion in Section 4), however the results of the velocity measurements are insensitive to minor changes in [Fe/H] or [α/Fe].

Radial velocities were measured on the wavelength intervals 4120–4320 and 4350–4600 Å, covering the region of the the best signal-to-noise ratio (S/N) while avoiding the Hγ line. The final velocities, obtained by taking the mean of these two measurements, are given in Tables 2 and 3, which list the heliocentric Julian dates at mid-exposure, velocities of the primary and secondary components, and the orbital phases of the observations calculated according to the ephemerides given by Equations (1) and (2). An estimate of the error in the velocity measurements is the rms of half of the differences between velocities obtained from the two wavelength intervals. The corresponding values are 0.63, 0.59, 0.23 and 0.52 km s⁻¹ for the primary and secondary of V40 and the primary and secondary of V41, respectively. They are consistent with the rms values of the residuals from orbital fits given in Table 4.

Table 4.  Orbital and Photometric Parameters^a of the Systems V40 and V41

Parameter	V40	V41
γ (km s−1)	−12.34(11)	−12.55(8)
Kp (km s−1)	70.18(17)	46.59(14)
Ks (km s−1)	73.63(25)	52.57(17)
σp (km s−1)	0.72	0.53
σs (km s−1)	1.07	0.65
A (R⊙)	15.037(32)	33.294(71)
e	0.05054(73)	0.3125(13)
ω (deg)	206.6(17)	321.23(33)
i (deg)	88.22(7)	89.547(20)
rp	0.08817(45)	0.03225(13)
rs	0.06627(87)	0.02195(12)
SV	0.9815(35)	0.7333(63)
SB	0.9686(18)	0.669(11)
(Lp/Ls)V	1.803(54)	2.964(27)
${({L}_{{\rm{p}}}/{L}_{{\rm{s}}})}_{B}$	1.8272(57 )	3.264(44)
σrms(V) (mmag)	11	12
σrms(B) (mmag)	13	16
Vp (mag)b	18.698(13)(20)	19.089(7)(17)
Vs (mag)b	19.338(22)(27)	20.274(10)(18)
${(B-V)}_{{\rm{p}}}$ (mag)b	0.542(18)(23)	0.550(10)(18)
${(B-V)}_{{\rm{s}}}$ (mag)b	0.556(31)(34)	0.650(16)(22)

Notes.

^aNumbers in parentheses are the errors of the last significant digit(s). ^bFor V and (B − V) both the internal error (from the photometric solution and profile photometry) and the total error is given, the latter including 0.015 mag uncertainty of the zero point of the magnitude scale.

Download table as:  ASCIITypeset image

We assume that reddening is uniform for NGC 6362. Harris (1996, 2010 edition) lists $E(B-V)=0.09$ as measured by Zinn & West (1984) based on the Q39 index. The map of total galactic extinction by Schlegel et al. (1998) predicts $E(B-V)=0.075$ mag for the position of NGC 6362. This value agrees with an estimate resulting from isochrone fitting applied to HST/ACS photometry of the cluster by Dotter et al. (2010). They obtained E(6–8) = 0.070, equivalent to $E(B-V)=0.071.$ This value is also compatible with values obtained by Kovács & Walker (2001) ($E(B-V)=0.073\pm 0.007$) and Olech et al. (2001) ($E(B-V)=0.08\pm 0.01$) from an analysis of RR Lyr stars in NGC 6362. We adopt $E(B-V)=0.075$ and neglect the differential reddening, which according to (Bonatto et al. 2013) is 0.025 mag on the average. (Bonatto et al. 2013) comment that this differential reddening may result from uncorrected zero-point variations in the original photometry.

We employed a quadratic formula for limb darkening, and calculated coefficients iteratively together with effective temperatures of the components of both systems. The temperatures were derived from dereddened B − V colors using the empirical calibration of Casagrande et al. (2010). The coefficients were then interpolated from the tables of Claret (2000) with the help of the JKTLD code.⁸ The updated colors in turn were computed from total V and B magnitudes at maximum light and light ratios derived from the light curve solution. We note that a slightly more recent calibration by Sousa et al. (2011) results in temperatures higher by 8 K on average, however the agreement between calibrations is excellent over the relevant range of colors.

We checked that using a square root or logarithmic formula for the limb darkening instead of a quadratic formula had a negligible impact on the photometric solution. In particular, relative radii derived from the three two-parameter formulas differed by less than ∼0.2σ), whereas using a linear limb darkening formula resulted in a slightly poorer fit, and caused the relative radii to differ by ∼1σ from those obtained with the quadratic formula. We conclude that, for the stars analyzed here, two-parameter formulate provide a better description of the limb darkening than the linear law.

Light and radial velocity curves of V40 and V41 were analyzed with the JKTEBOP v34 code, which allows for a simultaneous solution of one light curve and two radial velocity curves (one for the primary and one for the secondary component). The V-curve was solved for the following parameters: inclination i, sum of fractional radii r_p and r_s, ratio of the radii r_s/r_p; surface brightness ratio S (secondary to primary), luminosity ratio L_s/L_p, $e\mathrm{cos}(\omega )$ and $e\mathrm{sin}(\omega )$ where e is the orbital eccentricity and ω the longitude of the periastron, amplitudes of radial velocity curves K_p and K_s, and systemic velocity γ. Errors of stellar parameters were derived using a Monte Carlo approach implemented in mode 8 of the JKTLD code. We performed 10,000 Monte Carlo simulations for each of V40 and V41.

Since the B-data were of poorer quality than those collected in the V-band, we decided to use them solely to a measurement of the luminosity ratio (L_s/L_p)_B. This was accomplished by fixing all parameters but the surface brightness ratio S_B. However, the formal error of S_B calculated by JKTEBOP would have been seriously underestimated if uncertainties of other parameters had not been taken into account. To avoid this problem we used the PHOEBE package v0.31a (Prša & Zwitter 2005), which allows for a simultaneous solution of several light curves but can be unstable when attempting to solve light and velocity curves simultaneously. We initialized the package with the photometric solution found by JKTEBOP, loaded BV curves, and, allowing all geometrical parameters to vary, found the standard error of S_B using another Monte Carlo procedure written in PHOEBE-scripter. The procedure replaces the observed light curves B₀ and V₀ with the fitted values B_f and V_f, generates Gaussian perturbations δB_f and δV_f such that the standard deviation of the perturbation is equal to the standard deviation of the corresponding residuals shown in Figure 2, and performs PHOEBE iterations on ${B}_{{\rm{f}}}+\delta {B}_{{\rm{f}}}$ and ${V}_{{\rm{f}}}+\delta {V}_{{\rm{f}}}.$ As before, 10,000 Monte Carlo simulations were performed for each binary. In our opinion such a hybrid approach based on both the JKTEBOP and PHOEBE codes is optimal when analyzing data sets similar to those presented in this paper.

Zoom In Zoom Out Reset image size

The quality of the final fits is illustrated in Figures 2 and 3. Orbital and photometric parameters derived for both systems are presented in Table 4, and absolute parameters of the components are listed in Table 5. The errors of the effective temperatures include uncertainties of observed colors, reddening, and photometric zero points (0.015 mag for V and B − V). The bolometric luminosities are calculated from radii and temperatures adopting ${L}_{\odot }^{\mathrm{bol}}=3.827\times {10}^{33}$ erg.⁹ To derive absolute visual magnitudes, we used bolometric corrections based on atmosphere models of Castelli & Kurucz (2004) which ranged from −0.07 mag for the primary of V40 to −0.10 mag for the secondary of V41. The absolute visual magnitudes along with the V-magnitudes derived from light curve solutions allowed for a determination of the apparent distance moduli listed in Table 5. The location of the components of V40 and V41 on the color–magnitude diagram (CMD) of NGC 6362 is shown in Figure 4. All four stars are placed on the main sequence: three of them at the turnoff region, and the fourth one (the secondary of V41) well below the turnoff.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Absolute Parameters^a of the Components of V40 and V41

Parameter	V40	V41
Mp (M⊙)	0.8337(63)	0.8215(58)
Ms (M⊙)	0.7947(48)	0.7280(47)
Rp (R⊙)	1.3253(77)	1.0739(48)
Rs (R⊙)	0.997(13)	0.7307(46)
Tp (K)	6156(123)	6124(121)
Ts (K)	6100(158)	5747(122)
Lbolp (L⊙)	2.27(18)	1.46(11)
Lbols (L⊙)	1.24(13)	0.524(45)
${M}_{{V}_{{\rm{p}}}}$ (mag)	3.92(8)	4.40(8)
${M}_{{V}_{{\rm{s}}}}$ (mag)	4.59(11)	5.54(9)
${(m-M)}_{{V}_{{\rm{p}}}}$ (mag)	14.78(8)	14.69(8)
${(m-M)}_{{V}_{{\rm{s}}}}$ (mag)	14.75(9)	14.73(8)
Agep (Gyr)	11.7(3)	10.7(3)
Ages (Gyr)	11.5(4)	9.2(6)
Ageb (Gyr)	11.7(2)	10.4(3)

Note.

^aNumbers in parentheses are the errors of the last significant digit(s). The age of the binary Age_b is a weighted mean of component ages.

Download table as:  ASCIITypeset image

The fundamental parameters obtained from the analysis of V40 and V41 allow us to derive the ages of the components of these systems, using theoretical isochrones. However, before a model comparison is made it is necessary to review the information available on the chemical composition of NGC 6362 since the model-based ages are sensitive to the adopted values of helium abundance, [Fe/H] and [α/Fe].

Unfortunately, the metallicity of NGC 6362 is not well known. A literature search reveals that determinations of [Fe/H] from high resolution spectra are limited to the study by Carretta & Gratton (1997), who found [Fe/H] = −0.96 and [Fe/H] = −0.97 for two red giants. Carretta et al. (2009) used this determination along with some older estimates based on the integrated photometry of NGC 6362 to derive [Fe/H] = −1.07 ± 0.05 on their newly proposed metallicity scale. For the sake of completeness, we note that Harris (1996, 2010 edition) listed [Fe/H] = −0.99—a value obtained by combining that of Carretta et al. (2009) with [Fe/H] = −0.74 ± 0.05 derived by Geisler et al. (1997) based on Washington photometry of 10 cluster giants. Olech et al. (2001) obtained a helium content Y = 0.292 from the analysis of RRc stars. Finally, Dalessandro et al. (2014) presented convincing evidence of multiple populations in NGC 6362. We adopt a fiducial composition of [Fe/H] = −1.07, [α/Fe] = +0.4 and Y = 0.25, but will investigate the dependence on these parameters when deriving the ages of the component stars.

In a study of three binary systems in M4 (Kaluzny et al. 2013) we used two sets of theoretical isochrones: Dartmouth (Dotter et al. 2008, henceforth DSED) and Victoria-Regina (VandenBerg et al. 2012). Since the differences in the derived ages turned out to be smaller than the uncertainties imposed by the observational errors, the present analysis is based on DSED isochrones only. Figure 5 shows DSED isochrones for ages ranging from 9 Gyr to 13 Gyr plotted on mass–radius and mass–luminosity diagrams, together with values for V40 and V41. Both diagrams indicate that V41 may be younger than V40, and the mass–radius diagram suggests that the secondary of V41 is younger than the primary (or, if both components are of the same age, the secondary is too small compared to the primary).

Zoom In Zoom Out Reset image size

A detailed comparison of the ages of the individual components, which requires a quantitative age estimate from the data and the isochrones, was undertaken following the approach developed by Kaluzny et al. (2013) to study the binary systems in M4. Briefly, the (M−R) or (M−L) plane is densely sampled in both dimensions for each component separately, and the age of an isochrone matching the properties of the sampled point is found. The age ensembles generated in this way are visualized in the form of mass–age diagrams (MADs). The MAD of a given component is composed of two ellipses derived from (M−R) and (M−L) relations. Each ellipse defines a range of ages achieved by stars whose parameters differ from the corresponding entries (M, R) or (M, L) in Table 5 by at most _R or _L, where ${\varepsilon }_{R}^{2}=\sigma {(M/{M}_{\odot })}^{2}+\sigma {(R/{R}_{\odot })}^{2}$ for the (M−R) plane and likewise for the (M−L) plane; σ(M/M_⊙), σ(R/R_⊙) and σ(L/L_⊙) being the standard errors of M, R and L from Table 5. The MADs of the components of V40 and V41 are shown in Figure 6.

Zoom In Zoom Out Reset image size

The age of each component was calculated as the weighted mean of the ages derived from the two corresponding MADs. For V40 we obtained 11.7 ± 0.3 Gyr for the primary and 11.5 ± 0.4 Gyr for the secondary. The difference is 0.2 ± 0.5 Gyr. For V41 the ages are 10.7 ± 0.3 Gyr and 9.2 ± 0.6 Gyr, respectively, for the primary and the secondary. The difference is 1.5 ± 0.7 Gyr. Possible origins of this two-sigma effect are discussed in Section 6. We estimated the age of each system from the weighted mean of the ages derived from each pair of corresponding MADs. For V40 we obtained 11.7 ± 0.2 Gyr, in agreement with 12.5 ± 0.5 Gyr derived by Dotter et al. (2010) from CMD fitting. For V41 we obtained 10.4 ± 0.3 Gyr, which suggests that this system is significantly younger than V40, with a formal difference of 1.3 ± 0.4 Gyr.

To independently estimate the age of the cluster as a whole, we fitted the CMD of NGC 6362 with DSED isochrones (see Figure 7). Since uncertainties of stellar models grow larger with evolutionary time, when fitting we assigned the largest weight to main sequence and turnoff stars. This results in an age estimate of 12.5 ± 0.75 Gyr—again in agreement with the result of Dotter et al. (2010), which is based on an independent set of photometric data. The CMD fit yielded a distance modulus of 14.72 mag—larger than 14.55 mag derived by Dotter et al. (2010), but in good agreement with 14.68 mag listed by Harris (1996, 2010 edition). The reddening derived from the CMD-fit ($E(B-V)=0.061$ mag) is slightly lower than the value of 0.075 mag adopted in Section 4, but decreasing [Fe/H] from the fiducial −1.07 to −1.15 brings it back to the adopted value (we checked that such a small change in metallicity does not influence the parameters listed in Tables 4 and 5 in any significant way). The distance moduli of all four stars forming V40 and V41 agree with each other (see Table 5), and the mean weighted modulus of 14.74 ± 0.04 mag is consistent with that obtained from CMD fitting.

Zoom In Zoom Out Reset image size

To sum up, the age derived from the MADs of the V40 components is compatible with that derived from CMD fitting, however the age of the primary of V41 is 1.0 ± 0.4 Gyr younger than the age of V40, with the secondary of V41 being another 1.5 ± 0.7 Gyr younger. The problem looks much milder in the (M−L) diagram, in which all the four components are marginally coeval within 1σ limits. We stress, however, that the (M−L) diagram and associated MADs must be given a smaller weight because the luminosities in Table 5 are sensitive to adopted chemical composition (via T_eff, see Dotter et al. 2010) and possible differential reddening, whereas the radii are entirely reddening and composition-independent.

We derived absolute parameters of the components of V40 and V41, two DEBs located near the main sequence turnoff of NGC 6362. The accuracy of mass and radii measurements is generally better than 1%, with the only exception being the radius of the secondary in V40 which is determined with an accuracy of 1.3%. The luminosities of the components were found using effective temperatures estimated from (B − V)–T_eff calibrations compiled by Casagrande et al. (2010). While their errors are rather large (between 7% and 10%), all the four distance moduli have nearly the same value, and the formally calculated mean modulus of 14.74 ± 0.04 mag agrees very well with 14.72 mag derived from the comparison of DSED isochrones with the CMD of the cluster.

Based on the parameters of V40 and DSED isochrones for Y = 0.25, [Fe/H] = −1.07 and [α/Fe] = +0.4, we derived an age of 11.7 ± 0.2 Gyr. The good agreement of this result with the age of NGC 6362 derived from isochrone comparison with the CMD of the cluster is encouraging. This is particularly the case given the age of the cluster as obtained by Dotter et al. (2010) based on an independent photometry (12.5 ± 0.5 Gyr). However, one should keep in mind that the uncertainty of 0.2 Gyr is a purely statistical one, and, as pointed out by Kaluzny et al. (2013), it should be increased by ∼0.85 Gyr to account for the sensitivity of stellar age to uncertainties in the chemical composition.

The situation with V41 is less clear. Isochrone fits to (M−R) and (M−L) diagrams of this system suggest that it is significantly (by 1.3 ± 0.4 Gyr) younger than V40. However, given the recent discovery of two stellar populations in NGC 6362, such a possibility does not seem to be entirely implausible. Upon examining CMDs obtained from HST photometry, Dalessandro et al. (2014) found that the red giant branch (RGB) of the cluster splits into two sub-branches separated by about 0.1 mag in (U − V), while such effect is not observed in (V − I). Their Figure 1 explains the split by a difference in light-element content. While they use isochrones of the same age, which is reasonable given the relatively weak dependence of RGB location on time, they clearly identify the Na-poor populations with the first generation of stars, and the Na-rich population with the second one "formed during the first few ∼100 Myr of the cluster life from an intra-cluster medium polluted by the first generation." We observed similar effects (i.e., a large split in (U − V) and a much smaller one in (V − I)) using same-age DSED isochrones differing by ∼0.1 dex in [Fe/H], or ∼0.3 dex in [α/Fe], or ∼0.1 in Y. The RGB also splits when the populations have the same composition but differ in age by 4–5 Gyr, however an age difference this large causes the subgiant branch to split even more widely than the RGB, while the opposite is observed in NGC 6362. Thus, we corroborate the results of Dalessandro et al. (2014) and their scenario involving two star formation episodes, with the second-generation stars slightly enriched in metals (and possibly in helium).

The discrepancy of the ages of V41 components (1.5 ± 0.7 Gyr) may be regarded as statistically insignificant. However we decided to study the possible origin of the difference by varying chemical composition parameters Y, [Fe/H] and [α/Fe] (only one of them was varied at a time), and looking for isochrones joining locations of the primary and the secondary in the (M−R) and (M−L) diagrams. We found acceptable (albeit only marginally) fits with Y = 0.26 at an age of 9.5 Gyr; [α/Fe] = 0.2 at an age of 9.7 Gyr, and [Fe/H] = −0.90 at an age of 11.3 Gyr. This suggests that V41 may have formed from material enriched in products of stellar nucleosynthesis, and thus it may indeed be younger than V40. We note parenthetically that the results of our fitting experiments are in line with the findings of Kaluzny et al. (2013), who, while discussing (M−R) diagrams for binaries in M4, wrote "increasing Y decreases the age derived from the (M−R) plane while increasing [Fe/H] or [α/Fe] increases the age."

V40 is also somewhat peculiar because of the eccentricity of its orbit. At an age of ∼12 Gyr and a period of 5.3 days the orbit should long ago have been fully circularized by tidal friction (Mathieu et al. 2004; Mazeh 2008). There is no evidence for the presence of a third body in this system, which means that it must have undergone a close encounter during the last few Gyr. According to Dalessandro et al. (2014) NGC 6362 may have had a complicated dynamical history, losing perhaps as much as 80% of its original mass. Since close encounters account for a few percent of the total mass loss (M. Giersz 2015, private communication), the very existence of the nonzero eccentricity of V40 lends some support to this conjecture.

While V40 and V41 seem to belong to distinct populations, we feel it is too early to claim a direct detection of non-coevality among stars in GCs. It is clear, however, that DEBs offer the potential to test theories describing the origin of multiple populations in GCs, and to verify estimates of age differences between populations based on CMD-fitting. For the latter, both the theoretical predictions and fitting procedures yield values ranging from ∼10⁷ years to several hundred Myr (Nardiello et al. 2015, and references therein), and any independent constraints are highly desirable. The most obvious test for the reality of the age difference between our two systems would be to determine the chemical composition of their components using disentangling software (see e.g., Hadrava 2009). The existing spectra do not have adequate quality to measure abundances, but, given the brightnesses of the components and the orbital periods, it should be possible to improve the S/N.

There is also room for improving the accuracy of parameter measurements. As noted by Kaluzny et al. (2013), the errors in the masses of the components originate almost entirely from the orbital solution, which may only be improved by taking additional spectra (preferably with the same instrument). Admittedly, this would require a substantial observational effort, as doubling of the present set of radial velocity data would lead to an improvement of only 33% in the mass estimates (Thompson et al. 2010). The errors in the radii of the components are dominated in turn by the photometric solution, whose quality depends on the accuracy of the photometry. We estimate that an accuracy of 0.002–0.003 mag in V would reduce the errors of the radii by 50%. This goal could be easily achieved on a 6–8 m class telescope equipped with a camera capable of good PSF sampling. Finally, IR photometry should be obtained to reduce errors in luminosity with the help of existing accurate calibrations linking (V − K) color to V-band surface brightness.

This series of papers is dedicated to the memory of Bohdan Paczyński. J. K. and M. R. were supported by the grant DEC-2012/05/B/ST9/03931 from the Polish National Science Center. A. D. received support from the Australian Research Council under grant FL110100012.

Kwee & van Woerden (1956) proposed an analyis of the timing of eclipses by fitting them with symmetric curves, so that any asymmetries would be reflected in observed-minus-calculated (O−C) deviations. Their classical method does have its drawbacks. As far as the algorithm is concerned: (i) the interpolation of a light curve into a somewhat arbitrary evenly spaced time mesh yields a slight dependence of a final estimate of the central time on its initial approximation and (ii) the resulting estimate of the eclipse center depends on the initial selection of data points and on the initial guess of the eclipse center. There are also statistical drawbacks: (iii) interpolation effectively assigns uneven weights to points and (iv) the effective model, a polyline, suffers from excess number of parameters describing all nodes. Such a waste of degrees of freedom is known in statistics to decrease accuracy.

Let us adopt the zero point of the time scale as the current estimate of the eclipse center T_c⁽ⁿ⁾. Following Kwee & van Woerden (1956) we fit both the input light curve and its reflexion with respect to T_c⁽ⁿ⁾. The novelty here is that we employ the orthogonal polynomial approximation of the combined light curve (e.g., Ralston 1965), rather than on a poly line interpolation. In doing so we optimize the sensitivity of our method by keeping the number of model parameters at a statistically justifiable minimum (Occam razor). By virtue of the Fisher lemma (e.g., Brandt 1975), parameters of the orthogonal model are uncorrelated with noise hence they yield the minimum variance model estimate.

Next, again following Kwee & van Woerden (1956), we freeze the model curve and shift T_c⁽ⁿ⁾ by ±δ and calculate χ² of O−C for the time shifts $-\delta ,\;0,+\delta .$ These three values of χ² are fitted with a parabola

$$\begin{eqnarray}&&{\chi }^{2}(t)={{at}}^{2}+{bt}+c\end{eqnarray} \tag{ 3 }$$

and the location of its minimum ${T}_{{\rm{c}}}^{(n+1)}=-b/2a$ may be adopted as the next approximation of ${T}_{{\rm{c}}}^{(n+1)}.$ However, to assure convergence, we trim time steps ${{dT}}_{{\rm{c}}}={T}_{{\rm{c}}}^{(n+1)}-{T}_{{\rm{c}}}^{(n)}$ and carefully select δ to avoid going too far away from the solution, yet possibly ensuring bracketing it within ±δ. Such orthogonal polynomials retain only even powers and are optimal both in terms of degrees of freedom and numerical procedure relying on orthogonal projection. By the incorporation of weight coefficients in the scalar product our algorithm is able to account for uneven errors of data.

Complications, both conceptual and algorithmic, arise due to our ultimate aim of minimizing the dependence of the final result on somewhat accidental selection of the initial data points. Discrete selection/modification of observations may result in numerical instability and/or non-converging limit cycles. To minimize such danger we resort to smooth down-weighting of points near edge of eclipse by a symmetric Fermi–Dirac distribution function (e.g., Kittel 1996):

$$\begin{eqnarray}&&f(t)=\displaystyle \frac{1}{{e}^{\frac{| t| -{\rm{\Delta }}}{\delta }}+1}.\end{eqnarray} \tag{ 4 }$$

The final weight of data measured at time t_i with an error σ_i is ${w}_{i}=f({t}_{i})/{\sigma }_{i}^{2}.$ As our time zero-point moved during iterations toward the final solution, so our weights vary between iterations. We roughly followed principles of the minimization algorithm by simulated annealing (Otten & van Ginneken 1989). At the start of iterations we adopt a rather crude model curve for m = 2 and smooth weighting δ = Δ/10. Note, that in statistical physics δ would be proportional to the temperature. Only after successful progress of our algorithm do we tighten our data sample by decreasing δ (i.e., "lower temperature") and refine our model by increasing m. Also, we trimmed the time step dT_c ≤ Δ/10.

After each iteration we calculate two T_c error estimates: Σ and σ. The error Σ corresponds to such a model shift δ = Σ that χ² grows over its minimum value ${\chi }_{\mathrm{min}}^{2}$ by its variance $\sqrt{2d},$ where $d=n-m/2-1$ is the number of degrees of freedom. The error σ is defined similarly, for growth by 1. For this purpose χ² has to be renormalized by a factor C such, that $E\{C{\chi }^{2}\}=d.$ For an approximate normalization we substitute the minimum value of the parabola in Equation (3) for its expectation, so that $C=(4{ac}-{b}^{2})/4{ad}.$ Σ and σ are calculated by solving of Equation (3). In this way we obtain: $a{{\rm{\Sigma }}}^{2}=\sqrt{2d}{\chi }_{\mathrm{min}}^{2}/d$ and similarly for σ²:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}^{2}=\sqrt{2}\;\displaystyle \frac{4{ac}-{b}^{2}}{4{a}^{2}\sqrt{d}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }^{2}=\displaystyle \frac{4{ac}-{b}^{2}}{4{a}^{2}d}.\end{eqnarray} \tag{ 6 }$$

To force a positive result we apply Equation (3) to $\mathrm{log}\;{\chi }^{2}$ rather than to χ². In statistical terms, Σ corresponds to the maximum reasonable uncertainty of T_c, were other parameters strongly correlated with it and σ is its likely error, assuming weak correlation.

During iterations we keep δ as large as δ = 2Σ. Thus δ decreases for an improved model fit at a moderate rate. The polynomial order is raised in steps of 2, up to its maximum allowed value, only after decrease of corrections of T_c to ${{dT}}_{{\rm{c}}}\equiv {T}_{{\rm{c}}}^{(n+1)}-{T}_{{\rm{c}}}^{(n)}\lt 0.01{\rm{\Sigma }}.$ We terminate iterations when m reaches its selected maximum value and time steps become small ${{dT}}_{{\rm{c}}}\lt \epsilon {\rm{\Delta }},$ where is our tolerance parameter. Our result is the final eclipse center T_c and its likely error σ. Note, that for correlated observations errors need to be increased (Schwarzenberg-Czerny 1991).