© ESO, 2017

1. Introduction

HD 64315 (HIP 38430, CD − 26°5115, V402 Pup) is the main ionising source of the Galactic H ii region Sh2-311 (Sharpless 1959). This bright nebulosity (and hence the star itself) is apparently connected to a number of dark and bright clouds, extending over almost 1° on the sky in the region of Puppis, which is frequently referred to as NGC 2467. Originally, NGC 2467 was believed to be a large cluster, but several authors have concluded that it is actually the projection of a number of bright foreground B- and A-type stars on top of a distant star-forming association (Lodén 1966; Fenstein & Vázquez 1989). Within these clouds, two compact young open clusters, Haffner 18 and Haffner 19, were also found to be illuminating smaller H ii regions. Haffner 18 contains an O7 V star and probably three O9 stars (Moreno-Corral et al. 2005). Haffner 19 only contains stars up to B0–B1 (Munari & Carraro 1996). Examination of wide field images of the area very strongly conveys the impression of a single star-forming region, including the two clusters and NGC 2467 (now understood only as the area surrounding HD 64315), an idea also supported by analysis of Spitzer observations of the area (Snider et al. 2009). Several authors have investigated this hypothesis by deriving distances to the three clusters, obtaining discrepant results. FitzGerald & Moffat (1974) placed Haffner 18 and Haffner 19 at 6.9 kpc, a much larger distance than the 3.7−4.4 kpc calculated for HD 64315 and its associated H ii region (Georgelin & Georgelin 1970; Cruz-González et al. 1974; Pismis & Moreno 1976). More recent studies suggest that Haffner 18 is more distant than Haffner 19 (Munari & Carraro 1996; Munari et al. 1998; Moreno-Corral et al. 2002, 2005; Yadav et al. 2015), the latter having a distance more compatible, but still larger than NGC 2467.

In the case of HD 64315, the few distance determinations found in the literature (Cruz-González et al. 1974; Pismis & Moreno 1976) are based on photometry and the assumption that it is a single star. However, this star, initially classified as O6Vn (Walborn 1982), was found to be a double-lined spectroscopic binary by Solivella & Niemela (1986). By using medium dispersion spectrograms obtained at CTIO between 1982 and 1984, Solivella & Niemela identified the two components as ~O6 stars and derived preliminary orbital elements for the binary system, obtaining a period of 1.34 days and a mass ratio of ~0.83. In view of this binary nature, all the studies mentioned above must have underestimated the distance to the star.

With the aim of investigating the orbital and stellar properties of this binary system, we obtained extensive high-resolution, high signal-to-noise ratio (S/N) spectra of HD 64315 with the Fiber Extended Range Optical (FEROS; Kaufer et al. 1999) and the FIbre-fed Echelle Spectrograph (FIES; Telting et al. 2014) spectrographs. In the analysis process, we found strong signatures of more than two components present in the spectra, which complicated the spectroscopic analysis of the star, but allowed us to discover (in combination with other photometric and spatial information) that HD 64315 is in fact a multiple system comprising at least four components.

In this paper, we present observational evidence of the multiple nature of HD 64315 and its implications for the distance determination to this stellar system. The paper is structured as follows. The spectroscopic observations and the photometric data are presented in Sect. 2. The spectra and their spectral classification are discussed in Sect. 3. A preliminary discussion on the distance to the source is made in Sect. 4, based on the interstellar lines present in the spectra. We then present evidence in Sect. 5 showing that HD 64315 is composed of two spectroscopic binaries. In Sect. 6 we describe the procedure to extract the radial velocities of each component from the spectra and develop a comprehensive analysis to obtain the orbital and stellar parameters for both binary systems. We conclude with the direct estimation of the distance in Sect. 7. The discussion of results is presented in Sect. 8, and our main conclusions are presented in Sect. 9.

2. Observations

A total of 104 spectra grouped in 25 observing blocks (OBs) were obtained in service mode at random phases between 2006 October and 2007 March with the FEROS instrument at the ESO/MPG 2.2 m telescope¹ (see Table 1; spectra numbered from 1 to 25). This first set of spectra was complemented with seven further OBs (14 spectra; see Table 1; spectra numbered from 26 to 32) obtained by CJE observations on 5 nights in 2009 March. This time, they were observed at specific phases, with some of the OBs separated by only a few hours. The third set of spectra was obtained in the framework of the IACOB project (Simón-Díaz et al. 2015) with the high-resolution FIbre-fed Echelle Spectrograph (FIES) attached to the Nordic Optical Telescope (NOT), located at the Observatorio del Roque de Los Muchachos (La Palma, Spain) between 2013 January 29 and 31 (11 spectra; see Table 1; spectra numbered from 33 to 43). Finally, the last group (see Table 1; nine spectra numbered from 44 to 52) was also taken with FEROS by DM during a long run in 2014 May. All the FEROS spectra were reduced using the reduction pipeline that runs under the MIDAS environment (Kaufer et al. 1999). The spectra from each observing block were combined in order to have a higher signal-to-noise ratio (S/N) and eliminate possible cosmic ray contamination. Those taken with FIES were homogeneously reduced using the FIEStool² software in advanced mode. A complete set of bias, flat, and arc frames obtained on each night were used to this end. For wavelength calibration, we used arc spectra of a ThAr lamp. In most cases the S/N of the target spectra is in excess of 70 per resolution element; half of them have a S/N above 100. The spectra were normalised and heliocentric corrections were applied using our own code developed in IDL. The final set of spectra is summarised in Table 1, where we also show the exposure time and S/N for every spectrum.

The spectroscopic observations were complemented with photometric data from the All Sky Automated Survey³ (ASAS Pojmanski 2003) and observations taken by the Hipparcos satellite in the H_p band. There were a total of 544 photometric datapoints extracted from the All Sky Automated Survey catalogue (all in the V band). We chose the photometric data corresponding to aperture MAG 2, which has the smallest intrinsic error (σ = 0.034 mag). They are displayed in Fig. 1 (top panel). The H_p passband embraces the V and B passbands. The transformations from H_p to Johnson filters are comprehensively described in Harmanec (1998). The number of photometric points extracted from the Hipparcos catalogue is 149 for every passband (shown in both panels of Fig. 1). In total, we gathered 693 points in the V filter and 149 in the B filter. We observe the same photometric variability in the data from the ASAS and Hipparcos catalogues.

Fig. 1 Photometric data extracted from the All Sky Automated Survey (black dots) and Hipparcos catalogues (red dots). The V passband is shown in the top panel. The lower panel shows the B-band data.

3. Description of the spectra

Figure 2 illustrates four representative examples of the spectrum of HD 64315 at different phases. We plot the wavelength range between 4300 Å and 4710 Å, where the main lines used to define the spectral type (SpT) in mid-O stars are found (see e.g. Walborn & Fitzpatrick 1990). The first spectrum from the bottom (#13) is probably similar to the one analysed by Walborn (1982). Following the Morgan-Keenan (MK) system of spectral classification this spectrum can be classified as ~O6 Vn (where the suffix n indicates that the lines are broad). Continuing to the top of the figure, the third spectrum (#4) can be also classified as O6 V; however this time, the lines are narrower. In addition, the He ii 4686 line is now stronger than the He i 4471 Å and He ii 4542 Å lines. Spectra with this type of morphology have been given the “Vz” qualifier (see e.g. Walborn 2007), and hypothesised to correspond to lower (visual) luminosity and younger ages (but see discussion by Sabín-Sanjulián et al. 2014). In the two other spectra (#16 and #24), it becomes clear that the star is a double-lined spectroscopic binary (as seen most clearly in spectrum #24 in the He i 4471 Å line). Previous spectral classifications are thus the consequence of the morphological analysis of spectra obtained at phases in which the two components are blended together.

More recently, Arias et al. (2016), accepting the multiplicity of the system, classified HD 64315 as O5.5 V+O7 V, removing the system from the “Vz” category. At classification resolution, only two components are identified. From a first inspection of some of the spectra with maximum separation between both components (see e.g. spectra #16 and #24 in Fig. 2), the spectral types could be estimated as O6 and ~ O5−O5.5 for the components with stronger and fainter lines, respectively. Curiously, the component with fainter lines seems to have an earlier spectral type. In addition, the broadening of the lines from this component is noticeably different in some of the spectra with a similar separation between the two components. Indeed, a closer inspection of the global spectrum at certain phases already gives strong hints of the presence of more than two components, all moving rather quickly in radial velocity.

Fig. 2 Representative examples of the spectrum of HD 64315 at four different phases. Spectra are numbered for ease of identification.

4. Distance determination using interstellar lines

Before presenting the observational evidence for more than two components in HD 64315 and its physical characterisation using the available spectroscopic and photometric datasets, we discuss the kinematic distance to HD 64315. Using the interstellar Na i D lines (5890.0 Å, 5895.9 Å), we studied the radial velocity distribution of the interstellar material in the direction towards HD 64315 (; ). We calculated the velocity scale with respect to the local standard of rest (LSR) by assuming that the Sun’s motion with respect to the LSR corresponds to +16.6 km s^-1 towards Galactic coordinates l = 53°; b = + 25°. The interstellar Na i D lines are shown in Fig. 3 as normalised flux as a function of LSR velocities. The two lines have a very similar shape, showing two distinct components. Both components display only positive velocities from +1 to + 74 km s^-1. None of the components is saturated, and so we can determine their centres. The broader feature is centred at + 21 km s^-1, while the narrower feature is at + 58 km s^-1.

There are few stars along similar lines of sight with well-studied interstellar lines. HD 68761, with galactic coordinates ; , displays only one component in its interstellar lines, centred at + 8 km s^-1. This star is situated at a distance ≲ 1.5 kpc (Hunter et al. 2006). HD 58978 (; ), has several distance estimates, ranging from ~ 800 pc to ~ 1.7 kpc (Hunter et al. 2006). Its interstellar lines display a main component centred at + 11 km s^-1 and a very weak component with an edge velocity around + 29 km s^-1.

Gyulbudaghian & Akopian (2002) measured the velocities of all catalogued molecular clouds in this area, finding that all had velocities between +20 and +25 km s^-1 and typical distances ≲ 1 kpc. All this suggests that the main component in the interstellar lines of HD 64315 arises from relatively nearby clouds. The clear separation of the second component suggests that the extinction is very low at intermediate distances. The second component shows velocities corresponding to much higher distances.

Fig. 3 Interstellar lines in the spectrum of HD 64315. The main panel shows the components of the Na i D doublet (5889.95 Å solid line; 5895.92 Å dashed line) in velocity space. The inset shows the Galactic rotation curve along this line of sight. The velocities are measured with respect to the local standard of rest (LSR) assuming a solar motion of + 16.6 km s-1 towards Galactic coordinates ℓ = 53°; b = + 25°.

Pismis & Moreno (1976) found kinematic distances to the H ii regions around HD 64315, Haffner 18 and Haffner 19. Their observations were re-reduced by Moreno-Corral et al. (2002), who find LSR values of +50, +58 and +50 km s^-1, respectively. However, Pismis & Moreno (1976) warn that the H ii velocity close to HD 64315 seems to be affected by the expansion of gas around the ionising star, with smaller values of v_LSR in the immediate vicinity of the star. This interpretation is borne out by the presence of interstellar material with higher v_LSR along the line of sight. The data available can be interpreted as suggesting that the actual radial velocity of the complex is at least + 58 km s^-1, with the lower values measured at some points due to expanding shells around the ionising stars. This agrees very well with the shape of the weaker component of the interstellar lines, centred at v_LSR ≈ + 58 km s^-1. The inset in Fig. 3 shows the Galactic rotation curve in this direction, computed considering circular galactic rotation and adopting the rotation curve of Reid et al. (2014). Along the line of sight towards HD 64315, all radial velocity curves display only positive and monotonically increasing values, from a distance around 0.25 kpc. Assuming v_LSR ≈ + 58 km s^-1, the kinematic distance estimate to HD 64315 corresponds to d ≈ 5 kpc, which we adopt as a preliminary distance.

5. Evidence for more than two components

5.1. Period determination

We carried out a timing analysis of the photometric data described in Sect. 2. We used the period program inside the Starlink suite for every passband. The Lomb-Scargle algorithm (Lomb 1976; Scargle 1982), in a range of frequencies 0−35.7 d^-1, gives a photometric period of 1.0189650 ± 0.000112 d and 1.018958 ± 0.000018 d for the B- and V-filter data, respectively. The agreement between the two periods is excellent. These values are confirmed with the clean algorithm (Roberts et al. 1987), which removes spurious periods caused by the window function. The Lomb-Scargle periodogram in a range of frequencies 0−10 d^-1 and an inset of the clean periodogram, up to ν = 5 d^-1 are shown in Fig. 4. There is clearly only one peak above significance.

Fig. 4 Results of the Lomb-Scargle periodogram for the photometric data in the V filter (left). The result of applying the clean algorithm is shown (right). The frequency peak corresponds to half the orbital period.

Figure 5 shows all the photometric data for both passbands folded on the derived period. The corresponding light curves in the V and B filters show an amplitude of ≈ 0.18 mag and ≈ 0.05 mag, respectively. Error bars of the data are also displayed in Fig. 5, showing that most of the errors exceed the amplitude of the light curves (probably due to the process of transformation to the standard system). The shallow amplitudes and short-period observed might suggest that the variability is due to ellipsoidal light variations (Wilson & Sofia 1976). However, as we shall see, the light from the system is not dominated by the binary producing the photometric variability, and the variations represent an eclipse.

Fig. 5 Photometric light curves (V filter, top; B filter, bottom), including data and the corresponding error bars from ASAS (blue open dots) and Hipparcos (black solid dots).

Fig. 6 Comparison of several pairs of spectra taken on different nights, but corresponding to the same phase (blue numbers) of the 1.019 d photometric period. The presence of more than two components is obvious.

5.2. Spectroscopy and radial velocities

Solivella & Niemela (1986) determined a spectroscopic period of 1.34 d, based on around 20 medium dispersion spectrograms, presenting a large dispersion in radial velocities with respect to the model curve fitted. The difficulty in obtaining a good fit is evidenced by a difference of 40 km s^-1 between the systemic velocities of the two components. A more recent attempt to characterise HD 64315 as an eclipsing binary was carried out by our group (Lorenzo et al. 2010). We used the first 24 spectra described in Table 1 and obtained a binary period of 2.71 d (twice the period of Solivella & Niemela 1986), with the same systemic velocity for the two stars. Even so, the radial velocity curves showed an unacceptable standard deviation of 35 km s^-1; the residuals of the radial velocities increase at phases around zero and 0.5 (i.e. near the eclipses). Given the total disagreement between the photometric and spectroscopic period and the hints of the presence of more stellar components, we decided to schedule further observations so that some spectra were taken at the same phases with respect to the photometric period (1.01896 d), while other spectra were taken separated by just a couple of hours, a very small phase shift with respect to the spectroscopic period. With these constraints, we expected to reveal changes in the morphology (line profiles), while keeping the components of the eclipsing binary at the same radial velocity.

As a consequence of this successful strategy, observations between 2009 March and 2014 May have been much more useful to understand how many components the system has and how each one of contributes to the combined spectrum. As examples, in Fig. 6, we show a set of pairs of spectra, represented by the He i 4471 and He ii 4542 lines, which define spectral type for O-type stars. Every pair was observed at a very similar phase (φ), according to the photometric period, but they were obtained on a different date. For instance, the spectra #30 and #47 (top panel of Fig. 6) were acquired with a time difference of more than 5 years, but at phase φ ≈ 0.91 (according the period of the eclipsing binary, 1.019 d). First of all, they do not show the same morphology. Moreover, the line morphology is not typical of a binary close to eclipse, where the radial velocity of both components should be close to systemic velocity. Rather He ii 4542 shows a double line in spectrum #30. A comparable situation happens for the pairs of spectra #7, #51 (φ = 0.87) and #9, #52 (φ = 0.83). However, the spectra #37, #8 (φ = 0.81) have a time difference of more than 6 years, but their profiles are morphologically very similar.

We display a few more examples in Fig. 6 to emphasise the complexity of this stellar system. The only sensible conclusion after all these observations is that the stars forming the spectroscopic binary (or, at least dominating the spectrum) are not the same ones giving the photometric signal at 1.019 d. The 2.71 d period found from the radial velocity analysis (Lorenzo et al. 2010) corresponds to a different binary.

Once we have established that there is a non-eclipsing SB2 binary giving the prevailing spectroscopic period of 2.71 days and an eclipsing binary causing the photometric 1.019 d signal, we can guess that the components of the SB2 system should be more luminous than the stars in the eclipsing binary system, as they dominate the combined flux. Since two stars with spectral types not very far from O6 are seen, the spectral types of the components of the eclipsing binary must be later than O7 V, so that they are later than the secondary star in the non-eclipsing binary system.

5.3. Interferometry and parallax

Using Speckle interferometry obtained over a period of a decade, Mason et al. (2009) derived the angular separation of HD 64315, resolving at least two components with a separation of . In Sect. 4, we discussed the distance to HD 64315 and estimated d ≈ 5 kpc. At this distance, the geometric distance between the components would be 455 AU (≈ 100 000 R_⊙). This separation is certainly not consistent with the short orbital periods that we have found. Furthermore, Tokovinin et al. (2010) obtained Speckle interferometry images of HD 64315 (see Fig. 12 in the mentioned study) suggesting that the system may consist of three visual components in a linear configuration (Tokovinin et al. 2010). The fit of the x-axis scan with two components is marginal, but experiments with three component made little or no improvement (Aldoretta et al. 2015). Therefore the two visual components resolved represent two systems of comparable (but not equal) brightness that are either in a very wide orbit or are not bound. The simplest possibility is that each the visual components represents one of the two binaries that we have identified (one in the photometric signal and the other in the spectra). Hereafter, we name the non-eclipsing binary system A, composed of stars Aa and Ab, while the eclipsing binary is known as B, with components Ba and Bb. Both are spectroscopic binaries (SB2), but the lines corresponding to Ba and Bb are mostly hidden within the complex and broad profiles generated by the more luminous components of system A. In the following section (Sect. 6), we describe the technique developed to disentangle the spectra and the subsequent orbital analysis.

6. HD 64315, a double spectroscopic binary

6.1. Determination of radial velocities

In Sect. 5, we presented inescapable evidence that HD 64315 is a double spectroscopic binary, with the non-eclipsing binary (nEB = A) composed of two mid-O stars and the eclipsing binary (EB = B) containing two stars of lower temperature. The orbital period of B is obviously the photometric period derived in Sect. 5.1, while the orbital period of A must be around 2.71 d. Our spectroscopic monitoring (see Table 1) covers almost8 years, during which time A has completed more than 1000 orbital cycles.

Fig. 7 Final fit: the observed lines (blue), the sum of four Gaussian functions (black), four Gaussian functions in every spectral line (Aa: blue; Ab: red; Ba: green; Bb: magenta). Numbers represent the order of the spectrum (for every pair: left He i and right He ii).

Every one of the four components contributes to every spectrum with its flux. Since we see them in the He i, He ii and Balmer lines, all components are O-type stars. According to the classical criteria for O-type spectral classification by Walborn & Fitzpatrick (1990), the ratio between the He i 4471 Å and He ii 4542 Å lines is sensitive to temperature and thus we chose these lines to fit the spectral line shape with multiple functions. The spectral line of a single star can be approximated with a Gaussian function, characterised by three parameters: amplitude, width and centroid of the peak. Even though line shapes in isolated stars present more complex profiles, the many added complications that are discussed below call for the choice of a simple Gaussian shape. The amplitude and width are dependent on the spectral line, and the centroid is dependent on the time when the spectrum was taken (i.e. the radial velocity of the component at that time). Since we have analysed two spectral lines in 52 spectra and we require four Gaussian functions to fit each line profile, our model must have eight widths, eight amplitudes, and 52 positions of the centre of the peaks. Our unknowns are ten orbital parameters: period, zero point of ephemeris, systemic velocity, and semi-amplitude of the velocity curve for each binary (P_A, T_0,A, v_0,A, K_Aa, K_Ab, P_B, T_0,B, v_0,B, K_Ba, K_Bb). From the light curves, two of these parameters, P_B and T_0,B, are well determined. This means that we know the phases corresponding to the orbital period of EB for every spectrum. Moreover, given the short orbital periods and large stellar sizes, we can assume that both binary systems are circularised and synchronised, and so both eccentricities are zero. To derive the radial velocities, we solve the inverse problem: we vary the free orbital parameters until the sum of four Gaussian functions matches the shape of the line observed. The minimization function is constrained by the expression for the radial velocity due to the orbital motion derived from the Kepler laws (see Eq. (2.45) in Hilditch 2001). From an analytical point of view, if we consider that f_iα(v) is the normalised flux of a spectral line i in the spectrum numbered α, with a dependence on velocity v, and g_iαX(v) is the Gaussian function which represents the flux distribution for star X (Aa, Ab, Ba, Bb) of the spectral line i (where i can be one of two spectral lines; He i 4471 Å  or He ii 4542 Å) in the spectrum α, and finally M_iα(v) is the sum of the four Gaussian functions, representing the overall line, given by (1)the function to minimise F(v) would be equivalent to (2)The function F(v) was implemented via the Python interface called lmfit. This package builds complex fitting models for non-linear least-squares problems. Spectra were transformed to the velocity space, in a range from − 600 km s^-1 to + 500 km s^-1, and rebinned to 200 bins. The remaining free orbital parameters, P_A, T_0,A, v_0,A, K_Aa, K_Ab, v_0,B, K_Ba and K_Bb, were bounded. The sixteen free parameters to characterise the widths and heights of every Gaussian function were constrained using our knowledge of the system: given that P_A is longer than P_B and the orbits are synchronised, the rotational velocities of Ba and Bb will be higher than the rotational velocities of Aa and Ab, and so we expect that Gaussians corresponding to Ba and Bb are wider than those for Aa and Ab. We consider that the flux of the stars is constant during an orbital cycle and that the stars do not pulsate. Finally, we do not take into account the Struve-Sahade (Struve et al. 1958) and Rossiter-Mclaughlin (Rossiter 1924; McLaughlin 1924) effects.

The resulting models (104 in total = 2 spectral lines × 52 spectra) are displayed in Fig. 7. In each panel, we can see the four Gaussian functions, representing the contribution of every component to the spectral line. The positions of the peaks are the radial velocities of the different components. To evaluate the goodness of the fit, we present the Pearson product-moment correlation coefficient R² for every fit in Table 2. The coefficient R² is > 0.9 in ≈ 42% of the models for the He i line, and just 8% of them have R²< 0.7. In the case of the He ii line, ≈ 62% of the models have R²> 0.9 and only 2% present R²< 0.7. We did not find any correlation between the best fits and the orbital phases, or the observing campaigns. Radial velocities are shown in Table 3. In the next section, we use these four sets of radial velocities to derive orbital and stellar parameters.

6.2. Orbital analysis

6.2.1. Non-eclipsing binary (AaAb)

We determined the orbital parameters from the radial velocities derived (Fig. 7). We used the Spectroscopic Binary Orbit Program (SBOP) by Etzel (2004). The guess parameters are those obtained by the fit described in Sect. 6.1. The radial velocity curve and residuals are displayed in Fig. 8. The orbit is circular, as this was an initial assumption. In the case of system A, we cannot determine the absolute parameters, as we do not know the inclination of the system with respect to the line of sight, but we can obtain a minimum mass for every component. The spectroscopic period determined is very similar to that obtained in Lorenzo et al. (2010), but now the uncertainty is very small, and the residuals are very low, less than 1 km s^-1. These are not the real uncertainties that can be derived from the observations, because, as explained in Sect. 6.1, we imposed a constraint on the model to determine the radial velocities. Radial velocities ordered by phase are shown in Table 4. The parameters derived are shown in Table 5. The mass ratio is 0.9366, very close to unity. The systemic velocity is the same for both radial velocity curves. When transformed to the LSR, it is + 55 km s^-1, in very good agreement with the velocity of the H ii nebula within which HD 64315 is embedded (Sect. 4).

Fig. 8 Radial velocity curves for the non-eclipsing binary (A), fitted to the observational data and shown against orbital phase (blue line: star Aa; red line: star Ab). The red dotted line corresponds to the systemic velocity. The residuals are shown in the top panel.

Recently, Arias et al. (2016) classified HD 64315 as O5.5 V +O7 V. This classification does not take into account the contribution to the spectrum of the eclipsing binary (Ba and Bb). Our results suggest that Aa and Ab are indeed quite similar in mass and temperature. If we calculate the sizes of their Roche lobes, following Eggleton (1983), we find R_⊙  and R_⊙. We can estimate the projected rotational velocity by assuming that the two stars just fill their lobe radii. In this case, it would be ~ 150 km s^-1 for both components. As the stars may not fill their Roche lobes, their actual rotational velocity will be higher, in agreement with the expectation of synchronisation. The lack of eclipse implies a lower limit on the orbital inclination. If the stars fill their Roche lobes, it will be around 45°. For this inclination, the masses of the two components are 31 M_⊙  and 29 M_⊙. As these masses are consistent with calibrations for a spectral type O6 V (Martins et al. 2005), we do not expect the inclination to be much lower. This strongly suggests that all the assumptions are approximately correct. The semi-major axis of the binary system would then be 32 R_⊙. Given these orbital parameters, this system is not a contact binary (see Fig. 9), but a detached binary system. However, given estimates for synchronisation times in stars with radiative envelopes (Zahn 1975; Claret & Cunha 1997), the assumption of synchronisation is fully justified, and so our assumption that the lines of both components of A are narrower than those of Ba and Bb is correct. As mentioned in Sect. 3, the component with fainter lines seems to have an earlier spectral type in some spectra. With our solution, this can be understood as follows: the components both have approximately the same spectral type, around O6 V, but they appear slightly different in different spectra because changes in the width and height of spectral lines are caused by the many confusing effects inherent to this kind of early-type close binary (such as Struve-Sahade or tidal distortion) with the added complication of the contribution of the eclipsing binary sibling.

Fig. 9 Representative drawing of HD 64315 AaAb at different phases created with the phoebe 2.0-alpha code via the Python interface. The secondary is slightly greyer and smaller to differentiate it from the primary.

Fig. 10 Radial velocity curves fitted to the observational data and shown against orbital phase (blue line: Ba; red line: Bb). The red dotted line corresponds to the systemic velocity. The residuals are shown in the top panel.

6.2.2. Eclipsing binary (BaBb)

We derived the radial velocity curve (see Fig. 10) and the corresponding orbital parameters, by using SBOP and the same methodology as in Sect. 6.2.1. The mass ratio, semi-major axis and systemic velocities obtained were included as guess parameters in the generalised Wilson-Devinney (WD) code (Wilson & Devinney 1971) in its 2010 version. The period and zero point ephemeris derived from the spectroscopic data are 1.018965 d and HJD 2 452 550.6272, respectively. Radial velocities ordered according to the phase are shown in Table 6. The tiny residuals are a consequence of the same procedure discussed in Sect. 6.2.1.

We obtained the light curve and radial velocity curve models by computing the parameters via differential corrections until all free parameter adjustment of light curves and radial velocity curves is reached according to the least-squares criterion. We assumed that stars Ba and Bb form a contact system, in which case circularisation and synchronisation are acceptable approximations. We chose mode 1 of the WD code, corresponding to overcontact binaries. In this case, the surface potentials are the same for both stars (Ω₁ = Ω₂). The radiative model for both components of the binary system is an atmosphere model by Kurucz (1993). The surface is divided into a grid of 40 × 40 elements for each star. To improve the convergence of the solution, we chose symmetrical derivatives (Wilson & Sofia 1976). The code can apply the detailed reflection model of Wilson (1990), a treatment especially recommended with overcontact binaries. We have also considered proximity effects on both stars. A square root limb-darkening law was applied during the process, as it is an order of magnitude more precise than the linear law (van Hamme 1993). The bolometric albedos of both components were fixed at A₁ = A₂ = 1, because the atmospheres are expected to be in radiative equilibrium (von Zeipel 1924). Because of local energy conservation, this also implies gravity brightening exponents g₁ = g₂ = 1. Other constraints applied are described in mode 1 of the WD code. Both temperatures were fixed to 32 000 K, in agreement with the expectation of a spectral type not later than O9.5 V nor earlier than O8.5 V (as both components have to be O-type stars, and at the same time considerably fainter than the ~O6 V components that dominate the spectrum). The temperature ratio is not an adjustable parameter during the convergence process. This is a compulsory constraint due to the morphology of the binary star, where both components share a volume of their Roche lobes and thermal contact is assumed. The third light is included due to the presence of the binary system A.

Fig. 11 Light curve model fitted to the observational data and residuals in V filter (above) and B filter (below).

The process to convergence of all free parameters is iterative and simultaneous for all observables analysed, i.e. the radial velocity curves and light curves. The criterion for convergence adopted is as follows. For three consecutive iterations, all adjustable parameters must be within two standard deviations. Once convergence is reached, five solutions are derived by varying the parameters within the standard deviation and fitting the observations again. We choose the fit with the smallest dispersion as a final solution.

Light curve models are shown in Fig.11. Their shape phenomenologically corresponds to an eclipsing variable of the EW type; there is no plateau between the eclipses. The light curve exhibits a continuous and monotone shape along the cycle. This shape of the model light curve confirms that the two stars are overfilling and sharing their Roche lobes. Light curves show a significant dispersion. As a consequence, it is difficult to visually distinguish which one of the two minima is deeper. For a best estimate of the difference in depth, we averaged all photometric data points between φ = 0.99 and 0.01, and all points between φ = 0.49 and 0.51, obtaining a difference of 5 millimagnitudes for the V filter and 2 millimagnitudes for the B filter. The latter value is not significant because of the small number of points. Both differences are smaller than the intrinsic dispersion (σ_V = 0.02; σ_B = 0.006) of the corresponding light curves, and so we have to conclude that both minima are of equal depth. Residuals are under 0.08 mag for the V filter and under 0.02 mag for the B filter, showing a reasonably good fit in either case, in spite of the low quality of the photometric data.

The linear ephemeris equation, where the epoch of successive times of primary-eclipse minima (phase zero), T_min, is calculated from the period and zero-time ephemeris derived from the light curves, is All the stellar and orbital parameters are shown in Table 7. The difference between the spectroscopic and photometric period is about seven-tenths of one second, and the zero-time ephemeris derived from spectroscopic and photometric data differ by about 100 s, demonstrating excellent agreement, given the complex methodology and low quality of the data. The systemic velocity transformed to the LSR coordinates would be + 35 km s^-1, significantly different from the systemic velocity of the non-eclipsing binary (AaAb). Both components have the same mass, 14.6 M_⊙. In view of this, we assume that both components have spectral type O9.5 V. As is typical in this kind of early-type overcontact system, the surface effective gravity obtained is quite high, typical of zero-age main-sequence (ZAMS) stars.

The radii derived from the surface potential show that the two stars are sharing their atmosphere. We provide a representative drawing (see Fig. 12) where we can see how the atmosphere of the two stars overlap.

7. Direct distance estimation

Eclipsing binaries allow the derivation of geometrical distances to the systems, which can be very precise (e.g. Southworth et al. 2004; Vilardell et al. 2010). In overcontact binaries, there are many complications, owing to the interaction and geometrical distortion, but a direct distance estimate is still possible from the stellar parameters. In HD 64315, we find the added complication of a third body (the detached non-eclipsing binary) that is brighter than the eclipsing binary. For the calculation, we assumed effective temperatures of 32 000 K for components Ba and Bb, and 40 000 K for components Aa and Ab. This is an approximation based on qualitative criteria described in Sects. 6.2.1 and 6.2.2 and the SpT-T_eff calibrations by Martins et al. (2005). The apparent magnitudes of every star (for every passband) were calculated and included in Table 8. We assumed the uncertainties to be twice the standard deviation of the corresponding light curve. Then, we estimated a distance-dependent flux, in standard physical units, including those proximity effects supported by the model and accepting constraints relative to the radii of components of Aa and Ab (R_⊙  and R_⊙). This procedure, although approximate, is more accurate than the simple use of the mean radii. We follow the procedure described in Vilardell et al. (2010). The parameters used are shown in Table 8, together with our estimates of distance, 4.7 ± 0.6 kpc, for the non-eclipsing binary (AaAb) and 5.0 ± 0.8 kpc for the eclipsing binary (BaBb). Both values are fully consistent. The errors are mainly due to the relatively large uncertainties in the extinction and the bolometric correction.

Fig. 12 Representative drawing of HD 64315 BaBb to scale at quadrature phase, created with the phoebe 2.0-alpha code via the Python interface.

Fig. 13 Radial velocity curves of the four stars of HD 64315 in time (above; Aa: blue, Ab: red, Ba: green, Bb: magenta). The He ii line of each spectrum is shown (solid black line) compared to the model from all components (dashed line).

8. Discussion

HD 64315 presents evident difficulties in the determination of its orbital and stellar parameters. Firstly, all components have broad and shallow spectral lines, because of their high rotational velocities. As a consequence, cross-correlation techniques are difficult to apply because the peaks of the spectral lines for every component are not distinguishable, thereby increasing the uncertainties. Moreover, the short periods of both binaries add further complication to any disentangling process. As an example, in Fig. 13 we show radial velocity curves for all components as a function of time (HJD = 2 454 072−2 454 077). Spectra 8, 9, 10 and 11 are marked on the radial velocity curve with their corresponding times. The spectral lines of He ii corresponding to these four spectra are displayed together with the contribution of every component and the sum of these Gaussian functions. Apparently, radial velocities in spectra #8 and #11 are very similar quantitatively, but the morphology of the spectral lines is totally different: while #11 shows two peaks, #8 shows an erratic profile aggravated by its low S/N. In the case of spectra #9 and #10, we note the asymmetries of the spectral lines, and how a spectral line showing a single peak can hide four lines coming from the four components. The contributions of Ba and Bb to the profile are shallow and weak, and, given the higher orbital velocity of system B, on many occasions appear as very faint extended wings on the sides of the lines due to system A or, when their radial velocities are low, are completely lost inside the lines from A.

Despite these difficulties, our methodology is able to derive orbital and stellar parameters with a very limited set of assumptions. We find that HD 64315 is a quadruple star, consisting of two binary systems separated by about 100 000 R_⊙  (455 AU) at a distance of around 5 kpc. The non-eclipsing system (O6 V+O6 V) is a detached system, However, the inclination derived suggests that the two components are not far from filling their respective Roche lobes, and we have used this approximation for some estimates.

The eclipsing binary (O9.5 V+O9.5 V) is a contact system with a very short period. We consider this object a very strong merger progenitor candidate. To estimate when the merger will take place, we assume that the coalescence of the two components will occur when the outer Roche lobe radius is reached. Application of the equation derived by Yakut & Eggleton (2005) leads to an outer Roche lobe radius for HD 64315 Ba of 6.3 R_⊙. If we assume that the Kevin-Helmholtz time is the merging time scale, the duration of the merger will be 26 000 years. The two components of HD 64315 B are already overfilling the volume of their Roche lobe, presenting an overlapping volume between them. If we take, for comparison a Geneva track (Georgy et al. 2013) for a 15 M_⊙ with Z = 0.01 and Ω/Ω_crit = 0.568 (equivalent to v/v_crit = 0.4), we see that a polar radius ≈ 5 R_⊙ (as we find from our solution) corresponds to an age around 3 Myr (in good accord with the presence of an H ii region around the star). If the star was isolated and stellar evolution was the only driver of shape changes, an equatorial radius ≈ 6.3 R_⊙ would be reached by an age of ≲ 8 Myr. However, in the case of HD 64315 Ba, we see that the distortion from a spherical shape is much larger than caused by rotation alone. Alternatively, a decrease of 2000 K in its T_eff would also imply overfilling the outer Roche lobe. Again, the evolutionary timescale for an isolated star to achieve that stage is a few Myr.

In the case of HD 64315 Aa, the outer Roche lobe radius of HD 64315 Aa is 15.1 R_⊙, and so the stars are very far away from the merging condition, as expected, given the much longer period. For two stars of this size, the merging timescale would be around 10 000 years. In both cases, the merger happens in about of the lifetime of an individual star with the same spectral type.

The two spectroscopic binaries that we have observed must correspond to the two objects separated by speckle interferometry, as we do not see any stationary component in the spectra. There may be further components, but then they have to be fainter than binary B. With a separation of ≈ 455 AU, the two systems must orbit each other with a period of ~ 1000 yr. They are thus likely gravitationally bound.

8.1. Astrophysical context

It is difficult to find stellar objects with similar characteristics in the literature, as most of the quadruple system studied correspond to visual binaries or systems with much lower masses. A notable example is the quadruple system QZ Car (Mayer et al. 2001), composed of binaries A (O9.7 I+b2 v, P_A = 21 d) and B (O8 III+o9 v,P_B = 6 d). It presents several differences with HD 64315: the primaries in the two binaries within QZ Car are evolved stars, and the secondaries are rather less massive. Moreover, the common orbit has a period likely measured in decades rather than centuries.

During the last decade, several studies have found the multiple nature of systems previously believed to be binaries, for example, LY Aur (Mayer et al. 2013) or SZ Cam (Tamajo et al. 2012), the latter in the central dense region of the open cluster NGC 1502. We expect that many binary systems known to present spectral asymmetries without explanation or deviating radial velocities will be found to be multiple in the near future (cf. Rappaport et al. 2013; Lohr 2015), probably resulting in a need to review statistics about percentages of binary and multiple stars in young clusters.

The two binaries within HD 64315, on the other hand, present quite typical characteristics. A system very similar to HD 64315 A is DH Cep. It has a very similar orbital period of 2.8 d. The components are classified as O5.5 V and O6.5 V, and their masses, as derived by Hilditch et al. (1996), are around 33 M_⊙ for the primary and 30 M_⊙ for the secondary. As in HD 64315 A, DH Cep is not an eclipsing binary and the measured rotational velocities are around 150 km s^-1. Since DH Cep is not complicated by the presence of a second binary, Hilditch et al. (1996) can derive accurate radii, finding values ≈ 8.5 R_⊙ for both components, well within their Roche lobes.

HD 64315 B is an overcontact system, and one of the most massive ones known. Among systems containing O-type stars, only OGLE SMC-SC10 108086 has a shorter orbital period of 0.88 d (Hilditch et al. 2005). With masses of 14 and 17 M_⊙, this SMC system is quite similar to HD 64315 B. The more massive LMC system VFTS 352 has an orbital period of 1.12 d (Almeida et al. 2015). Another massive binary with a shorter period, GU Mon, contains later-type components (B1 V+B1 V) with lower masses of around 9 M_⊙ (Lorenzo et al. 2016). All these objects have been considered likely merger binaries. HD 64315 B is exceptional among them, because of its membership in a multiple system. Its eventual merger will lead to the formation of a hierarchical triple system where all the components have about equal masses.

8.2. Formation and environment

We present different estimates of the distance to HD 64315, all agreeing on ≈ 5 kpc. Previous studies had not considered the multiple nature of the system and had obtained distances ≲ 4 kpc, based on calibrations of the luminosity for the spectral type assumed. From the interstellar lines seen in our spectra, we have identified a local component with a kinematic distance of ~ 5 kpc. The direct distance estimation has also given distances of 4.7 ± 0.6 kpc and 5.0 ± 0.8 kpc for the two binary systems. All measurements are clearly compatible. This new value brings HD 64315 to the same distance as most recent estimates for Haffner 18 (Yadav et al. 2015, and references therein). However, several recent papers point to a rather higher distance for Haffner 19, which is visually closer to HD 64315 (only away) and, unlike Haffner 18, contains an O-type star. We must look with some scepticism to these long distance determinations for Haffner 19, not only because of the morphological reasons presented by Snider et al. (2009), but also because the ionised gas around the cluster has essentially the same radial velocity as that of Sh2-311. It would be extremely surprising to see a distant very young cluster without any associated nebulosity through the molecular cloud associated with Sh2-311, especially if we consider the low reddening, E(B−V) ≈ 0.6 mag, to Haffner 19. Pending accurate Gaia distances, we advance that Sh2-311, at a Galactocentric distance of ~ 12 kpc, has a subsolar composition, in agreement with the Galactic abundance gradient (García-Rojas et al. 2005; Rodríguez & García-Rojas 2010) and thus the use of solar abundance tracks may result in overestimated distances.

HD 64315 lies in isolation near the centre of Sh2-311. Observations by Yadav et al. (2015) show that most of the bright stars in its surroundings are foreground late-B objects. Indeed, using photometry of a circular area of radius ~ 3′ around HD 64315 they only find two objects with colours compatible with a young population associated w the H ii region until they reach ~ 4 mag fainter. Given our total M_V ≈ −6, this means that there are at most two O-type or early-B stars (earlier than ~B2) in the immediate vicinity of HD 64315. This is a very unusual configuration for a star sitting in the middle of an H ii region with active star formation (Snider et al. 2009), since generally the mass of the most massive star in a cluster correlates with cluster mass (Weidner et al. 2010). Against a deterministic interpretation of this correlation, Oey et al. (2013) presented a sample of 14 OB stars in the SMC that meet strong criteria for having formed under extremely sparse star-forming conditions in the field. HD 64315 can be interpreted as a Galactic equivalent to these objects, However, it also represents a cautionary tale about the meaning of an “isolated” star.

Was this complex multiple system born in isolation in the middle of Sh2-311? Though unusual, this scenario looks quite likely. The only possible alternative, if Haffner 19 is really at the same distance as HD 64315, against recent analyses, is that HD 64315 is a runaway from this cluster. However, it is extremely difficult to conceive a dynamical interaction that can result in the ejection of such a wide binary without disrupting it. Of course, if the two binaries are not physically bound, the ejection scenario also breaks down, because they should have been ejected individually in exactly the same direction with exactly the same velocity, an even more improbable occurrence. Finally, the systemic radial velocity of binary A, which contains most of the system mass, is quite similar to that of the surrounding medium, again suggesting in situ formation.

We cannot exclude the possibility that HD 64315 is really a compact cluster; most of the mass is concentrated in the two observed binaries. Lower-mass components would not be observable in the glare of this very bright system. However, in any case, it is highly unlikely that HD 64315 is surrounded by a cluster with ~ 1000 M_⊙, as is usual for other O6 V stars (Weidner et al. 2010).

Binary separation in multiple stars is a possible indicator to discern between two of the main mechanisms proposed for stellar formation. Turbulent fragmentation leads to initial separations >500 AU (Offner et al. 2010), while disk fragmentation produces initial separations <500 AU. Unfortunately, the separation between the binaries A and B is around 500 AU, just between the predictions of the two theories. Another criterion used to distinguish between these two formation scenarios is the alignment of the stellar spin (Offner et al. 2016). When formed via disk fragmentation, stars have common angular momenta and therefore aligned stellar spin. The fact that the inclinations of the two binaries are quite similar supports in this case the disk fragmentation model, even though this argument has no statistical significance on its own.

9. Summary and conclusions

By using a complex procedure to analyse 52 high-resolution spectra, we are able to confirm that HD 64315 contains two binary systems, one of which is an eclipsing binary. The non-eclipsing binary (system A) has a period of 2.7 d, and is quite similar to the well-studied DH Cep. Its components are hotter and more luminous than those of the eclipsing binary (system B), and dominate the appearance of the system. System A is a detached binary composed of two stars with spectral types around O6 V, with minimum masses of 10.8 M_⊙ and 10.2 M_⊙, and likely masses ≈ 30 M_⊙. The eclipsing binary has a shorter period of 1.0 d, and produces a weak, but observable effect in the system light curve. We have derived masses of 14.6 ± 2.3 M_⊙ for both components, which are late O-type stars. They are almost identical: they overfill their respective Roche lobes and share a common envelope. System B is thus one of the most massive overcontact binaries known, and a very likely merger progenitor. Its merger within such a complex system may lead to the formation of a hierarchical triple system with three stars of similar masses.

We are not able to rule out an accompanying low-mass compact cluster with current observations, but HD 64315 has few nearby OB-type companions and does not appear to have been ejected from a nearby open cluster. It thus seems likely that the system, with a total mass above 90 M_⊙, formed in relative isolation near the centre of the Sh2-311 H ii region. In summary, HD 64315 is potentially a massive hierarchical system that formed in a sparse environment, which nicely highlights the need for detailed studies of multiplicity in apparently “isolated” stars.

Acknowledgments

This research is partially supported by the Spanish Government Ministerio de Economía y Competitivad (MINECO/FEDER) under grants AYA2015-68012-C1-1-P and AYA2015-68012-C2-2-P. This research has made use of the Simbad, Vizier and Aladin services developed at the Centre de Données Astronomiques de Strasbourg, France. Francesc Villardel acknowledges the support of the Spanish Ministry for Economy and Competitiveness (MINECO) and the Fondo Europeo de Desarrollo Regional (FEDER) through grant ESP2016-80435-C2-1-R, as well as the support of the Generalitat de Catalunya/CERCA programme.

References

All Tables

All Figures

	Fig. 1 Photometric data extracted from the All Sky Automated Survey (black dots) and Hipparcos catalogues (red dots). The V passband is shown in the top panel. The lower panel shows the B-band data.
In the text

	Fig. 2 Representative examples of the spectrum of HD 64315 at four different phases. Spectra are numbered for ease of identification.
In the text

Fig. 3 Interstellar lines in the spectrum of HD 64315. The main panel shows the components of the Na i D doublet (5889.95 Å solid line; 5895.92 Å dashed line) in velocity space. The inset shows the Galactic rotation curve along this line of sight. The velocities are measured with respect to the local standard of rest (LSR) assuming a solar motion of + 16.6 km s-1 towards Galactic coordinates ℓ = 53°; b = + 25°.

In the text

	Fig. 4 Results of the Lomb-Scargle periodogram for the photometric data in the V filter (left). The result of applying the clean algorithm is shown (right). The frequency peak corresponds to half the orbital period.
In the text

	Fig. 5 Photometric light curves (V filter, top; B filter, bottom), including data and the corresponding error bars from ASAS (blue open dots) and Hipparcos (black solid dots).
In the text

	Fig. 6 Comparison of several pairs of spectra taken on different nights, but corresponding to the same phase (blue numbers) of the 1.019 d photometric period. The presence of more than two components is obvious.
In the text

	Fig. 7 Final fit: the observed lines (blue), the sum of four Gaussian functions (black), four Gaussian functions in every spectral line (Aa: blue; Ab: red; Ba: green; Bb: magenta). Numbers represent the order of the spectrum (for every pair: left He i and right He ii).
In the text

	Fig. 8 Radial velocity curves for the non-eclipsing binary (A), fitted to the observational data and shown against orbital phase (blue line: star Aa; red line: star Ab). The red dotted line corresponds to the systemic velocity. The residuals are shown in the top panel.
In the text

	Fig. 9 Representative drawing of HD 64315 AaAb at different phases created with the phoebe 2.0-alpha code via the Python interface. The secondary is slightly greyer and smaller to differentiate it from the primary.
In the text

	Fig. 10 Radial velocity curves fitted to the observational data and shown against orbital phase (blue line: Ba; red line: Bb). The red dotted line corresponds to the systemic velocity. The residuals are shown in the top panel.
In the text

	Fig. 11 Light curve model fitted to the observational data and residuals in V filter (above) and B filter (below).
In the text

	Fig. 12 Representative drawing of HD 64315 BaBb to scale at quadrature phase, created with the phoebe 2.0-alpha code via the Python interface.
In the text

	Fig. 13 Radial velocity curves of the four stars of HD 64315 in time (above; Aa: blue, Ab: red, Ba: green, Bb: magenta). The He ii line of each spectrum is shown (solid black line) compared to the model from all components (dashed line).
In the text