© A. Musso Barcucci et al. 2019

1. Introduction

Circumstellar discs are the natural by-products of the protostellar accretion process and they are the birthplaces of planetary systems. They evolve with time and undergo a series of processes. After the material that was in the original protoplanetary disc has been dissipated (usually within ∼10 Myr, see Ercolano & Pascucci 2017), a new generation of dust is created and continuously replenished via planetesimal collisions, forming a second generation debris disc (DD). These destructive encounters are triggered when the planetesimals are dynamically excited such that their relative velocities increase above a critical value (low-velocity collisions can happen in non-excited DDs as well, but they produce a different and recognisable emission spectrum, see Heng & Tremaine 2010). Three possible stirring processes have been proposed so far that could induce such an excitation in the disc: stellar encounters, self-stirring and companion-stirring. Of these three, the first scenario is the least likely one to be observed, since close stellar encounters are rare (particularly among field stars) and the disc brightness resulting from dust production drops too quickly to be detectable (Kenyon & Bromley 2002).

In the self-stirring scenario (Kenyon & Bromley 2008), planetesimals with low relative velocities form increasingly large bodies that in return dynamically excite smaller neighbours above the critical threshold for planetesimal destruction. The planetesimal growth scales with orbital period, resulting in an inside-out collisional cascade. Since a maximum growth speed is set by the host star and disc parameters, at any given time there is a maximum disc size that can be explained by self-stirring.

In the companion-stirring case (Mustill & Wyatt 2009), the planetesimals are excited by the companion’s secular perturbations, and the maximum disc size at a given time is a function of the physical properties of both the central star and the companion.

The optimal scenario to investigate these processes is therefore the one in which the disc and the companion(s) are observed and characterised at the same time. At the moment of writing, only a handful of such systems have been found: HR 8799 is one of the most extensively studied (Marois et al. 2008), alongside HD 95086 (Rameau et al. 2013) and β Pic (Lagrange et al. 2010). In addition, only two systems are currently known where the companion is in the stellar mass regime: HR 2562 (Konopacky et al. 2016) and HD 206893 (Milli et al. 2017).

The limited number of systems suitable to investigate the companion-disc interaction does not allow us to fully comprehend this phenomenon, and therefore augmenting this sample is our primary goal.

Detecting and characterising giant planets (GPs) around DD host stars is one of the scientific goals of the Imaging Survey for Planets around Young Stars (ISPY, Launhardt et al., in prep.), currently underway at the Very Large Telescope (VLT). It makes use of the NaCo instrument (Lenzen et al. 2003; Rousset et al. 2003) in L′ band, and the Angular Differential Imaging (ADI, Marois et al. 2006).

2. HD 193571

Within the NaCo-ISPY survey, we observed HD 193571 (HR 7779, GJ 969, κ 01 Sgr), an A0V field star at a distance of 68.45 pc (Gaia Collaboration 2018), which is part of a wide-separation (> 40″) three-component system¹ (WDS Catalogue, see Mason et al. 2014).

The age of this target is uncertain: David & Hillenbrand (2015) derived stellar parameters for more than 3000 nearby early-type (BAF) field stars, and compared them with stellar isochrones. They computed final ages and masses with both a Bayesian inference approach and classical isochrone interpolation, obtaining 161 Myr and 66 Myr, respectively. They presented criteria to decide between the two values, but for HD 193571 it is unclear which age or mass estimate should be preferred. Throughout this study we use a primary mass of M = 2.2 ± 0.1 M_⊙, which encompasses both the Bayesian inferred mass and the mass derived via interpolation.

The age estimates for HD 193571 are presented in Table 1, together with the main stellar properties.

HD 193571 is known to harbour a debris disc, inferred from its infrared excess (f = L_disc/L_⋆ = 2.3 × 10⁻⁵). We fit its spectral energy distribution (SED) to derive the stellar luminosity and effective temperature, and the debris belt radius. We fit simultaneously a stellar atmosphere (PHOENIX; Husser et al. 2013) plus a single black-body (BB) model to the observed photometry and the Spitzer IRS spectrum. The photometry includes a wide range of filters and wavelengths, from: “Heritage” Stromgren and UBV (Paunzen 2015), 2MASS (Skrutskie et al. 2006), HIPPARCOS/Tycho-2 (ESA 1997), AKARI (Ishihara et al. 2010), WISE (Wright et al. 2010), and Spitzer (Chen et al. 2014). The fitting method uses synthetic photometry of grids of models, and finds the best-fitting model with the MultiNest code (Feroz et al. 2009). The SED of HD 193571 is best fit by an A0 stellar model plus a one-temperature BB model locating the dust at a distance of R_BB = 62 ± 4 au, with a temperature of 81 ± 3 K. The best fit is shown in Fig. 2.

Fig. 1. Classically ADI reduced images for the two NaCo datasets (left and centre) and for the GPI dataset (right). The images are oriented with north up and east left, and the green cross indicates the position of the central star. The companion is clearly visible close to the centre in all three datasets. The images are normalised and the colour map was chosen for a better visualisation of the data.

Fig. 2. Flux density distribution of HD 193571, showing the photometric datapoints found in the literature (in blue) and the IRS spectrum (in black), together with the fitted stellar (green) and disc (red) fluxes.

The BB radius of the dust disc is given by (Pawellek & Krivov 2015)

$$\begin{array}{c}\hfill R_{\mathrm{BB}}=(\frac{278\mathrm{K}}{T_{\mathrm{dust}}}{)}^2(\frac{L}{L_{\odot }}{)}^{1/2}·\end{array}$$

An estimate of the “true” disc radius, R_disc, is then obtained by applying a stellar luminosity-dependent correction factor, Γ, which accounts for the radiation pressure blowout grain size

$$\begin{array}{c}\hfill \mathrm{\Gamma }=a{(L_{\ast }/L_{\odot })}^b,\end{array}$$

(Pawellek & Krivov 2015), using the new coefficients given in Pawellek (2016): a = 7.0 and b = −0.39. After applying this correction, the estimated disc size for HD 193571 is 120 ± 15 au. The disc has never been imaged in scattered light, and additional SPHERE/IRDIS observations were inconclusive in this respect (see Appendix A).

We used the fitted stellar spectrum to derive the stellar H and L′ magnitudes (reported in Table 1), integrating over the NaCo H- and L′-band filters. We used zero points of 1.139 × 10⁻¹⁰ erg cm⁻² s⁻¹ Å⁻¹ and 5.151 × 10⁻¹² erg cm⁻² s⁻¹ Å⁻¹, respectively².

3. Observations and data reduction

HD 193571 was observed at two different epochs with NaCo at the Very Large Telescope (VLT), and an additional third epoch was obtained with the Gemini Planet Imager (GPI, Macintosh et al. 2014) through the Fast Turnaround observing mode (Program ID: GS-2018A-FT-111).

3.1. VLT/NaCo

Coronagraphic ADI observations of HD 193571 were obtained in May 2016 and June 2018 in L′ band (see Table 2), making use of the Annular Groove Phase Mask (AGPM, Mawet et al. 2013) vector vortex coronagraph to suppress as much as possible the diffraction pattern from the host star. We used cube-mode, saving 100 frames per cube. The observations were interlaced with frequent sky observations for background subtraction (every ∼8 min) and bracketed with non-coronagraphic flux measurements to create an unsaturated point spread function (PSF) reference. The data was reduced with the ISPY end-to-end modular reduction pipeline GRAPHIC (Hagelberg et al. 2016). The main reduction steps comprise background subtraction, flat field correction, bad pixel cleaning, and centring. Each cosmetically reduced cube is then median combined. For a more detailed explanation on how the data reduction is performed we refer to the ISPY overview paper (Launhardt et al., in prep.). The observations are summarised in Table 2.

3.2. Gemini/GPI

HD 193571 was observed in the H band with GPI in coronagraphic ADI mode on 12 August 2018, obtaining 76 frames and achieving a total field rotation of 88°. The integration time for each exposure was 60 s.

The photometry of GPI data can be calibrated using the satellite spots, which are four reference spots created by diffraction of the central star light from a square grid superimposed on the pupil plane (Wang et al. 2014). They can be used to extract the photometry and spectroscopy of the central star. During the observations there was a misalignment of the grid that produces the satellite spots, resulting in a diffraction spike above two of the four satellite spots, thus rendering them unusable for photometric calibration. Therefore, in the following analysis, when referring to the satellite spots, we only refer to the two unbiased ones.

The data were reduced making use of the publicly available GPI Data Pipeline (Maire et al. 2010), with the following reduction steps:

• Calibration files were created using the “Dark” and “Wavelength Solution 2D” recipes, applied to the dark frame and the Argon lamp calibration snapshot taken as part of the observations;

• A bad pixel map was created combining the results of the “Hot Bad Pixel Map” and “Cold Bad Pixel Map” recipes, which have been applied respectively to a set of 15 dark frames and a set of 5 daytime Wollaston disperser flat frames for each filter (Y, J, H, K1, and K2). The calibration files were chosen from the Gemini Data Archive to be the closest in time to the observations;

• The data were reduced applying the “Calibrated Datacube Extraction” recipe, using the above-mentioned newly created calibration files. This recipe also includes an automatic search and characterisation of the four satellite spots, storing in the header the location and peak flux (in ADU) of all the spots, for each wavelength channel;

• The flux-calibrated cubes were oriented using the internal GPI recipe “Rotate North Up”.

4. Analysis and results

The final classically ADI reduced images for all the three epochs are shown in Fig. 1. A close-in companion is clearly visible in all three epochs south of the star.

4.1. Astrometry and photometry

To analyse the two NaCo datasets we used the ANDROMEDA (Cantalloube et al. 2015)³ package, which uses a maximum likelihood estimation approach together with negative fake signal injection to evaluate the astrometry and photometry of a companion in an ADI dataset. The algorithm needs as inputs the reduced frames (corrected for the AGPM throughput), the parallactic angles, and an unsaturated and exposure time-scaled image of the central star. Since we were interested in analysing only the known companion, we set the Inner Working Angle and Outer Working Angle keywords to 0.2 λ/D and 20 λ/D, respectively (we refer to Cantalloube et al. 2015 for a detailed explanation of the ANDROMEDA package). The final x and y offsets (and relative 3σ uncertainties) were converted into separation and position angle using a platescale for NaCo of 27.2 mas pix⁻¹, assuming a conservative error of 0.5 pixels on the centring of the frames, and correcting for the true north offset of $0\stackrel{°}{.}486\pm 0\stackrel{°}{.}180$ (Launhardt et al., in prep.). Given the target’s distance and L′ band magnitude (see Table 1), we converted the flux evaluated with ANDROMEDA, and relative 3σ uncertainties, into an absolute L′ magnitude for both epochs accounting for the uncertainties on the host star magnitude and distance from the system. The final astrometry and photometry values for the two NaCo epochs, as well as the GPI epoch, are given in Table 3.

For the GPI dataset we evaluated astrometry and photometry of the companion in a slightly different way since no unsaturated exposure of the central star was obtained.

For the astrometry, we made use of the satellite spots (visible in all the reduced frames) to create a PSF reference: we first averaged the two satellite spots in each frame, and then we averaged over the 76 frames, obtaining a PSF for each spectral channel. We use this PSF, together with the ANDROMEDA package, to obtain the astrometry of the companion (as was done for the NaCo datasets) in each spectral cube. The final astrometry is the weighted mean of the astrometric positions at each wavelength, and is given in Table 3 taking into account the GPI pixel scale of 14.166 mas pix⁻¹, the additional true north offset of 0.10 ± 0.13° as reported in Rosa et al. (2015), and a conservative error on the centring of 0.5 pixels.

To obtain the photometry of the companion we calibrated the cubes extracted in Sect. 3.2 in the following way:

• For each spectral channel, we averaged the satellite spots peak flux (stored in the header), obtaining a mean satellite flux in ADU, and relative standard deviation;

• We then converted the frame from ADU to physical units, using the following equation (as detailed on the GPI website⁴): $\mathrm{frame}[\mathrm{units}]=\frac{\mathrm{frame}[\mathrm{ADU}]}{\mathrm{Satellite}\mathrm{spectrum}[\mathrm{ADU}]}\times \frac{\mathrm{Star}\mathrm{Spectrum}[\mathrm{units}]}{\mathrm{Star}-\mathrm{to}-\mathrm{Satellite}\mathrm{Flux}\mathrm{ratio}}$ The “Star-to-Satellite Flux ratio” was calibrated by the GPI team⁵, and it is =(2 × 10⁻⁴)⁻¹. The “Star Spectrum” (in the desired flux units) is obtained from the stellar spectrum fitted in Sect. 2. We accounted for the uncertainty on the “Star-to-Satellite Flux ratio”, the uncertainties on the stellar spectrum, and the standard deviation of the satellite spots flux;

• To account for possible contamination from the stellar halo, we median combined all the frames in each spectral channel, and then subtracted this median from each photometrically calibrated cube;

• We then extracted a spectrum for the companion from each median-subtracted, photometrically calibrated cube, fitting a Gaussian to the companion to get the peak flux. The final spectrum is the weighted average of the spectra in all cubes.

The final spectrum of the companion is shown in Fig. 4. We integrated this spectrum over the NaCo H-band filter, obtaining a NaCo H-band apparent magnitude of 11.07 ± 0.06. This corresponds to an absolute magnitude of 6.89 ± 0.06. The final astrometry and photometry for the companion is given in Table 3.

The close separation makes it unlikely for the companion to be a background star. Nevertheless, we evaluated the position that the object would have on the sky at epoch 2018, starting from its position in epoch 2016, if it were a background object with no significant proper motion. The results are shown in Fig. 3. The object is clearly co-moving with the host star, at a projected separation of ∼11 au.

Fig. 3. Proper motion analysis of the companion showing the astrometry for the three epochs. The black data point is the position that the companion would have at the epoch of the GPI observation if it were a background star with no motion, using its position in 2016 as starting point and considering the proper motion of the host star. The companion is clearly co-moving with the star (shown in yellow).

4.2. Physical properties

We compared the GPI H-band spectrum with observed spectra of early M dwarfs from the stellar spectral library⁶ of the CARMENES survey (Reiners et al. 2018), which is the first large library of M dwarfs with high-resolution spectra in the infrared. We plot three of the best matching spectra (binned to the GPI H-band resolution) in Fig. 4, a non-matching spectrum (dotted grey line) for comparison, and the H-band spectrum of HD 193571 B. From the comparison, we can infer a surface gravity of log g ∼ 4.9, a temperature of ∼3500 K, and a spectral type between M3 and M2, which seem to fit the data reasonably well. However, a high-resolution and/or broader band spectrum would be needed to properly constrain the surface gravity and spectral type of the companion.

Fig. 4. Comparison between the spectrum of the companion and observed spectra of early M dwarfs. The blue shaded area is the flux density of the companion in the GPI H-band, in Jansky. The spectrum is the weighted average of the spectra extracted from the 76 GPI datacubes and the area encompass the uncertainties (derived from the uncertainty on the flux of the host star). The solid lines are three spectra from the CARMENES stellar spectral library, for various Teff and log g values (evaluated in Hintz et al. 2019) and the dotted grey line is an additional spectrum of an M1 object.

We estimated the mass of the companion using the BT-Settl evolutionary tracks (Allard et al. 2012)⁷, by comparing them with the observed L′- and H-band photometry. In the colour-magnitude diagram of Fig. 5 we show the companion L′-band absolute photometry of 6.19 ± 0.08 mag (evaluated as the weighted mean of the two NaCo epochs), as well as evolutionary tracks for two representative ages of 60 Myr (dashed line) and 150 Myr (solid line). As shown in Fig. 5, the photometry does not allow us to distinguish between the two age estimates, so we use both age values in the rest of the paper. We interpolated the BT-Settl models to estimate the mass of the companion for both L′- and H-band photometry, in mass steps of 0.034 dex. Taking into account the photometric uncertainty in both bands, we obtained a weighted mass of 0.395 ± 0.007 M_⊙ for an age of 161 Myr, and 0.305 ± 0.025 M_⊙ for an age of 66 Myr.

Fig. 5. Colour-magnitude diagram showing the weighted mean L′-band magnitude derived from the 2016 and 2018 NaCo datasets, together with the H-band magnitude derived from the GPI dataset. We plot the evolutionary tracks for the BT-Settl models from Allard et al. (2012), for ages of 60 and 150 Myr. The photometry does not allow us to distinguish between the two age estimates.

4.3. Orbital motion

The astrometry of the companion between the three epochs shows signs of orbital motion. Following the prescription in Pearce et al. (2015), we can explore the possible orbital solutions for a companion imaged over a short orbital arc, using the dimensionless parameter B ($\sqrt{B}=V_{\mathrm{sky}}/V_{\mathrm{esc}}$ is the sky-plane velocity of the companion divided by the escape velocity), and the direction of motion φ, where φ = 0° is motion along a vector from the primary to the companion.

We assumed a total system mass of 2.6 ± 0.1 M_⊙ (for an age of 161 Myr) and 2.55 ± 0.1 M_⊙ (for an age of 66 Myr) and we derived⁸B and φ for the three epochs (NaCo 2016, NaCo 2018, and GPI 2018). For both age estimates the values agree within the uncertainties, and we obtain $B=0.{25}_{-0.11}^{+0.16}$ and φ = 100 ± 15°, which leads to a minimum semi-major axis of a_min = 8.20 ± 1.77 au (see Eq. (5) in Pearce et al. 2015). Following Pearce et al. (2015), we can draw the following conclusions:

• Even considering the uncertainties, the B value is < 1, so the companion’s sky-plane motion is below the escape velocity. While the object could be unbound if the line of sight velocity (or separation) is high, this is unlikely;

• We cannot place constraints on the eccentricity of the orbit, meaning that a circular orbit cannot be ruled out (this will have an impact on our stirring mechanisms study in Sect. 4);

• We can place a loose upper limit of ∼80° on the inclination.

We also explored the possible orbital motion parameters using the python package orbitize⁹ with the Orbit For The Impatient (OFTI) algorithm detailed in Blunt et al. (2017) (see Appendix C). While the uncertainties on the astrometry and the limited amount of datapoints do not place any meaningful constraints on the orbital elements, the periastron distance is restricted to ≲15 au. This result is confirmed by exploring the possible orbital parameters using the method of Pearce et al. (2015). Therefore, if the companion’s orbit is nearly coplanar with the disc, the entire orbit should be interior to the disc, otherwise the companion would have disrupted the disc on a dynamical timescale. Assuming a circular orbit and a semi-major axis of 11 au, the companion would have a minimum period of ∼23 years, implying that a baseline of several years would be needed before any additional astrometric datapoint could provide better constraints on the orbital elements. The companion is massive enough that even in the unlucky case of an almost face-on orbit (i ∼ 1°) it would produce a radial velocity signal strong enough to be detected (semi-amplitude K ≳ 120 m s⁻¹); however, this would also require a time baseline of many years.

5. Stirring mechanisms

The relative importance of self- and companion-stirring mechanisms is a non-trivial problem. It depends on the companion’s physical and orbital parameters, the host star age and mass, and the disc mass in solids. The equations used in this section are from Wyatt (2008) and Mustill & Wyatt (2009), and are summarised in Appendix B. We note that to be consistent with the underlying assumptions of these two papers, we use the black-body disc radius of 62 au while working with equations from Wyatt (2008), and the corrected disc radius of 120 au for the Mustill & Wyatt (2009) equations (see Appendix B). That is, the model in Wyatt (2008) uses parameters derived by fitting to black-body radii, while the model of Mustill & Wyatt (2009) uses orbital dynamics, so is based on physical disc radii.

Assuming that the mutual inclination between the plane of the orbit and the disc is not too large, there are two conditions that need to be satisfied for a companion to dominate the stirring process at a certain distance from the star, and at a given time: (a) the companion must be able to stir planetesimals, at that location, to relative destructive velocities and (b) the timescale for companion-stirring at that distance must be greater than the self-stirring timescale.

The first condition is encapsulated by Eqs. (B.2) and (B.3), which give the maximum distance at which a companion with a given semi-major axis a_pl and eccentricity e_pl can stir planetesimals above the disruption threshold velocity v_rel. This velocity is a function of the planetesimal size R and, as shown by Eq. (B.2), has a minimum at R ∼ 80 m. We set this maximum distance equal to the estimated true disc radius of 120 au, and we plotted the a_pl–e_pl relationship in Fig. 6 for the R = 80 m case (solid light blue curve). The companion would not be able to stir planetesimals at that distance if its semi-major axis and eccentricity were below this curve. The planetesimals might be smaller or larger than 80 m, and this would increase v_rel and push the light blue curve rightwards and upwards. While R has a definite minimum (particles smaller than a certain size, typically around few μm, would be blown away by radiation pressure from the central star) it is not straightforward to define a maximum R value. We proceeded as follows:

• At any given time, there is a maximum size of planetesimals that participate in the collisional cascade (because larger objects will have collision timescales longer than the stellar age). This maximum size R_max can be evaluated by inverting Eq. (B.1). For a disc size of 62 au, and with a fractional luminosity of the disc f, stellar mass and stellar luminosity as in Table 1, we have R_max = 132 m. This is the maximum value for R, assuming that the disc has been stirred for all of its life (t_stir = t_age = 66 Myr. In the 161 Myr case we obtain R_max = 790 m);

• An internal perturber can influence the timescale of orbit crossings for planetesimals, and thus t_stir might be less than the stellar age (i.e. the disc was stirred more recently). We use Eq. (B.4) to calculate this orbit crossing timescale t_cross as a function of the perturber properties (eccentricity, semi-major axis, and mass);

• We now have a revised value for the total time the disc has been stirred as t_stir = t_age − t_cross, and consequently a revised R_max value as a function of the perturber properties (i.e. we have a relationship between R_max, a_pl, and e_pl);

• Combining this relationship with Eqs. (B.2) and (B.3), we can trace R_max in the (a_pl, e_pl) parameter space.

Fig. 6. Boundaries between a self-stirring and companion-stirring dominated disc. The light blue lines mark the (apl, epl) parameter space in which the companion would be able to stir planetesimals of size R to destruction velocities at a distance of 120 au. The shaded area around the solid light blue (R = 80 m) line takes into account the errors on the disc size and the stellar mass. The dashed purple line shows the Rmax for 66 Myr (close to the solid light blue line) and the dashed green line shows the Rmax value for the 161 Myr case. The shaded red areas indicate the boundaries between the self-stirring and companion-stirring dominated cases, for a fixed distance and companion mass, and for two representative xm values; accounting for errors on disc size, stellar mass, and companion mass (the areas encompass both age estimates). The horizontal dotted black line is the lowermost boundary of the minimum possible companion semi-major axis calculated in Sect. 4.3. The companion dominates the stirring process only for combinations of apl and epl lying above the light blue curve (the companion can stir planetesimals at the disc distance) and the red curve (the companion stirs the disc faster than the disc stirs itself).

As can be seen in Fig. 6, when we plot this for the 66 Myr case, R_max is relatively small (∼132 m along the curve) and almost overlaps with the R = 80 m case. The R_max in the 161 Myr case is plotted with a dashed grey curve. The companion can stir the disc over most of the shown parameter space.

The second condition requires that, at a given time and distance, the companion-stirring timescale is shorter than the self-stirring timescale. Mustill & Wyatt (2009) made such a study and defined the parameter Φ as the distance at which self and companion-stirring times are equal (see Appendix B). It is a function of the companion’s properties (mass m_pl, semi-major axis a_pl, and eccentricity e_pl), the central star’s mass, and the disc’s mass in solids (expressed by the dimensionless parameter x_m, see Appendix B). Since we are interested in which stirring process is dominant at the location of the debris belt, we set Φ = 120 au and obtain the equilibrium relationship between self- and planetary-stirring. Tracing this line in the (a_pl, e_pl) parameter space marks the boundary between the domination of the two stirring processes, thus allowing us to investigate the combination of a_pl and e_pl for which the disc is dominated by self-stirring. Since there is a dependence on the x_m value as well, in Fig. 6 we plotted two representative values for x_m of 1 and 10 (solid red lines). The curve for x_m = 10 lies above the x_m = 1 case because a more massive disc forms large planetesimals more quickly, and can thus self-stir earlier. As discussed in Mustill & Wyatt (2009), x_m ≳ 10 discs may be problematic as their high masses imply gravitationally unstable discs at earlier times when the gas was present. Thus, it is likely that the x_m = 10 line in Fig. 6 represents an upper limit to where the disc could be self-stirred. Given an x_m value and fixing the companion mass to 0.25 M_⊙, any combination of eccentricity and semi-major axis above the curve would imply that companion-stirring is quicker than self-stirring at the distance of the disc, hence the companion-stirring would dominate the stirring process. An additional constraint can be placed on the minimum semi-major axis, as discussed in Sect. 4.3, which is shown by the dashed black line in Fig. 6.

It is important to note that both conditions must be satisfied for the companion to dominate the stirring process, and this is true only for certain combinations of eccentricity and semi-major axis. In the plot it is clear how, given an eccentricity ≳0.1, any semi-major axis places the companion above both curves, and thus the companion would dominate. For eccentricities ≳0.002, any a_pl would lie above the light blue curves (both for the R ∼ 80 m and for the R_max case), but only certain a_pl would satisfy the criterion for companion-induced stirring dominating over self-stirring (depending on the x_m value), so low-eccentricity companions must be closer to the disc to dominate the stirring. Finally, for extremely low eccentricities (≲0.002) and small semi-major axes, the companion would not be able to stir planetesimals at the distance of the disc (below the light blue curve), and in any case the self-stirring would be dominant at that distance (below the red curve).

As shown in Fig. 6, it is most likely that the companion is dominating the stirring process, and self-stirring is relevant only when the companion has a very low eccentricity (in combination with a small semi-major axis).

6. Conclusions

We presented the first detection of a close low-mass stellar companion around the A0 star HD 193571. The three epochs obtained with VLT/NaCo and GPI confirm that the companion is co-moving with the host star, showing the potential of multi-band/multi-instrument follow-up to confirm direct imaging candidates. Comparing M_H and M_L′ band photometry to evolutionary tracks suggests a mass of ∼0.305 ± 0.025 M_⊙ for an age of 66 Myr (∼0.395 ± 0.007 M_⊙ for the 161 Myr case), which would make it an M2-2.5 dwarf. Comparison to observed spectra seems to suggest a surface gravity of ∼4.9 and a temperature of ∼3500 K. The orbital motion detected in the three epochs is not enough to place solid constraints on the orbital parameters, but allows us to confirm the co-motion with the host star and to exclude an edge-on orbit.

Given the projected separation of ∼11 au and a maximum periastron of ∼15 au, the companion appears to orbit interior to the circumstellar debris belt (inferred via SED IR-excess to be at ∼120 au). We investigated the plausibility that both self- and companion-stirring mechanisms are responsible for the currently observed debris belt radius. Since no constraints can be put on the eccentricity, we cannot exclude a fully self-stirring scenario for the disc. However, a small deviation from a circular orbit would result in the disc being dominated by companion-stirring (as shown in Fig. 6) and if the orbit is sufficiently eccentric the disc will appear eccentric as well. The companion is likely responsible for the stirring of a disc that appears to be an order of magnitude further away, showing how a massive companion can influence a debris disc at large distances.

At the moment, only a handful of systems are suited for a study of stirring mechanisms, and the HD 193571 system represents an important addition, containing the third known M-dwarf companion to a young star discovered to be orbiting within the primary’s circumstellar disc, and the first one found around an A0-type star. In the future, radial velocity observations as well as a resolved image of the disc could be useful in deepening our understanding of this system.

Acknowledgments

A.M.B. thanks Christian Ginsky for the help in reducing the SPHERE/IRDIS data and the useful discussion afterwards, and Wolfgang Brandner for the help in reducing the GPI data. A.M.B. would also like to thank the anonymous referee for the useful and constructive comments. G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. G.C. and S.P.Q. thank the Swiss National Science Foundation for the financial support under the grant number 200021_169131. J.O. and N.G. acknowledge financial support from the ICM (Iniciativa Científica Milenio) via the Núcleo Milenio de Formación Planetaria grant. J.O. acknowledges financial support from the Universidad de Valparaíso, and from Fondecyt (grant 1180395). N.G. acknowledges support from project CONICYT-PFCHA/Doctorado Nacional/2017 folio 21170650. A.M. and A.Q. acknowledge the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (MU 4172/1-1 and QU 113/6-1). This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory, and of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. 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. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: IRDIS disc non-detection

We observed the target with SPHERE/IRDIS at the VLT in coronagraphic Differential Polarisation Imaging (DPI) mode on 26 September 2018, using the H broad-band filter.

We took eight polarimetric cycles, each consisting of four data cubes, one per half wave plate (HWP) position. Each data cube consisted of four individual exposures with exposure times of 32 s. The science observations were bracketed with 2 s exposures, to create an unsaturated PSF reference for the central star.

The data were reduced following the prescription in Ginski et al. (2016), obtaining the radial Stokes components Q_Φ and U_Φ (see Schmid et al. 2006), where Q_Φ would contain any polarisation signal coming from dust scattered light, and it is shown in Fig. B.1. No emission is visible at the expected location of the disc ($\sim 1\stackrel{″}{.}75$) or anywhere else. The faint emission from the centre is due to the stellar halo, and the spider is vaguely visible extending approximately in the north-south direction.

Appendix B: Stirring mechanisms

B.1. Self-stirring

From Wyatt (2008) the maximum fractional luminosity f_max of a planetesimal belt at distance r around a star of mass m_⋆, luminosity L_⋆, and age t_age is

$$\begin{array}{c}\hfill f_{\max }=0.58\times {10}^{-9}{r}^{7/3}(\mathrm{d}r/r)R_{\max }^{0.5}{Q}_{\mathrm{D}}^{\star 5/6}{e}^{-5/3}{m}_{\star }^{-5/6}{L}_{\star }^{-0.5}{t}_{\mathrm{age}}^{-1},\end{array}$$(B.1)

where R_max is the maximum size of the planetesimals that are participating in the cascade at that given time (called D_c in Wyatt 2008), $Q_{\text{D}}^{\star }$ is the planetesimal strength in Jkg⁻¹, e is the mean planetesimal eccentricity, and dr/r is the relative width of the planetesimal belt. It was found (see Wyatt 2008) that the population of debris discs around A stars can be fitted assuming $Q_{\text{D}}^{\star }=150\text{J}{\text{kg}}^{\text{-1}}$, e = 0.05, and dr/r = 0.5. All of this assumes that the disc has been stirred for its whole lifetime (i.e. t_stir = t_age). The disc evolution model developed in Wyatt (2008) is SED-based, and therefore the planetesimal belt distance r refers to the black-body radius R_BB of 62 au, inferred via SED fitting.

Fig. B.1. DPI data taken with SPHERE/IRDIS, with a total field of view of ∼10″ × 10″, given a pixel scale for IRDIS of 12.25 mas pix−1. No polarised signal from the disc scattered light is visible. The image is oriented with north up and east left.

B.2. Companion-stirring

From Mustill & Wyatt (2009), the threshold velocity above which collisions between planetesimal of size R become destructive is

$$\begin{array}{c}\hfill v_{\mathrm{rel}}^{\star }(R)={[0.8{\left(\frac{R}{80\mathrm{m}}\right)}^{-0.33}+0.2{\left(\frac{R}{80\mathrm{m}}\right)}^{1.2}]}^{0.83}\mathrm{m}{s}^{-1}.\end{array}$$(B.2)

A companion of mass m_pl internal to the disc on an orbit of semi-major axis a_pl and eccentricity e_pl, around a primary of mass m_⋆, would be able to stir planetesimals to catastrophic collisions only up to a maximum distance a^⋆:

$$\begin{array}{c}\hfill a^{\star }(R)=3.8\mathrm{au}{\left(\frac{e_{\mathrm{pl}}}{0.1}\right)}^{2/3}{\left(\frac{m_{\star }}{1M_{\odot }}\right)}^{1/3}{\left(\frac{a_{\mathrm{pl}}}{1\mathrm{au}}\right)}^{2/3}{\left(\frac{v_{\mathrm{rel}}^{\star }(R)}{1\mathrm{k}m{s}^{-1}}\right)}^{-2/3}·\end{array}$$(B.3)

In addition, it is possible to calculate the timescale for orbit crossing of planetesimals at a distance a as

$$\begin{array}{cc}\hfill t_{\mathrm{cross}}\sim & 1.53\times {10}^3\frac{{(1-e_{\mathrm{pl}}^2)}^{3/2}}{e_{\mathrm{pl}}}{\left(\frac{a}{10\mathrm{au}}\right)}^{9/2}\hfill \\ \hfill & \times {\left(\frac{m_{\star }}{M_{\odot }}\right)}^{1/2}{\left(\frac{m_{\mathrm{pl}}}{M_{\odot }}\right)}^{-1}{\left(\frac{a_{\mathrm{pl}}}{1\mathrm{au}}\right)}^{-3}\mathrm{yr}.\hfill \end{array}$$(B.4)

B.3. Companion-stirring versus self-stirring

Mustill & Wyatt (2009) also defined the parameter Φ as the distance boundary between self-stirring and companion-stirring at a fixed age as

$$\begin{array}{cc}\hfill \mathrm{\Phi }=& 630\mathrm{au}{(1-e_{\mathrm{pl}}^2)}^{-1}{e}_{\mathrm{pl}}^{2/3}{\left(\frac{m_{\mathrm{pl}}}{M_{\odot }}\right)}^{2/3}\hfill \\ \hfill & \times {\left(\frac{a_{\mathrm{pl}}}{1\mathrm{au}}\right)}^2{\left(\frac{m_{\star }}{M_{\odot }}\right)}^{-4/3}{x}_{\mathrm{m}}^{-0.77},\hfill \end{array}$$(B.5)

where the dimensionless parameter x_m is a scaling factor relating the disc surface density to the minimum mass solar nebula density (see Mustill & Wyatt 2009 and Kenyon & Bromley 2008).

The model developed in Mustill & Wyatt (2009) is a dynamic model that depends on the physical structure of the disc, and therefore on the real disc size of 120 au (see Sect. 2).

Appendix C: Orbital constraints with OFTI

We explored the possible orbital motion parameters using the python package orbitize with the Orbit For The Impatient (OFTI) algorithm detailed in Blunt et al. (2017). We used two total mass estimates: 2.6 ± 0.1 M_⊙ (for an age of 161 Myr) and 2.5 ± 0.1 M_⊙ (for an age of 66 Myr). We used a uniform prior for the semi-major axis, and in the epoch of periastron passage and argument of periastron. We used a sin(i) prior for the inclination angle, and a linearly descending prior for the eccentricity, with a slope of −2.18. For both age estimates, the results agree within the error bar, and in Fig. C.1 we show the posterior distribution function for the 161 Myr case. As shown in the figure, the uncertainties on the astrometry and the limited number of datapoints do not allow us to place any meaningful constraints on the orbital elements, but the periastron distance q is restricted to ≲15 au.

Fig. C.1. Posterior distribution function for the orbital parameters derived with the orbitize package using the OFTI implementation.

All Tables

All Figures

Fig. 1. Classically ADI reduced images for the two NaCo datasets (left and centre) and for the GPI dataset (right). The images are oriented with north up and east left, and the green cross indicates the position of the central star. The companion is clearly visible close to the centre in all three datasets. The images are normalised and the colour map was chosen for a better visualisation of the data.

In the text

	Fig. 2. Flux density distribution of HD 193571, showing the photometric datapoints found in the literature (in blue) and the IRS spectrum (in black), together with the fitted stellar (green) and disc (red) fluxes.
In the text

Fig. 3. Proper motion analysis of the companion showing the astrometry for the three epochs. The black data point is the position that the companion would have at the epoch of the GPI observation if it were a background star with no motion, using its position in 2016 as starting point and considering the proper motion of the host star. The companion is clearly co-moving with the star (shown in yellow).

In the text

Fig. 4. Comparison between the spectrum of the companion and observed spectra of early M dwarfs. The blue shaded area is the flux density of the companion in the GPI H-band, in Jansky. The spectrum is the weighted average of the spectra extracted from the 76 GPI datacubes and the area encompass the uncertainties (derived from the uncertainty on the flux of the host star). The solid lines are three spectra from the CARMENES stellar spectral library, for various Teff and log g values (evaluated in Hintz et al. 2019) and the dotted grey line is an additional spectrum of an M1 object.

In the text

	Fig. 5. Colour-magnitude diagram showing the weighted mean L′-band magnitude derived from the 2016 and 2018 NaCo datasets, together with the H-band magnitude derived from the GPI dataset. We plot the evolutionary tracks for the BT-Settl models from Allard et al. (2012), for ages of 60 and 150 Myr. The photometry does not allow us to distinguish between the two age estimates.
In the text

Fig. 6. Boundaries between a self-stirring and companion-stirring dominated disc. The light blue lines mark the (apl, epl) parameter space in which the companion would be able to stir planetesimals of size R to destruction velocities at a distance of 120 au. The shaded area around the solid light blue (R = 80 m) line takes into account the errors on the disc size and the stellar mass. The dashed purple line shows the Rmax for 66 Myr (close to the solid light blue line) and the dashed green line shows the Rmax value for the 161 Myr case. The shaded red areas indicate the boundaries between the self-stirring and companion-stirring dominated cases, for a fixed distance and companion mass, and for two representative xm values; accounting for errors on disc size, stellar mass, and companion mass (the areas encompass both age estimates). The horizontal dotted black line is the lowermost boundary of the minimum possible companion semi-major axis calculated in Sect. 4.3. The companion dominates the stirring process only for combinations of apl and epl lying above the light blue curve (the companion can stir planetesimals at the disc distance) and the red curve (the companion stirs the disc faster than the disc stirs itself).

In the text

	Fig. B.1. DPI data taken with SPHERE/IRDIS, with a total field of view of ∼10″ × 10″, given a pixel scale for IRDIS of 12.25 mas pix−1. No polarised signal from the disc scattered light is visible. The image is oriented with north up and east left.
In the text

	Fig. C.1. Posterior distribution function for the orbital parameters derived with the orbitize package using the OFTI implementation.
In the text