© ESO 2020

1. Introduction

The formation of massive stars remains one of the most important open questions in astronomy today (e.g. Zinnecker & Yorke 2007; Tan et al. 2014). Observing the early phases of massive stars formation remains challenging at best: massive star formation is rare and found at large distances, in addition, their formation timescale is short and they are born in an environment strongly obscured by gas and dust.

Several formation scenarios have been proposed, among them: formation through stellar collisions and merging (Bonnell et al. 1998), competitive accretion (Bonnell et al. 2001; Bonnell & Bate 2006), and monolithic collapse (McKee & Tan 2003; Krumholz 2009). Except for the merger process, most theories are in agreement regarding the need for dense and massive accretion disks to overcome the radiation barrier. These disks likely fragment under gravitational instabilities (Kratter et al. 2010), which may result in the formation of companions, however model predictions are still scarce. Studying the correlations between multiplicity characteristics may provide crucial observational constraints for distinguishing between the different scenarios of massive star formation and assistance in the development of future theoretical models.

The multiplicity properties of massive stars have already been the subject of recent surveys in the Milky Way and nearby galaxies (for a recent overview, see e.g. Sana et al. 2017). Some studies focused on the spectroscopic analysis of young massive stellar clusters (Sana et al. 2012; Almeida et al. 2017) and OB associations (Kiminki & Kobulnicky 2012) and others on the high-angular astrometric observations of massive stars (Mason et al. 2009; Sana et al. 2014; Aldoretta et al. 2015; Gravity & Karl 2018) in order to determine the binary fraction of massive stars in these regions.

Among those previous studies, the Southern MAssive Stars at High angular resolution survey (SMaSH+, Sana et al. 2014) was an ESO Very Large Telescope (VLT) Large Program (P89-P91) that combined optical interferometry (VLTI/PIONIER) and aperture masking (NACO/SAM) to search for mostly bright companions (ΔH <  4) in the angular separation regime 0.001″ <  ρ <  0.2″, around a large sample of O-type stars. The entire NACO field of view was further analysed to search for fainter (ΔH <  8) companions up to 8″.

The SMaSH+ results show the importance of such studies for the understanding of massive star formation. They conclude that almost all massive stars in their sample have at least one companion and that over 60% have two or more. In addition, a larger number of faint companions are seen at large separations, corresponding roughly to the outer edge of the accretion disk. This is in agreement with expectations from the theory of disk fragmentation. These companions may correspond to outwardly migrating clumps resulting from the fragmented accretion disk or from tidal capture. Investigating whether low-mass companions exist at closer separations or if there is a characteristic length at which the flux versus separation distribution changes is therefore critical.

Nevertheless, large areas in the parameter space that have not been probed by these surveys remain, mostly due to instrumental limitations. In particular, while a rather complete view of companions down to mass ratio of about 0.3 has now been achieved, the existing surveys have so far failed to probe the lower-mass end of the companion mass-function. In the last few years, extreme adaptive optics (XAO) instruments have gone online, peering far deeper and more accurately than previously possible. The XAO, implemented at the VLT through the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2019), provides the necessary spatial resolution and dynamics to search for faint companions to nearby massive stars.

In this context, the Carina High-contrast Imaging Project of massive Stars (CHIPS) aims to characterise the immediate environment of a large sample of massive stars within 3° from η Car. A total of 93 O- & Wolf-Rayet type stars were selected from the Galactic O-Star Catalogue (GOSC, Maíz Apellániz et al. 2013) and the Galactic Wolf-Rayet catalogue (Rosslowe & Crowther 2015). So far about half of the potential targets have been observed, which ought to be sufficient to obtain constraints on the occurrence rate of companions in the SPHERE separation range with a precision better than 7%. SPHERE will allow us to investigate the presence and properties of massive star companions in the angular separation range of to (approx. 350–12 500 au) and Δmag ≈ 12 (mass-ratios > 0.03 on the main sequence). The range below a couple thousand au is particularly important as it corresponds to the approximated size of the accretion disk, where faint companions formed from the remnant of the fragmented disk could be found.

The present paper is the first in a short series. Here, we aim to establish a proof-of-concept using the first VLT/SPHERE observations of the QZ Car multiple system. QZ Car (≡ HD 93206) is a high-order multiple system composed of two spectroscopic binaries (Aa & Ac) and three previously resolved companions within 7″ (Ab, E & B). The pair (Aa1, Aa2) has a spectral type O9.7 I + B2 V, and an orbital period of 20.7 days. The pair (Ac1, Ac2) has a spectral type O8 III + O9 V, and a period of six days. These two binaries make up the central system (Aa, Ac), separated by roughly 30 milli-arcsec (mas) (Sana et al. 2014; Sanchez-Bermudez et al. 2017), and have a combined H and K_s-band magnitudes of 5.393 and 5.252, respectively. The companions Ab, E & B were detected by the SMaSH+ survey at separations of , and from the central system, respectively.

This paper is organised as follows. Section 2 presents the observations and data reduction. Section 3 describes the image post-processing algorithms as well as additional functionalities developed for our current studies. Our results are discussed in Sect. 4 and conclusions are presented in Sect. 5

2. Observations and data reduction

2.1. Observations

The QZ Car observations were obtained on Jan 25th, 2016 using the second-generation VLT instrument SPHERE, situated on the Unit Telescope 3 at the Paranal observatory in Chile. SPHERE is a high-contrast imaging instrument combining an extreme adaptive optics system, coronagraphic masks, and three different sub-systems with specific science goals. Our observations were executed in the IRDIFS extended mode (IRDIFS_EXT) mode using the Integral Field Spectrograph (IFS, Claudi et al. 2008) and the Infra-Red Dual-beam Imaging and Spectroscopy (IRDIS, Dolhen et al. 2008) sub-systems.

IFS images have a size of 290 × 290 pixels and a pixel size of 7.4 mas, hence corresponding to a field-of-view (FoV) of on the sky. The IRDIS camera has 1024 × 1024 pixels, covering a 12″ × 12″ FoV with a pixel size of 12.25 mas. The IRDIFS_EXT mode was chosen to allow combining the YJH-band observations with IFS to dual-band K-band observations with IRDIS. With its small FoV and spectroscopic capabilities, IFS allowed us to both detect and characterise companions at short separations. The larger FoV of IRDIS provided additional information on the local density of faint objects.

The observation sequence was composed of three types of observations: (i) CENTRE (C), allowing us to compute the centroid location of the coronagraph; (ii) FLUX (F), to obtain a reference flux point-spread function (PSF) of the central objects; and (iii) OBJECT (O), with the central star blocked by the coronagraph, hence delivering the scientific images that will be scrutinised to search for faint, nearby companions. FLUX observations were performed with the central star outside the coronagraph. The F-C-O sequence was repeated three times. Due to the brightness of QZ Car (H <  5.5), we used the neutral density filter ND2.0 – delivering a transmission of the order of 10⁻² – for the F and C observations of both instruments. The telescope was set to pupil tracking, that is, the centroid of the field is fixed on the science object and the sky rotates around it, as required for the later post-processing algorithm which uses the angular information of the movement of companions on an image, or angular differential imaging (Marois et al. 2006).

For the IFS OBJECT observations, we used 16 DITs of 4 s, for a total exposure time of ∼64 s. With IRDIS we chose 4 DITs of 4 s, giving an integration time of 16 s. Similarly we adopted NDIT × DIT of 4 × 16 s and NDIT × DIT of 4 × 8 s for the FLUX exposures with IFS and IRDIS, respectively. The observing setup and atmospheric conditions are detailed in Table 1. As a result of the OBJECT observations, we have obtained four-dimensional (4D) IFS and IRDIS data cubes of the QZ Car system. The IFS cubes are composed of spatial (2D) images of the IFS FoV for each of the 39 wavelengths channels (from 0.9 to 1.6 μm) and 48 sky rotations (due to the pupil tracking). The IRDIS data cubes contain 2D pixel images at each of the two wavelengths channels (K₁ and K₂) and 48 sky rotations, covering a total parallactic angle variation of about 3.5°.

2.2. Science data reduction

The data reduction of IRDIS and IFS images was processed by the SPHERE Data centre (Delorme et al. 2017, DC) at the Institut de Planetologie et d’Astrophysique de Grenoble (IPAG)¹. The SPHERE-DC process is standardised in terms of astronomical data reduction: removing bad pixels, dark and flat frames and estimating the bias in each exposure. They also calibrate the astrometry associated to the science frames using the on-sky calibrations from Maire et al. (2016), i.e. a True North correction value of 1.75 ± 0.08° and a plate scale of 7.46 ± 0.02 mas pixel⁻¹ for IFS and 12.255 ± 0.009 mas pixel⁻¹ for IRDIS. The system uses a modified version of the SPHERE Esorex pipeline² that is functional, automated and can be accessed by the user if requested. The end result from this data reduction is the reduced 4D science data cube, tables containing the wavelengths and rotational angles, and the 3D PSF cubes (see Sect. 2.3).

2.3. The point spread function of QZ Car

The FLUX observations are images of the central star taken without the coronagraph to estimate the PSF and flux of the central star. Three such PSF observations were taken during our observing sequence and delivered data cubes that contained the 2D images of the field at each of the instruments’ respective wavelength channels. To increase the signal-to-noise ratio (S/N) of the PSF frames and to average out the observing conditions, the median of the three PSF frames at each wavelength was computed so only a single 3D data cube was left. The obtained PSF is to be used for the companion modelling and characterisation techniques introduced in Sect. 3. However, in the case of QZ Car, additional complication arose in the IFS PSF frames.

The pair of close binaries at the core of QZ Car’s multiple system is separated by roughly 30 mas. This is just about the diffraction limit of an 8.2 m telescope in the Y to H band and sufficient for QZ Car’s point spread function (PSF) in IFS to display an elongated shape in the FLUX images. This had unintended consequences in the data processing as the reference PSF obtained from the IFS FLUX images is used to create a normalised PSF needed to inject artificial companions at different stages in the analysis, hence propagating the PSF’s deformation and leading to a number of artefacts. Therefore we adopted another reference PSF from an IFS observation of HD 93129A taken on the night of February 10, 2016. While HD 93129A is itself a long period binary, it was unresolved at the time of the observations (Maíz Apellániz et al. 2017). The flux of HD 93129A’s PSF is not the same as the original PSF of QZ Car. It was therefore scaled so that the new reference IFS PSF images have the same integrated flux as the QZ Car IFS PSF frames. This allows us to retain the original flux information while adopting a more representative PSF shape. Error estimation of the total flux measured from the PSF was accomplished by computing the standard deviation of the flux for all the FLUX images.

3. Data analysis

Once the data were reduced, image post-processing algorithms for high-contrast imaging were used on the target frames. In this study, we made use of the Vortex Imaging Processing package³ (VIP, Gomez Gonzalez et al. 2017) and a PSF-fitting technique (Bodensteiner et al. 2020). The VIP is an open-source python package for high-contrast imaging data processing that is instrument-agnostic. It was developed for exoplanet research, disk detection, and characterization. It is capable of performing angular, reference, and spectral differential imaging (ADI, RDI, SDI, respectively) and ADI-SDI simultaneously based on matrix approximation with the Principal Component Analysis (Amara & Quanz 2012; Soummer et al. 2012). We contributed to the implementation of 4D data analysis and IFS support. For comparison, a PSF fitting approach was also used to obtain the photometry of sources on IRDIS images and on companions which the VIP could not characterise because they were located too close to the edges of the detector (see Sect. 3.2.2).

Post-processed PCA/ADI images for the IRDIS and IFS science cubes are presented in Fig. 1. From the IRDIS image, previously known companions, Ab and E, are clearly distinguishable. However, these are not the only companions detected, as over a dozen of other point sources are now clearly visible on the image. The IFS image gives a close-up view of the system. The Ab companion is still present at the top of the image though cropped and cannot be analysed. To the south of the IFS image, a previously unknown companion (Ad) is clearly visible and delivering enough flux for a first spectral analysis (see Sect. 4.2).

Fig. 1. Post-processed science frames from IRDIS (left) and IFS (right) instruments. Angular and spectral differential imaging was applied with one principal component used for modelling and subtraction of the residual stellar light (see Sect. 3). The IRDIS image combines the two wavelengths channels K1 and K2. Dark circles on the images show the -diameter size of the coronagraph. The blue circles indicate the previously resolved companions Ab and E. White circles indicate the IRDIS sources detected at the 5σ level while yellow circles indicate the new companion QZ Car Ad. The IFS image is ; the IRDIS is 12″ × 12″.

3.1. Companion detection

A visual inspection of the final IRDIS and IFS PCA images displayed in Fig. 1 reveals a handful of rather bright companions and a larger number of much fainter point sources. To evaluate which ones are true detections, we first estimated their (S/N) using the S/N map function implemented in VIP. This module computes the S/N at every pixel of the frame as defined in Mawet et al. (2014). From this map, we set our detection limit to S/N = 5. As expected, companions Ab and E, as well as the new companion Ad (detected in both IFS and IRDIS), have a high S/N. There are 16 other companion candidates at separations beyond the IFS FoV and detected in the IRDIS image with S/N >  5. They are marked S1 to S16 in Fig. 1, yielding a total of 19 individual sources detected within from the QZ Car central quadruple system.

3.2. Source characterisation

Once we identified the true sources, we retrieved their position and contrast with respect to the central star. Starting from guess positions estimated from the post-processed frames, we measured accurate angular separations (ρ), position angles (PA), and flux contrast in each wavelength channel and for each companion with three different methods described below. Two methods included in the VIP package were used for Ad, Ab, and E companions and PSF fitting for all S sources. Final results are provided in Table 2 and Appendix A for IRDIS and IFS, respectively.

3.2.1. VIP

We first measured the flux of the stars using the aperture photometry spectral channel by spectral channel, which provided us with a first estimate of the spectrum. This initial guess was then passed onto the Simplex Nelder-Mead optimisation (hereafter referred to as the Simplex method) of VIP which estimated position and flux parameters by applying a NEGative Fake Companion technique (NEGFC). This method consists in inserting negative artificial sources in the individual frames, varying at the same time their brightness and position (starting from the values measured by aperture photometry). The artificial companions are obtained from the unsaturated PSF of the central star, measured in the FLUX observations. The residuals in the final images are then computed and compared to the background noise, measured in an annulus at the same radial distance. The combination of brightness and location that minimises the residuals are estimated through a Nelder-Mead minimization algorithm. This provides reference fluxes for each spectral channel. Dividing the fluxes of the different companions by the reference fluxes from the PSF data cubes, we obtained flux or magnitude contrasts at each wavelength channel.

Results from this simplex method were then injected into the MCMC engine (also available in VIP) to estimate confidence intervals for the angular separations, PAs and fluxes for each wavelength channel and each target. In this way, we obtained an (uncalibrated) IFS low-resolution spectrum for companion Ad and two flux values for all IRDIS sources. For sources with small radial distance difference between each other (Δρ <  200 mas), masks were applied on sources other than the target. This was to prevent issues with parameter estimation when multiple sources are in the same annulus when applying the simplex optimisation routine. For this purpose, circular masks were taken at the same radial distance as the source we wished to mask, but at a different position angle to preserve the noise properties. A rotation was then applied on the mask in order to preserve the radial dependence of the noise. For Ab, the MCMC algorithm failed to retrieve the K₁ uncertainties. The MCMC failed to converge for S13 to S16 while nor Simplex nor MCMC could be used for S15 and S16. In either case, this is caused by the fact that the field-of-view is too small to contain the full annulus at the target separation, meaning that the sources are too close to the edge of the field.

In order to test the accuracy of the MCMC confidence intervals for IRDIS sources, we further implemented a Monte Carlo (MC) method. Firstly, all detected sources are masked. Secondly, a number of artificial sources (25 in our case) are injected at the same radial distance and with the same flux as a given companion, but at varying position angles. The flux and position of the artificial companions are then measured using the Simplex algorithm and compared to the input values. The corresponding standard deviations finally yield an estimate of the 1σ error on the flux and positions of the considered companion. The process in then repeated for all sources. The results will be discussed in the next section.

The previous methods only take into account the statistical uncertainties from the image processing. For the final photometric errors we took into account the flux variations of the unsaturated images of the central star described in Sect. 2.3. For the calculation of the astrometric uncertainties, we adopted the plate scale and astrometric calibration precision given by Maire et al. (2016) and the ESO SPHERE user manual. The final astrometric errors are obtained by a quadratic sum of the Simplex-MC measurement errors, the star centre position uncertainty (1.2 mas, from Zurlo et al. 2016), the plate scale precision of 0.021 mas pix⁻¹ for IRDIS and 0.02 mas pix⁻¹ for IFS, the true north uncertainty (±0.08deg), and the dithering procedure accuracy (0.74 mas, Zurlo et al. 2016).

3.2.2. PSF fitting

For all IRDIS candidates beyond 2″, the central star’s influence is limited and the background noise dominates (see Sect. 3.4). Therefore, the use of ADI and SDI techniques is not necessarily needed to derive precise astrometry and photometry. A more widely used strategy in astronomical imaging is to use a PSF-fitting technique which provides accurate position and flux values for the S sources (Ad, Ab & E ). Our PSF-fitting method is based on the photutils⁴ python package along with an effective PSF model developed by Anderson & King (2000) and is described in Bodensteiner et al. (2020).

For this, we used the derotated and collapsed images which maximise the S/N in both K₁ and K₂. The FLUX observations with IRDIS are used to establish an accurate PSF model. This is then fitted to each source individually in order to obtain accurate positions and flux estimates. This technique is very useful for sources that are detected in the collapsed images but are too close to the edges of the frames for ADI techniques, that is, S15 and S16. However, in the K₂ band, PSF fitting could not converge to a good solution for sources S10, S13, and S14 as they are too close to the detection limit and do strongly benefit from the ADI post-processing.

A comparison plot of the magnitude contrasts obtained between the three methods used in this paper is shown in Fig. 2 for the K bands. A comparison between the X, Y coordinates can be found in Appendix B. These figures show that the positions and (in most cases) contrast values from the three methods are in excellent agreement, for sources S1 to S8, beyond which discrepancies arise. Sources too close to the edge could not be fitted with Simplex (S15-S16) nor MCMC (S13-S16), while the PSF-fitting approach failed for the faintest object in the K₂-band (S10, S13, S14). Aside from these differences in the best-fit values, MCMC usually yields magnitude contrast errors that are a factor of few larger than those obtained with the Simplex+MC approach and with PSF-fitting. It is likely that the true errors lay somewhat in the middle. From the model atmosphere-fit in the next section, we have strong evidence that the MCMC contrast errors are likely overestimated by a factor of three to four while the Simplex+MC and PSF-fitting uncertainties are too small by a factor of about two. The order of magnitudes is however correct. For future work and given the fact that the MCMC is much more computationally intensive, the present comparison certainly favours the use of PSF-fitting for most sources, but the one that are the closest to the detection limit. Final results adopted from the Simplex+MC or from the PSF fitting techniques are given in Table 2.

Fig. 2. Left: K1. Right: K2. IRDIS contrast magnitudes and associated errors of S sources obtained with three different methods: MCMC, Simplex-MC and PSF-fitting. Some errors for Simplex-MC and PSF fitting are small and not seen on these plots.

3.3. Flux calibration

To obtain the absolute fluxes of the companions, we would require a flux calibrated spectrum of the central QZ Car system in the same wavelength range as that of our SPHERE observations (Y to K); unfortunately, no such spectrum is available. To circumvent this issue, we modelled the spectral energy distribution of QZ Car’s central quadruple system using the non-local thermodynamic equilibrium (non-LTE) atmosphere code FASTWIND (Puls et al. 2005; Rivero-González et al. 2011). Each component of the QZ Car’s central system was modelled separately and their contribution within the PSF was then combined. The parameters for the computation were obtained using the spectral types from Sanchez-Bermudez et al. (2017) and the corresponding O-star calibration tables from Martins et al. (2005). Parameters for the Aa2 component were found using a combination from Trundle et al. (2007) and Parkin et al. (2011). We also calculated the mass-loss rate (Ṁ) and terminal wind velocities (v_∞) for each stellar component following Vink et al. (2001) as these are needed input for FASTWIND. Our results are summarised in Table 3.

Once the spectra for all the four central components were calculated and combined, we multiplied the contrast fluxes calculated previously by the model spectrum of QZ Car to obtain the absolute fluxes of each companion in the different wavelength bands. In particular, the IFS+IRDIS flux-calibrated spectrum of the Ad companion is displayed in Fig. 3. Throughout this process and later on in Sect. 4, it’s necessary to adopt a reference radius for the sphere at the surface of which the flux is computed. Without loss of generality, we arbitrarily adopted a value of 100 R_⊙. We emphasise that this value has no physical meaning.

Fig. 3. IFS+IRDIS flux-calibrated spectrum of QZ Car Ad at a reference distance of 100 R⊙. The plain line gives the best-fit ATLAS9 model with Teff = 8896 K, log g = 4.27 and R = 1.72 R⊙. Shaded area represents the 1-σ uncertainties on the observed spectra.

3.4. Detection limits

In this section, we estimate the sensitivity of our observations in terms of magnitude difference as a function of the angular separation ρ to the central object. Using the VIP contrast curve modules we computed the contrast limits for a chosen σ level by injecting artificial stars (based on the scaled PSF of HD 93129A, see Sect. 2.3) and calculating the noise at different radial distances from the centre. This implementation takes into account the small sample statistics correction proposed in Mawet et al. (2014). In order to avoid interference from the bright companions, all sources were masked (see Sect. 3.2.1). Although this significantly increases the quality of the contrast curves, small artefacts with a 0.2 mag amplitude remain visible in the contrast curves at a radial separation of 3″. The 5-σ sensitivity curves that we obtained are presented in Fig. 4.

Fig. 4. IRDIS (left) and IFS (right) detection limits (blue). The contrasts are given in the H-band for IFS and a mass scale is also provided on the right-hand axis. These masses were estimated using pre-main sequence and main-sequence tracks from Siess et al. (2000). Bumps at in the IFS contrast curves results from the trace of the deformable mirror of the instrument at this radial distance.

A contrast better than 8 mag is achieved at 200 mas with IFS, and as large as 11 mag at ρ >  600 mas. These magnitude differences correspond to flux contrasts of 1.5 × 10⁻⁴ and 4 × 10⁻⁵, respectively. Such detection limits are in line with past SPHERE observations (Zurlo et al. 2016; Mesa et al. 2019) in IRDIFS_EXT mode, if we consider that our total exposure time was roughly a minute with both IFS and IRDIS. Using a mass scale from Siess et al. (2000)’s evolutionary tracks, stars with masses < 1 M_⊙ could be easily detected by this system. Similarly, IRDIS delivers contrast better than 9 mag at and of better than 13 mag at separations larger than . This is about 5 mag deeper than previously achieved with SMaSH+, and up to 8 mag deeper in the poorly mapped region around 400 mas demonstrating the complementing capabilities of SPHERE with respect to previously obtained high-angular resolution observations of massive stars.

3.5. Distance of QZ Car

Knowledge of the distance of QZ Car with respect to the observer is crucial to convert the angular separation in physical (projected) separations. In an attempt to improve on the distance of QZ Car, we retrieved its astrometric information from the Gaia DR2 catalogue (Gaia Collaboration 2016, 2018; Lindegren et al. 2018), including positions, parallax, proper motion, their uncertainties and correlations. We used a galactic prior with a length scale of 2.5 kpc (Walborn 2012) and performed an MCMC fit to the distance and proper-motion vector as described by Bailer-Jones (2017). We obtained a distance of kpc. Similar results are obtained with a flat prior. This is in contrast with prior estimates of about 2.3 kpc (Walborn 1995; Smith 2006). We also computed a spectral distance modulus yielding a distance of about 2.0–2.1 kpc. Given these discrepancies, we scrutinised further the Gaia measurements. We computed the quality of the astrometric fit the so-called RUWE indicator of Lindegren et al. (2018), leading to RUWE = 2.54, which reveals a poor astrometric fit despite the reasonable relative uncertainty of 11% on the parallax, and the 299 good astrometric measurements. This poor RUWE is likely caused by the multiplicity of the central object, which biases the Gaia parallax. This issue could be resolved in the next Gaia data release, where binarity will be taken into account when deriving the parallax.

In the remainder of the paper, we adopt 2.3 kpc as the reference distance to convert the measured angular separations to (projected) physical distances. The results are given in Table 2. We note that the only results in this work that depend on the adopted distance are the projected physical separations, so that there is no impact on the presented results if the adopted distance to QZ Car would turn out to be incorrect.

4. Results and discussion

4.1. Probabilities of spurious association

While over a dozen faint sources are clearly resolved around QZ Car in Fig. 1, additional information, such as common proper motion, would be needed to confirm their physical connection to QZ Car. Unfortunately, we only have one observation of QZ Car so far and the closest object to QZ Car that was also detected by Gaia is at 7.2″, i.e. outside the IRDIS field-of-view. In the absence of such information, we resorted to a statistical argument. To this aim, we define the probability of spurious association (P_spur(ρ_i|Σ(K_i))) as the probability that at least one source is found by chance at a separation ρ equal or closer to QZ Car than that of the companion i (ρ ≤ ρ_i) given the local source density Σ of stars at least as bright as i (K ≤ K_i)⁵.

A query of the VISTA Carina Nebula catalogue (Preibisch et al. 2014) yielded N_obj = 1864 within a r = 2′ radius around QZ Car down to magnitudes of K_s ≈ 19. To compute P_spur, we first estimate the local source density Σ(K_i) = N_obj(K ≤ K_i)/(πr²) of objects at least as bright as the companion i. We then use a Monte Carlo approach and randomly generate 10 000 populations of N_obj(K ≤ K_i) stars uniformly distributed in πr². The probability of spurious association is finally obtained as the fraction of populations in which at least one star is to be found at ρ ≤ ρ_i.

This simple exercise confirms the very low probability of spurious association (P_spur <  0.02) – hence, the high confidence regarding physical association – of companions Ab, Ad, and E while the presence of most of the fainter ‘S’ sources is best explained by chance alignment given the overall surface density of sources in QZ Car’s surroundings.

In addition to computing the probability of spurious associations, we ran the Besançon model of the Galaxy (Robin et al. 2003) in the direction of QZ Car. The model predicts about 20 stars with K-band magnitude brighter than 19 mag in an area corresponding to our IRDIS field of view. All of them have K >  15.5 and the vast majority (18/21) are background stars (distance > 3 kpc). To the first order, this is compatible with the properties of the S sources in the IRDIS field of view and provides an additional argument to consider (most of) them as chance alignments. It also supports the fact that Ad, Ab, and E are physically connected as the Besançon model cannot explain the presence of such bright objects around QZ Car. Sources S1 to S6 have 0.25 <  P_spur <  0.50, so that, statistically, some of these could still be physically linked. In summary, and accounting for the four inner objects, QZ Car consists of seven likely physical companions within an ≈25-radius, as well as additional three to four fainter candidate companions within a similar angular separation.

Confirmation of common proper motion and characterization of orbital motion are of course crucial to definitely prove any physical association. For the closest companion, Ad, given the precision of our astrometry (see Table 2) and the proper motion of the central star (Gaia Collaboration 2018), it should be possible to prove common proper motion and measure a significant orbital rotation (excluding contamination by high proper motion background objects, see Nielsen et al. 2017) with observations separated by 1 and 7 yr, respectively (assuming a circular orbit).

4.2. Spectral modelling of Ad

The flux calibrated low-resolution IFS+IRDIS spectrum of Ad provides us with information about its spectral energy distribution (SED). Here we attempt to use that information to constrain the stellar parameters of the Ad companion for the first time.

The uncalibrated IFS spectrum is mostly flat, except for two broad absorption features at 1120 nm and 1372 nm (also seen on the calibrated spectrum on Fig. 3). These correspond to the expected location of Earth telluric bands. To check this, we adjusted a synthetic telluric spectrum to the Ad data using MolecFit (Smette et al. 2015; Kausch et al. 2015). We used the atmospheric conditions at the time of the observations and only consider telluric lines resulting from water. Once corrected for the telluric bands, and aside from (unphysical) edge effects, the Ad spectrum is feature-less as can be expected from its very-low spectral resolution. The present exercise is useful for confirming the origin of the broad absorption in the IFS Ad spectrum and to verify spectral ranges that can in principle be well corrected from telluric absorption. The MolecFit telluric corrected spectrum is, however, not used in the following as the division by the spectrum of the central object actually take care of the telluric correction automatically.

To constrain the stellar parameters from the spectrum, we used ATLAS9 LTE atmosphere models (Castelli et al. 2004). We associated an ATLAS9 model to each time step in the pre-main sequence (PMS) evolutionary tracks of Siess et al. (2000) and quantitatively compared it to Ad’s SED. For numerical reasons, in the computation of the χ²-contour curves, we also interpolated the ATLAS9 grid along the log g-axis.

Each ATLAS9 model was converted into flux at the reference distance of 100 R_⊙ and integrated over the width of each IFS and IRDIS wavelength channels, allowing us to compute the corresponding χ² accounting for the error-bars on Ad’s SED. Adopting the MCMC errors, the best-fit χ² has a reduced value of 0.1, suggesting that the error bars are heavily overestimated. Adopting the Simplex+MC values yields a best-fit reduced- χ² of about 4.4 suggesting that the latter are underestimated by a factor of about two. In the following, we rescale the Simplex+MC errors so that the best-fit reduced- χ² is equal to unity. The obtained χ²-map is displayed in the Hertzsprung-Russell diagram (HRD) of Fig. 5.

Fig. 5. χ2-surface resulting from ATLAS model fit to Ad’s IFS+IRDIS SED projected onto the HRD plane. The thick black and red lines are ZAMS and early-MS according to Siess et al. (2000) definitions (see footnote 5). Thin black lines are the evolutionary PMS tracks for stars with masses increasing from bottom to top from 1.1 to 3.5 M⊙. The 4 to 8 Myr parts of the respective evolutionary tracks are displayed with a thicker line. Finally, the red dashed lines give, from top to bottom, the 1, 5 and 10 Myr isochrones. The best-fit model is indicated with a white diamond.

As expected, our fit results in multiple possible combinations but the allowed mass remains limited in the range of 1.8 to 2.2 M_⊙. The best fit model is obtained for a 2.0 M_⊙ star with T_eff = 8896 K, log g = 4.27, L = 17.7 L_⊙ and R = 1.72 R_⊙. The best fit spectrum from ATLAS9 shows a nice fit with the measured spectrum of Ad in Fig. 3. However, the precision of the model retrieval is limited by two factors: the anticipated degeneracy between physical surface properties and evolutionary stage and the density of the model grid. The first is clearly illustrated in Fig. 5 by the elongated χ²-valley in the HRD. Focusing on the location of the best-fit model, there are statistically significant differences in the goodness-of-fit between the 2.0 M_⊙ PMS-tracks and the neighbouring 1.9 and 2.2 M_⊙-tracks. All other things being kept equal, these translates into statistical uncertainties of the order of 300 K and 0.2 dex in T_eff and log L/L_⊙, respectively. Had we adopted the MCMC errors, T_eff and L/L_⊙ values of 7500 to 9600 K and 11 to 23 L_⊙ would have been obtained within a 68% confidence interval, respectively.

According to our best-fit solution, the star’s evolutionary stage seems to be somewhere between Siess et al. (2000)’s ZAMS and early-main sequence⁶. Its age of 9.7 Myr is probably in fair agreement with QZ Car age estimates (Walker et al. 2017). The best-fit ATLAS9 model is displayed in Fig. 5 and should correspond to stars of spectral type A3, according to Siess et al. (2000) calibrations.

4.3. Physical properties of the IRDIS companions

IRDIS observations only provide us with two independent wavelength channels, K₁ and K₂ with central wavelengths of 2.110 and 2.251 μm, respectively. While this is insufficient to constrain the shape of the SED, it provides an important anchor point to assess the objects absolute K-band magnitude, assuming that the objects are located at the same distance and suffer from the same K-band reddening as QZ Car.

As for the data of QZ Car Ad, we computed χ² maps by comparing the K₁ and K₂ absolute fluxes of each companion sources to ATLAS9 LTE models. Good fits were obtained for most sources. The Ab companion seems to be more massive and we used the Brott et al. (2011) main-sequence evolutionary tracks rather than the PMS-tracks from Siess et al. (2000). The resulting χ²-maps, over-plotted in the HRD together with Siess et al. (2000)’s PMS evolutionary tracks and isochrones are displayed in Fig. 5 and 6. In this exercise, there exists, of course, a degeneracy between the age and the mass.

Fig. 6. Same as Fig. 5 for the IRDIS sources Ab and E.

While more wavelength channels would be desirable, our results show that the companions Ab, Ad, and E are compatible with the hypothesis of co-eval formation together with the central OB quadruple system in QZ Car. Adopting the QZ Car age range, that is, 4 to 8 Myr, first-order constraints on the masses of the individual companions can be obtained. With only 4.4 mag contrast with the QZ Car central system, the companion Ab is the most massive object and has a mass estimate of 10 to 12 M_⊙. The solution for E is fully degenerate with multiple minima on the χ²-map. We emphasise two of the extremer solutions: the best-fit is a low-mass (0.5 M_⊙) very young (2.5 Myr), cool (T_eff = 3770 K), and rather large (R = 3.7 R_⊙) PMS star. This solution is at the limit of the grid towards the low T_eff so that it is possible that the cooler model may provide an even better fit. The other minimum is a 2.5 M_⊙, 5.5 Myr-old star with T_eff = 10 kK, which is also a valid solution within 3-sigma. The latter is compatible with co-evality.

We also performed this exercise for the faintest S sources, implicitly assuming that they are located in the Carina region. Under this hypothesis, sources S1 to S16 are compatible with being low-mass pre-main sequence stars with masses in the range of 0.3 to 0.5 M_⊙. But for S6 that could be as young as 10 to 15 Myr if its mass is on the low-side of the confidence interval, all the other sources seem to be older than 20 Myr. Alternatively, they are even younger, lower-mass PMS star falling outside of the ATLAS9 model-grid or they are foreground or background sources unconnected to the Carina region, as discussed in Sect. 4.1.

4.4. Comparison with previous high-angular resolution surveys

Our results can be directly compared to previous high-angular resolution campaigns, such as those provided by the SMaSH+ and HST-FGS surveys (Fig. 7, Sana et al. 2014; Aldoretta et al. 2015). Clearly, by enabling the detection of sources at least 5 mag fainter than previously possible, SPHERE opens a new discovery space to investigate the low-mass end of the companion mass function of massive stars.

Fig. 7. QZ Car’s sources detected with SPHERE (star symbol) overlaid on the SMaSH+ (Sana et al. 2014) and HST-FGS (Aldoretta et al. 2015) companion detections in the magnitude contrast vs. angular separation plane. The thick lines give the limiting contrast curves of the different instruments (see legend).

In this context, we note a clear clustering of the ‘S’ sources in Fig. 7, that are all located at angular separations larger than 2″ (i.e. projected physical separation > 4.6 × 10³ au). Furthermore, there is a clear gap of 4- to 5 mag between the more massive, very likely physical and probably co-eval companions Ab, Ad, and E and the rest of the ‘S’ sources, most of which are either older or unconnected, as suggested by the large spurious alignment probabilities that we derived.

4.5. Future prospects

The spectrum of QZ Car Ad that we obtained with IFS in the IRDIFS_EXT mode has a spectral resolving power (λ/Δλ) of ∼50 only. At such low resolution, all spectral features but the telluric bands are smeared out (see Fig. 3). To better characterise the companion physical properties and age, a higher-resolution spectrum is desirable. However, seeing-limited spectrographs such as VLT/X-shooter will not be able to resolve the Ad companion. Very few AO-assisted spectrographs exist and, among those, almost none can deliver the require flux contrast. In Fig. 8, we investigate the resolving power of the Long Slit Spectroscopic (LSS) mode of IRDIS, delivering a spectral resolving power of 350, which would provide a valuable improvement to the current IFS SED and help us to estimate the parameters of the companion with greater accuracy.

Fig. 8. Spectrum of QZ Car Ad (blue). Same figure as Fig. 3 with the ATLAS9 (Castelli et al. 2004) models downgraded to the resolution of IRDIS-LSS (orange; λ/Δλ = 350).

5. Conclusions

In this work, we present the first SPHERE observations of QZ Car, a known quadruple system in the Carina region. Using the IRDIFS_EXT mode, we detected 19 sources in a 12″ × 12″ field-of-view; all but two (Ab and E) are newly detected. We used the high-contrast imaging software VIP, which is used for planet detection, to characterise the detected sources. Three of our sources (Ab, Ad, and E) are moderately bright, with K-band magnitude contrasts in the range of ∼4–7.5. The remaining sources have a magnitude contrast from 10 to 13. Most of the latter can be explained by spurious alignment given the source number density around QZ Car.

We further determined the limiting contrast curves, showing that SPHERE detection capabilities can reach contrasts better than 9 mag at only 200 mas and better than 13 mag at angular separations larger 2″. Our observations are sensitive to sub-solar mass companions over most of the angular separation range provided by SPHERE.

Finally, we used the known distance to QZ Car, a grid of ATLAS9 models and pre-main sequence evolutionary tracks to obtain a first estimate of the physical properties of the detected objects. This determination implicitly assumes that the companions are located at the same distance and suffer from the same NIR reddening than QZ Car central system. We found masses across the entire mass range, from a fraction of a solar mass up to 12 M_⊙, including a ∼2.0 M_⊙ companion at a (projected) physical separation less than 1700 au. While there is a degeneracy in the physical parameter versus age determination given the limited constraints, the three most massive, likely physical companions (Ab, Ad, and E) can be fitted with ages of 4 to 9 Myr, meaning that their formation is potentially contemporaneous to that of the inner quadruple system making QZ Car one of the highest order multiple systems known.

Future work can follow two directions. On the one hand, a better characterisation of the detected companions is desirable and will ultimately provide an independent age diagnostic. This will help to confirm physical connection of the companion through proper motions as well as high-resolution spectroscopy and may be possible with the SPHERE LSS for the brightest companions of QZ Car.

On the other hand, and based on the present results, it is clear that SPHERE is opening a new parameter space to investigate the presence and physical properties of faint companions within only a few 1000 au from massive stars. Performing similar observations of the entire sample of massive stars may allow us to investigate the outcome of the massive star formation process as well as to investigate the pairing mechanism of these faint companions.

Acknowledgments

This work is based on observations collected at the European Southern Observatory under programmes ID 096.C-0510(A). We thank the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris) and supported by a grant from Labex OSUG@2020 (Investissements d’avenir a ANR10 LABX56). We especially thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their help during the data reduction process. We acknowledge support from the FWO-Odysseus program under project G0F8H6N. This project has further received funding from the European Research Council under European Union’s Horizon 2020 research programme (grant agreement No 772225) and under the European Union’s Seventh Framework Program (ERC Grant Agreement n. 337569). LAA acknowledges financial support by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior and Fundação de Amparo à Pesquisa do Estado de São Paulo. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. JDR acknowledges the BELgian federal Science Policy Office (BELSPO) through PRODEX grants Gaia and PLATO. VC acknowledges funding from the Australian Research Council via DP180104235. Facilities: VLT UT3 (SPHERE).

References

Appendix A: Spectrum of QZ Car Ad

Appendix B: Comparison between parameter estimation techniques: Simplex, MCMC and PSF-fitting

Fig. B.1. Comparison of the X coordinates and associated errors. Most errors for Simplex and PSF fitting are small and not seen on these plots. Left: K1. Right: K2.

Fig. B.2. Comparison of the Y coordinates and associated errors. Most errors for Simplex and PSF fitting are small and not seen on these plots. Left: K1. Right: K2.

All Tables

All Figures

Fig. 1. Post-processed science frames from IRDIS (left) and IFS (right) instruments. Angular and spectral differential imaging was applied with one principal component used for modelling and subtraction of the residual stellar light (see Sect. 3). The IRDIS image combines the two wavelengths channels K1 and K2. Dark circles on the images show the -diameter size of the coronagraph. The blue circles indicate the previously resolved companions Ab and E. White circles indicate the IRDIS sources detected at the 5σ level while yellow circles indicate the new companion QZ Car Ad. The IFS image is ; the IRDIS is 12″ × 12″.

In the text

	Fig. 2. Left: K1. Right: K2. IRDIS contrast magnitudes and associated errors of S sources obtained with three different methods: MCMC, Simplex-MC and PSF-fitting. Some errors for Simplex-MC and PSF fitting are small and not seen on these plots.
In the text

	Fig. 3. IFS+IRDIS flux-calibrated spectrum of QZ Car Ad at a reference distance of 100 R⊙. The plain line gives the best-fit ATLAS9 model with Teff = 8896 K, log g = 4.27 and R = 1.72 R⊙. Shaded area represents the 1-σ uncertainties on the observed spectra.
In the text

	Fig. 4. IRDIS (left) and IFS (right) detection limits (blue). The contrasts are given in the H-band for IFS and a mass scale is also provided on the right-hand axis. These masses were estimated using pre-main sequence and main-sequence tracks from Siess et al. (2000). Bumps at in the IFS contrast curves results from the trace of the deformable mirror of the instrument at this radial distance.
In the text

Fig. 5. χ2-surface resulting from ATLAS model fit to Ad’s IFS+IRDIS SED projected onto the HRD plane. The thick black and red lines are ZAMS and early-MS according to Siess et al. (2000) definitions (see footnote 5). Thin black lines are the evolutionary PMS tracks for stars with masses increasing from bottom to top from 1.1 to 3.5 M⊙. The 4 to 8 Myr parts of the respective evolutionary tracks are displayed with a thicker line. Finally, the red dashed lines give, from top to bottom, the 1, 5 and 10 Myr isochrones. The best-fit model is indicated with a white diamond.

In the text

	Fig. 6. Same as Fig. 5 for the IRDIS sources Ab and E.
In the text

	Fig. 7. QZ Car’s sources detected with SPHERE (star symbol) overlaid on the SMaSH+ (Sana et al. 2014) and HST-FGS (Aldoretta et al. 2015) companion detections in the magnitude contrast vs. angular separation plane. The thick lines give the limiting contrast curves of the different instruments (see legend).
In the text

	Fig. 8. Spectrum of QZ Car Ad (blue). Same figure as Fig. 3 with the ATLAS9 (Castelli et al. 2004) models downgraded to the resolution of IRDIS-LSS (orange; λ/Δλ = 350).
In the text

	Fig. B.1. Comparison of the X coordinates and associated errors. Most errors for Simplex and PSF fitting are small and not seen on these plots. Left: K1. Right: K2.
In the text

	Fig. B.2. Comparison of the Y coordinates and associated errors. Most errors for Simplex and PSF fitting are small and not seen on these plots. Left: K1. Right: K2.
In the text