© The Authors 2024

1 Introduction

During the early stages of star formation, a significant fraction of stars are found to be part of multiple systems (Reipurth et al. 2014; Offner et al. 2023). Multiple systems naturally evolve to hierarchical configurations (Duchêne & Kraus 2013), but allow a large diversity of possible configurations for protoplanetary discs. Circumstellar discs can in principle form around any hierarchical level of the system. As a consequence, discs in multiple systems are shaped by the disc-disc and star-disc interactions that depend on the geometrical configuration and physical properties of the system (Bate 2018). By means of those interactions, a large cavity is carved in circumbinary discs up to several times the binary semi-major axis, while the outer edge of the disc is set by tidal interactions with the outer stars of the system (Artymowicz & Lubow 1994; Miranda & Lai 2015). Spiral arms are usually launched in the process (Rafikov 2002). Gravitational torques from misaligned stars tend to incline discs and to make them precess (Papaloizou & Terquem 1995), resulting in misaligned geometries of the disc plane with respect to the stellar orbital plane and even leading to a disc warp or break in some cases (e.g. Nixon et al. 2013; Facchini et al. 2013; Rabago et al. 2023). The investigation of these dynamical behaviours is key to understand disc dynamics in multiple systems and their resulting planetary architectures.

Hydrodynamical simulations have proven to be a powerful tool to infer orbits from disc morphology and kinematics in multiple systems. Multiplicity-induced substructures depend on the orbital parameters of companions. By running a grid of hydro-dynamical models and comparing them to disc observations, one can now constrain orbits in multiple systems (Price et al. 2018a; Gonzalez et al. 2020; Nowak et al. 2024). Even if this method is limited by observations in snapshots in time and the cost of hydrodynamical simulations, it provides constrains independent from astrometric measurements, and it can help to discriminate between available orbits.

V892 Tau is a young triple system located in a clustered sub-region of the Taurus star-forming cloud (Fig. 1). The system is composed of two central 3 M_⊙ stars separated by 0.06″ surrounded by a large circumbinary disc (CBD), which extends up to approximately 1.5″ from the central stars (e.g. Monnier et al. 2008; Long et al. 2021). The M3 star V892 Tau NE is projected at 4″ north-east from the inner binary and completes the system (Esplin & Luhman 2019). The Gaia DR3 (Gaia Collaboration 2023) parallax measurements result in distances of d = 134.5 ± 1.5 pc and d = 131.3 ± 5.0 pc, respectively, for V892 Tau and V892 Tau NE, which is consistent with a bound triple system within the error bars.

Observations with the Atacama Large Millime-ter/submillimeter Array (ALMA) have unveiled the dust ring surrounding the inner binary of V892 Tau at a high spatial resolution (Pinilla et al. 2018). More recently, a detailed study of the system has been conducted by Long et al. (2021, hereafter L21) using that ALMA Band 6 data combined with additional Very Large Array (VLA) observations: The eccentric inner binary has been resolved and its orbit tightly constrained, with two families of solutions for the longitude of the ascending node Ω_in and the argument of the periapsis ω_in. The dust emission was found to have an 0.17″ -large inner gap, while the gas emission appears to extend to the innermost regions of the system. Both the dust and gas show hints of an inclined geometry with respect to the inner binary. The tilt of the disc with respect to the inner binary is measured at Δ = 8.0 ± 4.2° or Δ = 113.3 ± 3.0° depending of the orbital solution chosen. The size of the gas disc is consistent with a tidal interaction with V892 Tau NE and tentative spiral arms are distinguished in the redshifted side of the disc. A warp is suggested by deviations from Keplerian rotation in the outer disc and could be caused by interactions between the disc and the inner or outer binary. In addition, recent near-IR interferometric observations have refined the orbit of the binary and detected a narrow circumstellar disc wide of 0.016″ around the primary star of V892 Tau (Vides et al. 2023). In that work, the inner CBD was also modelled and tentative evidence of a warp was invoked to explain residual patterns from geometric models.

As a hierarchical triple system, V892 Tau is a perfect laboratory to study the formation channels of planets from discs in multiple systems, which depends on the binary-disc interactions and dynamics. Understanding the disc dynamics in such a system will help to constrain the planet formation regions in multiple stellar systems in general, but it would require the orbital configuration of the system to be known. Nonetheless, the orbit of the companion star V892 Tau NE remains unconstrained at the moment. Projected at a distance of ~520 au of the central binary, the orbit of V892 Tau NE is critical to understand the dynamics of the CBD and the subsequent planet formation.

Following that effort, here we present new ALMA observations of V892 Tau that we combined with archival data. The observational setup and the resulting data are detailed in Sect. 2. We fitted the observations with a Keplerian disc model before modelling the V892 Tau system as a triple star system with a CBD. We performed hydrodynamical simulations of that model, testing different orbits for the outer companion, and post-processed the simulation outputs to build synthetic observations. Methods are presented in Sect. 3, while an analysis of the observations, of the synthetic observations, and their relative comparison are described in Sect. 4. In Sect. 5, we give an overview of the system and discuss our results in the more general context of multiple systems. Our conclusions are summarised in Sect. 6.

Fig. 1 V892 Tau is located in a clustered sub-region of the Taurus molecular cloud. The background image from the Herschel/SPIRE survey (Griffin et al. 2010; André et al. 2010) is shown on a linear intensity scale. Magenta points show positions of pre-main-sequence stars in Taurus taken from Kenyon et al. (2008).

2 Observations

2.1 Data calibration, reduction, and combination

The triple system V892 Tau was observed with ALMA in the context of the ALMA programme 2021.1.01137.S (PI: J. Miley) in Band 6. The observations were performed in two sessions, with the first session covered shorter baselines spanning 15–1997 m (12–1536 kλ) with 42 antennas. We call the resulting dataset of that session Short Baselines (SB hereafter). The second session used a more extended configuration also with 42 antennas spanning baselines from 64 m to 8283 m (49–6371 kλ) and resulted in a dataset we call Long Baselines (LB hereafter). Each dataset contains four spectral windows (SPWs) with two of them set up for continuum detection centred at 217.875 and 233.000 GHz, with 2 GHz of bandwidth distributed over 128 channels. The other two SPWs were centred at 230.538 and 220.000 GHz targeting the ¹²CO (2–1) and ¹³CO (2–1) transitions, each containing 3840 channels of 122 kHz and 488 kHz widths, respectively.

The raw visibility data were downloaded from the ALMA archive and subsequently calibrated using the available scripts provided by ALMA staff using the required CASA (McMullin et al. 2007) 6.2.1.7 version. The band-pass and amplitude calibrator for the observations was J0435 + 2532, and J0438 + 3004 was used as the phase calibrator. The total on-source time reached approximately 18.5 min.

The data were self-calibrated after the line-free channels were combined altogether to create continuum datasets. We performed several phase self-calibration rounds on each continuum dataset until the solution interval was tuned down to the integration time of 6.05 s for the SB dataset and of 2.02 s for the LB dataset. This process resulted in an increase of the peak signal-to-noise ratio (S/N) of about 6% for SB and 36% for LB, giving a peak S/N of 582 and 233, respectively.

We downloaded the data used in L21 from the ALMA archive and re-calibrated it thanks to the ALMA reduction pipeline. The resulting dataset is called the Original Long et. al dataset (OLD hereafter). Then, we reduced the data and performed phase self-calibration rounds on the continuum data with solution intervals down to the integration time of 6.05 s. The peak S/N increased by 9% doing so to reach a value of 483, which is a significant improvement given the already high S/N of the data.

The SB and LB visibility amplitudes were scaled with respect to the one of OLD to make sure the fluxes of all the datasets were consistent altogether before combination. To do so, we used the functions estimate_scale_flux and rescale_flux of the reduction_utils.py python script from the DSHARP project¹. We then created a combined continuum model with tclean in CASA v6.5.3 from the three datasets, using a Briggs weighting with a robust parameter of 0.5 (Briggs 1995). On the basis of that common model, the individual datasets were self-calibrated in phase a second time. The resulting datasets were combined to create a final continuum image with tclean using a Briggs weighting with a robust of −0.5, which resulted in a RMS of 90 µJy. The resulting beam size is 0.11″ × 0.06″ at a central frequency of 224 GHz. Additional details regarding the data combination process are available in Appendix A.

We also applied the resulting self-calibration tables to the gas observations, from which the continuum emission was subtracted by fitting a polynomial function to line-free channels using the uvcontsub task. An image cube of the ¹²CO (2–1) line emission was then created combining the three datasets using a Briggs weighting with a robust of 0.5 and a channel width of 0.5 km s⁻¹, giving a resulting beam size of 0.17″ × 0.11″. We summarise the details of each of the continuum and ¹²CO (2–1) images in Table 1. This paper presents the continuum observations but mainly focusses on the analysis and discussion of the gas emission, as a detailed description of the continuum data will be addressed in a future study.

2.2 Dust disc

Figure 2 shows the ALMA 1 .3 mm continuum image and a deprojected brightness radial profile. The observed structure emerges as a Gaussian-like ring located at a separation of approximately 0.22″ from the centre of the disc, following the position of the peak intensity radius. The inner binary has carved a cavity in the innermost parts of the system. We found a cavity radius of R_cav = 0.14 ± 0.02″ that is defined as the radius at which the intensity first reaches half of its peak level with an incertitude taken as the difference with quarter peak radius. The dusty disc has a deprojected radius of R_90% = 0.38 ± 0.01″ measured as the radius encircling 90% of the total flux.

The inclination i and the position angle PA of the disc are respectively best-fitted to i = 54.0 ± 1.8° and PA= 51.5 ± 0.5°. Thanks to an improved angular resolution, the 1.3 mm emission is found to be more compact than in previous observations (Pinilla et al. 2018). The cavity is spatially resolved along with an azimuthal asymmetry in the NW side of the disc. Indeed the northern side of the disc along the minor axis is 18% brighter than its southern side. This trend was also found in L21 and is consistent with the reported value of 20%. This kind of horseshoe asymmetry is a known consequence of binarity (Ragusa et al. 2017). However optically thick warm dust emission from the inner rim of the disc may also produce similar patterns (Ribas et al. 2024). This emission coming from the far side of the disc, it would indicate that the closest side to the observer is the SE side of the disc. Like all the previous ALMA observations of V892 Tau, no emission is detected from the vicinity of V892 Tau NE in the 1.3 mm continuum. The sensitivity and RMS of the observations presented in this paper being similar to the previous data, this result is not surprising.

We confirm the detection of unresolved emission in the cavity with a 30σ significance. This emission could be explained by circumstellar dusty material and could trace the circumstel-lar disc detected in the inner system (Vides et al. 2023). Future multi-wavelengths observations could allow a precise characterisation of the dust properties, while multi-epochs observations could allow the clump’s dynamical behaviour to be prescribed. Those questions are left for future investigations.

2.3 CO emission

The ¹²CO emission of the disc is detected at a confidence level above 3σ from channel −6 km s⁻¹ to channel 23 km s⁻¹ in the kinematic local standard of rest (LSRK) velocity. The data suffer from absorption by cloud material at the systemic velocity around 8 km s⁻¹. The complete channel maps can be found in Appendix B. 99% of the flux is contained in a radius of 1.73″, with a maximum recoverable scale of the observational setup estimated at 2.27″ using the 5th percentile of the shortest baseline lengths. The gas disc has a radius of R_90% = 1.45″ and a flux of 15.8 Jy km s⁻¹, which is in good agreement with L21. The north-east side of the disc appears blueshifted compared to the systemic velocity of ~8 km s⁻¹, which means the disc rotates in an anti-clockwise way. The morphology of the emission in the velocity channels is similar to the observations reported in L21: the faint lower emission surface of the disc lacks clear visibility, indicating potential blending with the upper emission surface due to spatial resolution constraints or possible non-detection. This makes the true orientation of the disc ambiguous and the closest side to the observer difficult to define. Irregularities at the edge of the redshifted side are recovered as well as in the northern side of the disc. However these irregularities are not found to be part of larger structures that could have traced potential interactions with the companion star V892 Tau NE. In the following, we model V892 Tau as a triple system and try to constrain the orbit of V892 Tau NE based on hydrodynamical simulations. Our methods to model the ¹²CO emission in detail and the hydro-dynamical setups used are described in Sect. 3 below. From this modelling, we discuss in more detail the observations in Sect. 4.

Fig. 2 ALMA band 6 continuum emission map (top) and depro-jected azimuthally averaged intensity profile (bottom) at 224 GHz of V892 Tau. In the top image, contour levels are at 10σ, 50σ, and 150σ. In the bottom plot Rcav, Rpeak, and R90,cort correspond to the radius of the cavity, of the peak of the ring, and of the disc, respectively. The shaded orange area represents to the standard deviation of each annulus divided by the square root of the beam numbers along the annulus at each radial bin. The shaded grey area represents the major axis of the synthesised beam.

3 Methods

3.1 Discminer

In order to understand at a deeper level the kinematical information contained in the structured disc of V892 Tau, we built a quantitative Keplerian model to be compared with the observations. To do so, we used the python package DISCMINER (Izquierdo et al. 2021, 2023) to fit the ¹²CO (2–1) line emission and kinematics channel by channel by a Keplarian disc model. DISCMINER first builds a Keplerian disc model from a set of parameters described hereafter. The position offset (x_c, y_c), the inclination i and the position angle PA define the orientation of the disc. In this work, the PA was defined as the angle from the northern axis to the blueshifted semi-major axis. The velocity profile is set by the central stellar mass M_* and the systemic velocity υ_sys. The disc upper and lower surfaces are both defined by exponentially tapered power laws (Eq. (1)). The line profile width at half maximum L_w and the line profile slope L_s are defined as power laws as well (Eq. (2), valid for both L_s and L_w). The reference intensity is taken from a peak intensity power law I_p in a region 𝓓(R_out) that extends from r = 0 to r = R_out, where R_out is the outer radius of the disc (Eq. (3)). This whole set of parameters builds the final Keplerian model of intensity I_m using a Bell function kernel with a projected velocity υ_k,los along the line of sight corresponding to the channel velocity υ_ch (Eq. (4)). In the equation below, r₀ is the reference radius taken as r₀ = 100 au, while z₀, L₀ and I₀ are values of reference for the vertical height, linewidth/lineslope and intensity respectively. The critical radius r_b describes the exponential cutoff of the disc surface, with q describing the strength of this tapering. Finally, p is the power-law index of the disc surface. $$ɀ(r)={ɀ}_0{\left(\frac{r}{r_0}\right)}^p\exp {\left(\frac{-r}{r_b}\right)}^q,$$(1) $$L(r,ɀ)=L_0{\left(\frac{r}{r_0}\right)}^{p_L}{\left(\frac{\left|ɀ\right|}{{ɀ}_0}\right)}^{q_L},$$(2) $$I_p(r,ɀ)=I_0{\left(\frac{r}{r_0}\right)}^{p_I}{\left(\frac{\left|ɀ\right|}{{ɀ}_0}\right)}^{q_I}\mathcal{D}\left(R_{\text{out }}\right),$$(3) $$I_m\left(r,ɀ,{\upsilon }_{\text{ch}}\right)=I_p{\left(1+{\left|\frac{{\upsilon }_{\text{ch}}-{\upsilon }_{k,\mathrm{los}}}{L_w}\right|}^{2L_s}\right)}^{-1},$$(4)

DISCMINER then fits the previous model to the data in the image plane channel by channel thanks to MCMC sampling using EMCEE (Foreman-Mackey et al. 2013). More details about the initial conditions and the results of that procedure can be found in Appendix C.

From the best-fit model channels, we then built moment maps from the data and model channels. First, a moment 0 map was computed by integrating the intensity in each pixel along the velocity axis. We then collapsed the cube by selecting the brightest pixel along the velocity axis for each pixel to create a peak temperature moment map (moment 8). Picking up the centroid velocity corresponding to that peak allowed for a velocity map of the disc to be built (moment 9). Finally, we built a moment map corresponding to the line-width around that cen-troid velocity. We applied a 3σ clipping in the building of data moment maps. The same mask delimiting the detected disc area was applied to the model moment maps. Residual maps were built by subtracting the masked model moment maps to the data moment maps.

3.2 Hydrodynamical simulations

Several hints of interaction with the external star V892 Tau NE are detected in the V892 Tau CBD (tidal truncation, tentative spiral arms, disc tilt, tentative disc warp). In order to confirm that V892 Tau NE is bound to the V892 system and to constrain its orbit, we performed 3D hydrodynamical simulations using the PHANTOM (Price et al. 2018b) smoothed particle hydrodynamics (SPH, e.g. Monaghan 1992) code. Our model consisted in a central binary star surrounded by a CBD and an outer stellar companion orbiting further out. The model used of 10⁶ SPH particles distributed in the disc according to a surface density distribution Σ ∝ r^−p with p = 1 that extends initially from 14 au to 200 au for a total disc mass of 0.06 M_⊙ (L21). The true orientation of the V892 Tau CBD is unconstrained due to an uncertainty on the inclination sign of the disc. We assumed that inclination to be positive and of 54.6°, meaning the closest side to the observer was the SE side. The PA of the disc was set to 53° in agreement with the observations of L21. The inner border of the simulated CBD was set to 14.2 au, which is twice the semi-major axis of the inner binary and in line with theoretical predictions of the cavity size (Miranda et al. 2017). The aspect ratio H/r = 0.055 at the reference radius R₀ = 100 au and the sound speed exponent q = 0.185 were set to match the fitted temperature profile in L21.

We modelled the stars as sink particles (Bate et al. 1995) with accretion radii of 1 au for the two components of the inner binary and of 10 au for the outer star. The inner binary was initialised coplanar in a prograde configuration with respect to the disc while its other orbital parameters and its mass were set in line with the observations (L21), which allows us to disentangle more easily the effects triggered by the companion solely. Moreover the observed inclination of the disc to the inner binary plane is measured to be down to 4°, which is close to coplanarity (L21). The orbit of the companion star V892 Tau NE remains unconstrained so far. One of the goals of this work is to explore if this companion is bound to the inner binary and to constrain its eccentricity and inclination parameters e_out and i_out respectively. From there we modelled the V892 Tau system as a triple system in five numerical setups, with a different orbit for the outer companion each time:

• ref (reference case): e_out = 0.2 & Δi_out = 0°,

• e05 (eccentric case): e_out = 0.5 & Δi_out = 0°,

• i30 (inclined case 1): e_out = 0.2 & Δi_out = 30°,

• i60 (inclined case 2): e_out = 0.2 & Δi_out = 60°,

• ei60 (eccentric inclined case): e_out = 0.5 & Δi_out = 60°,

where Δi_out represents the mutual inclination between the companion and the inner binary orbital plane. It resulted in inclination parameters with respect to the sky plane of i_out = 54.6° for ref and e05, i_out = 24.6° for i30 and i_out = −5.4° for i60 and ei60. The maximal eccentricity of the outer binary has been constrained to e_out < 0.2 considering a coplanar orbit and a disc extent of R_out ~ 200 au (L21). This estimation assumed a coplanar companion, but a misaligned orbit could allow for more eccentric orbits in line with the truncation of the disc. We explored values up to e_out = 0.5. We expect different eccentricity and inclination parameters to have consequences on the disc truncation and on the disc orientation. From dynamical considerations, it is more likely that V892 Tau NE is observed close to apoastron (van Albada 1968). We accordingly set the argument of the periapsis to ω_out = 180°. Given the projected separation of the companion of ~ 520 au (Gaia Collaboration 2023), we assumed the semi-major axis to be a_out = 500 au. The longitude of the ascending node Ω_out = 53° was set in line with the observed disc PA. We assumed the mass of V892 Tau NE to be 0.5 M_⊙, which is in agreement with the mass of similar M3 stars (Luhman 2007; Flores et al. 2022).

Table 2 summarises the parameters defining the orbital configuration of the setups. In that table, the reference plane is the plane of the sky and the origin is the centre of mass of the inner binary. Figure 3 shows the initialised orbits for each setup projected in the sky-plane and a plane perpendicular to it. We evolved the hydrodynamical simulations for 50 P_out, with P_out = 4375 yr the period of the outer binary for a semi-major axis of a_out = 500 au. It corresponds to 28 400 periods of the inner binary. At such timescales, the disc had time to relax from its initial condition and to evolve significantly.

Fig. 3 Projection of the initialised orbits of each tested setup in the sky-plane xy and in a plane perpendicular to it, namely xz. The coordinates are centred on the primary. The CBD is represented by the shaded grey area. The projected position of V892 Tau NE in the skyplane is represented by a star and by a black line in the xz plane. The orbit of the inner binary (IB) is plotted in blue and the other colours correspond to the tested orbits for V892 Tau NE.

3.3 Radiative transfer and post-processing

The simulation outputs were post-processed via radiative transfer to produce intensity datacubes of the system. Because the disc oscillates and precesses during the simulations, its inclination and position angle are evolving with time. In order to have a proper comparison with the observations, the disc was manually moved back to an average inclination of 54.6° and an average PA of 53° before the post-processing. We report the inclination and PA values of the simulated discs before this rotation in Appendix D. We used the radiative transfer code MCFOST (Pinte et al. 2006, 2009) that employs a Voronoi tesselation to build a grid based on the SPH particle distribution, with one cell per SPH particle. Since our SPH simulations contains gas particles only, the dust spatial distribution is assumed to follow the gas distribution. Dust grains were assumed to be at local thermal equilibrium and subject to passive heating. They were assumed compact and following opacity laws of astrosilicates (Weingartner & Draine 2001) with their scattering properties computed following the Mie theory framework. We used 100 grain sizes from 0.03 µm to 1000 µm logarithmically spaced in scale. The global distribution was normalised by integrating over all grain sizes, assuming a typical power-law of exponent −3.5 (Mathis et al. 1977), and all over the grid to have a dust to gas mass ratio of 0.01. We adopted a uniform value of 1 × 10⁻⁴ for the ¹²CO abundance compared to H₂. The turbulent velocity was set to zero and the freeze-out temperature below which the grid cells do not emit, to 20 K. The sink particles were considered as spherical stars that radiate isotropically as black-bodies with temperature and luminosity chosen in agreement with their mass at an age of 3 Myr (Siess et al. 2000) consistent with an estimated age of the system of ~2 Myr (Küçük & Akkaya 2010). This way each 3 M_⊙ component of the inner binary had a temperature of 10 745 K and a luminosity of 72 L_⊙ while the 0.5 M_⊙ companion had a surface temperature of 3758 K and a luminosity of 0.3 L_⊙. 1.28 × 10⁸ photons packets were used to compute the 3D temperature using a Monte-Carlo approach and 1.28 × 10⁶ other photons packets were used to compute images with a ray-tracing method based on the temperature structure.

¹²CO (2–1) channel maps were produced with a 0.075 km s⁻¹ resolution from −15 km s⁻¹ to 15 km s⁻¹. The velocity axis was then shifted from the systemic velocity υ_sys = 7.86 km s⁻¹. The produced datacube was then post-processed with CASA (McMullin et al. 2007) to produce synthetic observations. Using the simalma task, Measurement Sets (MS) were generated from the radiatively processed data with the right sampling of the uv plane. The sampling in the uv plane was taken from the antenna configurations used in the observations. Considering that the observations are made from 3 observational datasets, we produced synthetic data corresponding to each individual dataset. The integration times were set to 9.25 min for the synthetic SB and LB sets, and to 8 min for the synthetic dataset matching the OLD set. The precipitable water vapour index PWV = 0.6 mm was also chosen in line with the observing conditions. It resulted in 3 synthetic MS each matching an actual MS. The simulated MS were then spectrally Hanning-smoothed to match the spectral resolution of the observations. The 3 synthetic datasets had their continuum subtracted before being imaged together with the exact same tclean parameters used in the imaging process of the observations (see Sect. 2.1). Moment 0, moment 8 (peak temperature) and moment 9 (peak velocity) maps were created from the treated synthetic cubes using the same clipping as in the observations. The central channels where the data signal is absorbed were ignored in the process.

3.4 N-body simulations

Additionally, we performed N-body simulations using the code REBOUND (Rein & Liu 2012). The setup reproduced the V892 Tau system with 3 stars and a CBD. Stars were modelled as to massive point particles and the disc was modelled using 10 test particles with the same inclination and position angle as the disc in the SPH setups. The test particles semi-major axis extended from 14.2 au to 200 au, corresponding to the borders of the disc. Their eccentricity and mass were set to zero. In agreement with the data, the test particles start the simulation with a mutual inclination of 8° with the inner binary. Five simulations were performed with the same parameters as the SPH setups ref, e05, i30, i60, and ei60, respectively. The system is then integrated with an ias 15 integrator (Rein & Spiegel 2015) during 2000 periods of the outer binary, corresponding to ~9 Myr.

Fig. 4 Channel maps of the 12CO (2–1) emission. The DISCMINER best-fit model is plotted in colours and the data is represented by the white contour lines at [3,7]σ. Negative emission at the same significance is indicated by dashed lines. The synthesised beam is plotted at the bottom left of each channel map, while the corresponding velocity is plotted at the top left. On the bottom left channel, an arrow indicates the direction to the companion star V892 Tau NE.

4 Results

In the following, we investigate the structure of the gaseous disc and its deviations from our Keplerian model (see Sect. 3.1). We focus our analysis on these deviations and to what extent the simulations are able to capture those non-Keplerian patterns.

4.1 Channel maps

Figure 4 shows the best-fit model channels from DISCMINER and their comparison to the data. The modelled emission in the channels fits well to the bulk emission of the data and the geometry of the disc is well captured by the Keplerian model. The best-fit inclination i = 54.7°, position angle PA= 53.6° and dynamical mass M_* = 5.99 M_⊙ are in line with previous observations (L21). The systemic velocity is found to be around the channels that are experiencing absorption from the cloud at υ_sys = 7.86 km s⁻¹. The disc is thin and sharply truncated at the outer edge given the best-fit lower and upper emission surfaces. The full details of the best-fit parameters and the priors used can be found in Table C.1.

Looking at the emission in individual velocity channels allows the deviations to be tracked from the model at specific velocities. From 3 km s⁻¹ to 4 km s⁻¹, the data emission extends further than the model emission. This excess in emission is directed towards V892 Tau NE. The 3σ emission contour shows irregular features at the top of the channel emission up to 5.5 km s⁻¹. From 5 km s⁻¹ to 11.5 km s⁻¹, the emission of the lower wing is truncated with respect to its model counterpart while the emission extent of the upper wing of each channel is well reproduced by the model. Small deviations are also observed in the south-west part of the 12.5–14 km s⁻¹ channels and irregular patterns are seen especially in the 12.5 km s⁻¹ channel. These deviations could be part of a larger non-detected structure (L21), as for instance large spiral arms commonly launched by external companion at the edge of discs. We discuss that possibility in more details in Sect. 4.4. Other possible explanations for these non-Keplerian features include a warped disc or a discrepancy in vertical height between the data and the model (Law et al. 2023).

Data channels at 3.5 km s⁻¹, 11 km s⁻¹ and 12.5 km s⁻¹ show clear deviations from our Keplerian model. In Fig. 5, we look at these particular channels in the simulations and check if the deviations are reproduced or not. The lower surface of the disc is also visible in the simulations reƒ, i30, and i60 at 3.5 km s⁻¹ and 12.5 km s⁻¹, which is not the case for the other simulations e05 and ei60 and the observations. The emission in the 11 km s⁻¹ synthetic channel extends further out in reƒ, i30, i60, and ei60. In general the signal is higher in the synthetic channels than in the data. The closest match is made with the synthetic channels built from e05: the emission pattern is comparable to the data in terms of size, morphology, and intensity. The e05 disc is sharply truncated and very thin due to most of the initial disc material being ejected by the companion at the beginning of the simulation. Some deviations are also recovered in the outer disc, suggesting that the kinematics observed in the data and e05 are similar.

4.2 Disc morphology

The moment maps built from the best-fit modelled channels by DISCMINER are shown in Fig. 6, which includes the data moment maps and the residual maps in the detected disc area following the subtraction of that data and best-fit. In this section we discuss the disc morphology based on the intensity and temperature moment maps that we compare to their corresponding synthetic maps. The main features of the residual maps are highlighted in Fig. 7. The velocity and linewidth moment maps will be discussed in Sect. 4.4 and Appendix E.

On the moment 0 map, the residual patterns show two main arcs at the disc edges: one to the north of the blueshifted side and one other to the south of the redshifted side. The first arc represents a brightness excess while the second indicates a brightness deficit. The two arcs are quite symmetric with respect to the inner binary. Such residuals have already been identified in L21. The possible explanation invoked in that study was the presence of a warp in the disc and a different PA between the inner and outer disc. The arcs are also recovered in the peak temperature residual map. The latter map also shows a temperature excess oriented towards the south-west at the centre of the redshifted side of the disc (see Fig. 7) which had also been detected in L21. Thanks to a better angular resolution, we are able to see this asymmetry heading towards the western direction as it reaches the outer disc. Globally, the V892 Tau CBD is found to be a warmer disc than most T Tauri stars (e.g. Wölfer et al. 2023), heated by the two A-type stars at its heart and with a peak brightness temperatures found up to a few 110 K.

The centre of the disc is not well reproduced by our Keple-rian model according to the high residual levels in the middle of the residual maps. The modelling of the inner system of V892 Tau is complicated by the absorption of channels around the systemic velocity at the centre of the system which puts it beyond the scope of this work.

Each tested orbit of the companion leads to a different morphology of the CBD, as shown in the snapshots of Fig. 8. While the discs in the ref and e05 simulations are oriented close to their initial orientation, the inclination and the PA of the discs in i30, i60, and ei60 have significantly changed. We comment on the disc orientation in Sect. 4.5 below. The disc in e05 shows a limited radial extent due the repeated passages of the companion star close to the disc, which also spread out material that forms streamers between the companion and the CBD in all the simulations. The external companion also triggers two large spiral arms at the edge of the disc while the inner binary carves a central cavity in every simulated discs. Lastly, a thin circumstellar disc is sometimes captured by the companion at the periastron but dissipates shortly after due to the large accretion radius used for the sink particle representing the companion. One can see that the i30, i60, and ei60 discs are mildly eccentric at the end of the simulation.

Figures 9, 10 and 11 show the moment 0, moment 8 (peak temperature) and moment 9 (velocity of the brightest pixel along the spectral axis) maps respectively of the post-processed SPH simulation outputs as described in Sect. 3.3. The synthetic moment 0 maps have generally a higher flux than the data by approximately a factor 3 but have similar levels of peak intensity except for e05 of which the flux is comparable. This is mainly due to the more extended emission in the simulations than in the observations. The CBDs of ref and e05 exhibit a morphology comparable to the observations on the moment 0 map. The moment 0 maps of i30 and i60 show visible spirals in the outer disc and material forming streamers to the companion. Since such extended structures are not recovered in the data, it indicates that the companion triggers mild perturbations in the disc or an unbound scenario. The irregularities in the outer disc (see Fig. 7) are recovered on the synthetic images of every simulation and trace faint material spread out by the companion. The synthetic moment 8 maps are generally brighter than the data (by a factor 3 in brightness flux and by 50% in peak value), but the average temperature in the outer disc is well reproduced. The NW side of the disc is seen warmer (by 10% in average) and the warm twisted pattern in the redshifted side of the disc is recovered in all the simulations. The synthetic maps are more structured than the data: ref, i30, i60, and ei60 present visible perturbations towards the companion that correspond to streamers linking the disc and the companion. The e05 simulation also shows subtle perturbations at the same locations, which better matches the data. As seen in e05, material could stand beyond the outer disc and not being bright enough to be detected (see bottom second snapshot from the right on Fig. 8 in comparison to the top right image on Fig. 9). In the same fashion, it is not excluded that the V892 Tau disc exhibits non-detected large scale structures which could be probed with deeper observations. Moreover the irregularities in the outer disc are similar between the data and e05 (see Figs. 9 and 10) and could trace similar structures.

Fig. 5 Comparison of channel maps at 3.5 km s−1 (left), 11 km s−1 (middle), and 12.5 km s−1 (right) of the observed emission (top row) and of the synthetic emission (following rows, in order: ref, e05, i30, i60, and ei60). The contours represent the [3,7]σ emission levels. The beam size is represented by the grey ellipses.

4.3 Disc extent

The disc radial extent is measured in the ¹²CO observational data up to R_90% = 1.45″ that corresponds to 194 au with the adopted distance of d = 134.5 pc. Figure 12 shows the intensity profiles of the observations and of the simulated discs and the resulting R_90% for these discs. The intensity profile drops inside the dust ring but the imaging of a potential inner cavity is made difficult by the absorption around the systemic velocity. The disc size and the radial intensity profile measured in our observations are comparable to those in previous works (see Table 1 and L21).

In our hydrodynamical simulations, the inner binary quickly carves a cavity in the inner gas disc. This cavity is 4–5a_in large, where α_in is the semi-major axis of the inner binary. Given the eccentricity and mass ratio of the V892 Tau inner binary, the inner rim of the cavity should lie at ~3.5 a_in in the case of a coplanar disc (Hirsh et al. 2020) which is in rough agreement with our simulations. The discrepancy could be explained by the large accretion radii used for the sink particles of the inner binary components. Due to the long duration of our simulations and the high computational cost of small accretion radii, modelling the inner disc precisely can not be achieved precisely with our numerical setup and we rather focus on the truncation of the disc by the external companion star.

Tidal truncation by the close outer companion swiftly takes place early in the simulations and material from the disc is spread out in the outer system. After 50 simulated orbits of the outer companion, the spread material has been cleared in ref, e05, and i30. In i60 and ei60, the size of the bulk of the disc is set after a few orbits of the companion even if particles are still ejected from the disc at the end of the simulation. If the companion follows an eccentric orbit, its periastron is closer to the disc and the truncation is stronger, resulting in a more compact disc. When the eccentricity is lower, the disc is allowed to spread radially. In the simulations ref, i30, and i60, the disc is much larger due to the truncation not being efficient enough or due to material being spread out by the companion. In the ei60 simulation, the eccentricity of the companion is higher but as is the inclination which result in a larger closest distance approach. Added to the material populating the outer system and contributing significantly to the disc emission, the estimations of the disc radius are higher in ei60. In the end, the observations are in better agreement with the e05 simulation as both show comparable disc sizes.

Fig. 6 12CO (2–1) moment maps of V892 Tau. The integrated intensity, the peak temperature, the peak velocity, and the line width correspond to the first, second, third, and last rows, respectively. The first column shows moment maps computed from the data emission channels, the second column from the model and the last column shows the residuals in the detected disc area computed from the two previous maps. A mask defined by the 3σ contour of the data has been applied to the model before comparison with the data.

Fig. 7 Sketch of the features identified in the CBD of V892 Tau from the residual maps shown in Fig. 6. The black line corresponds to the 3σ contour of the disc emission. The indicated high brightness residuals come from the moment 0 residuals. The spiral arms are identified on the velocity residual map, from which only the red part is plotted and the centre masked for clarity reasons. The temperature asymmetry coming from the peak temperature residual map is overlaid in dark brown. The arrow in the top left indicates the direction of V892 Tau NE.

4.4 Spiral arms and kinematics

The launch of spiral arms in a disc by an external companion is a well established result (Rafikov 2002) and has been observed in various systems already (e.g. HD 100453, Benisty et al. 2017; UX Tau, Ménard et al. 2020). A binary star at the heart of a CBD can produce spiral features as well (Poblete et al. 2019), that are most visible at the cavity edge that propagate towards the outer disc (e.g. AB Aur, Poblete et al. 2020; HD 142527; Price et al. 2018a) whereas an external companion triggers two wide open spiral arms at the Lindblad resonances, where one arm is pointing towards the companion and the other lays in the symmetrical side of the disc. For an orbiting companion, the scale of the spirals is reduced as the disc is truncated with time (see Fig. 8 and simulations in Ménard et al. 2020).

The velocity residual map on Fig. 6 shows the disc as dominated by the Keplerian rotation. The central regions of the system exhibit pronounced residuals, which may stem from cloud absorption within the corresponding channel emission. The residual patterns suggest that two opposite spiral arms are traversing the disc. One arm delineated by red residuals extends from the inner blueshifted region to the outer redshifted side, passing through the northern section of the disc. Symmetrically the other arm is marked by blue residuals in the southern regions of the disc. We highlight these spiral patterns on the sketch shown in Fig. 7. Spirals have been tentatively identified in the outer redshifted part of the disc (L21): a region where the disc is irregular and where brightness residuals are high. Our observations support the presence of a spiral arm at that location and unveil its potential symmetrical component. In the light of our simulations that we describe below, we present this given geometry of trailing spirals in the disc as a tentative evidence of CBD-companion interactions. The presence of these interactions is also supported by the rotation curve derived from the gas emission model, and that we detail in Appendix C.2. Due to the scale and prominence of the spirals, the best match is done with the e05 setup. Non-Keplerian features in the 9–11 .5 km s⁻¹ channel maps could also trace a large scale spiral structure and the truncated lower wing of the channel map emission mentioned in Sect. 4.1 could be consistent with a spiral perturbation at that location.

In every simulations performed, the outer companion triggers two large spiral arms in the disc when passing by the periastron. The spirals are sustained throughout the simulation runs and are visible in the column density rendered plots of the SPH setups of Fig. 8. After the post-processing done, hints of spiral structures remain visible in the synthetic moment maps. In ref, i30, i60 and ei60 the two spiral arms are traced by the gas emission on the moment 0 map where patterns are mainly visible to the south and north-west of the CBD. The peak velocity maps also show bright emission at the same locations in these simulations. However no clear spiral structure is seen in the moment 0 map of the e05 simulation and no clear deviation is observed in the corresponding peak velocity map of that simulation. In the small disc of the e05 simulation, the spirals disperse more rapidly and are of lower strength making their detection difficult when the companion comes close to apoastron. It suggests that faint spiral arms could be present in the V892 Tau disc.

Gravitational instability and embedded planets are also sources of spiral patterns in discs. However the triggering of the gravitational instability requires a disc massive enough for its self-gravity to dominate the gas pressure and rotation (Goldreich & Lynden-Bell 1965). This is expected when M_disc > 0.1 M_* (Kratter & Lodato 2016), which is not the case of the V892 Tau CBD where M_disc ~ 0.01 M_*. Embedded planets interactions with their host disc can result in the formation of two spiral arms with one of them pointing to the planet (e.g. Dong et al. 2015). It has also been showed that internal massive planets could reproduce spiral features in discs (Calcino et al. 2020). In our N-body simulations, planets in the inner system are ejected of the system by the inner binary (see Sect. 4.7) making the scenario of a planet inside the dust cavity unlikely. Planets could survive at larger separations and trigger the observed spiral patterns, but our data does not allow us to comment on that possibility.

Fig. 8 Rendered plots of the density spatial distribution integrated along the line of sight of the different simulated setups. The reƒ, e05, i30, i60, and ei60 setups are shown from left to right. The top row shows the system after one orbit of the outer companion. The bottom row displays the state of the system after approximately 50 orbits of the outer companion, corresponding to 0.2 Myr and 28 400 orbits of the inner binary. The white dots indicate the sink particles representing the stars.

Fig. 9 12CO (2–1) integrated intensity map of the observational data (top left) compared to the synthetic maps built from the simulations reƒ (top middle), e05 (top right), i30 (bottom left), i60 (bottom middle), and ei60 (bottom right). The discs were rotated back in the observed disc plane to allow for a proper comparison to the observations. The beam sizes are plotted in the bottom left corners.

Fig. 10 Same as Fig. 9 but for the 12CO (2–1) peak temperature map.

Fig. 11 Same as Fig. 9 but for the 12CO (2–1) emission velocity map.

Fig. 12 Azimuthally averaged radial profile of the 12CO intensity for the observations and the simulations. The vertical dashed lines indicate the R90% radius for the continuum, the gas emission, and the synthetic observations. The shaded grey area corresponds to the major axis of the synthesised beam.

4.5 Inclination and orientation of the disc

The tilt Δ is defined as the angle between the disc and the inner binary orbital plane. From our fitting process, we recovered a disc orientation (i, PA) very similar to what was derived in previous works (L21). Thus, we assumed the disc to have the same tilt as previously calculated. Due to an ambiguity of 180° on the inner binary position angle, the two possible tilt values are Δ = 8.0 ± 4.2° and Δ = 113.2 ± 3.0°. We believe the disc to be close to coplanarity and rotating prograde with the companion star. In this way, we assumed a tilt of Δ ~ 8.0° to be more likely. We further motivate this assumption in Sect. 5.1. We measure the tilt Δ and the position angle (PA or Ω) of the simulated discs by averaging the tilt and PA of all the SPH particles tilt and PA having a semi-major axis between the inner border of the disc and 200 au, which corresponds to an upper limit of the disc extent.

On the one hand, in the setups where the companion is on a coplanar orbit with the disc and the inner binary (ref and e05) the inclination and the PA of the disc remained the same as the initialised values during the entire simulation. The tilt Δ and PA are then of 0° and 53° respectively at the end of the simulation for these setups. On the other hand, in the setups where the companion is on an inclined orbit with respect to the disc and the inner binary, namely i30, i60, and ei60, both the inclination and the PA of the disc oscillate periodically with time, as seen in Fig. 13. The tilt oscillations damp with time and tend to a nonzero value with an amplitude of the oscillations starting at about 14°, 17° and 18° for i30, i60, and ei60 respectively. The PA of the disc decreases in average while oscillating with an amplitude of the oscillations appearing to be maintained.

For the inclined setups i30, i60, and ei60, we fitted the inclination of the disc with a damped oscillations model (see Eq. (5)). The PA is fitted by a linear model on the top of which are sinusoidal oscillations (see Eq. (6)). The model is adjusted to the simulation data using a non-linear least squares method. The disc relaxes from the initial conditions during the first orbits of the outer binary. We exclude at least the first 15 orbits of the outer companion of the fitting process to neglect that effect. The best-fit results can be found in Table 3. $$\Delta (t)=a_1+a_2{e}^{-a_{3}t}\cos \left(2\pi a_{4}t+a_5\right),$$(5) $$\Omega (t)=b_1+b_{2}t+b_3\cos \left(2\pi b_{4}t+b_5\right).$$(6)

The best-fit asymptotic tilt value is around 6° and 8° for i30 and i60 respectively. We notice that the i60 value is in excellent agreement with the value derived in L21. Yet the i30 value is consistent as well within the errorbars. In the case of ei60, the model indicates an asymptotic tilt value of 3.5°. The b₂ parameter of the PA model can be converted into a precession timescale of 1619 P_out, 923 P_out and 1003 P_out for i30, i60, and ei60, respectively. The uncertainty of the fit for these values is around 100 P_out. The a₃ parameter of the tilt model can be converted into a damping time that is 85 P_out, 226 P_out and 19 P_out, which corresponds to approximately 0.3 Myr, 0.9 Myr and 0.08 Myr for i30, i60, and ei60, respectively. The incertitude on these values is <20 P_out. The a₄ and b₄ parameters can be converted into an oscillation period of the disc that is fitted on the tilt profile to 14.8 P_out for i30, 19.3 P_out for i60 and 10.1 P_out for ei60. We applied the same process to the PA oscillations to find a period of 15.2 P_out for i30, 20.9 P_out for i60 and 10.2 P_out for ei60. The uncertainty of the fit for these values is of < 0.2 P_out. The oscillation period is shorter in ei60 than in i30 and i60. We note that the oscillation period of the tilt and the PA are similar for a given simulation. We also note that our model does not account well for the ei60 disc dynamics where the oscillations are quickly damped and the disc ends up close to coplanarity with the inner binary. We discuss in more detail those dynamical behaviours in Sects. 5.1 and 5.2.

The previous results indicate that an inclined companion with respect to the inner binary is able to misalign a disc initially in the inner binary plane. It proves that a misaligned geometry between the inner binary and the outer companion is needed to explain the observed tilt value of Δ_obs = 8° derived in L21, which is best reproduced by an inclination of 60° of the companion and a low eccentricity of it (i60). Considering the oscillations of the disc tilt, the disc could be caught in a middle of an oscillation and be measured misaligned with the inner binary. However when comparing the damping time of these oscillations to the estimated age of the system of ~2 Myr (Küçük & Akkaya 2010), it is unlikely that the disc is still undergoing significant oscillations.

Fig. 13 Tilt (top) and PA (bottom) evolution of the disc for the i30 (left), i60 (middle), and ei60 (right) setups. The data are represented by the coloured points while the fit is shown by the red curve. The tilt evolution was modelled by damped oscillations. The PA evolution was modelled by oscillations on the top of a linear precession. The best-fit parameters can be found in Table 3.

4.6 Disc warp

A warp is a non-constant inclination radial profile appearing preferentially in thinner, less viscous and larger discs (Papaloizou & Pringle 1983; Young et al. 2023). Disc warping is known to occur in CBDs around eccentric and inclined binaries (e.g. Nixon et al. 2013; Lodato & Facchini 2013) and in circumstellar disc perturbed by an external misaligned companion (Gonzalez et al. 2020). In non-coplanar triple systems, these two effects occurs simultaneously. Given the previous results, a misaligned orbit of V892 Tau NE could explain the observed disc tilt (see Sect. 4.5). If such an orbit is able to misalign the disc, it could be able in principle to trigger a warp in the disc too. In the following we investigate the presence of a warp in the CBD of V892 Tau.

In order to search for the presence of a warp in our simulated discs, the disc was binned in semi-major axis from the inner cavity edge r_in to r = 200 au and the tilt Δ was measured in each bin taking the standard deviation of the tilt of the particles in a bin as the corresponding error. In the SPH data, the cavity semi-major axis a_cav was defined as ρ(a_cav) = ρ_max/2 with ρ_max defined as the maximum of the radial density profile (Artymowicz & Lubow 1994). Since the disc eccentricity stays below 0.1 during the simulation, the cavity semi-major axis can be identified as the cavity radius. At t = 50 P_out, a_cav = 45 au for i30 and a_cav = 44 au for i60.

We inspect and plot the radial distribution of the tilt Δ(r) of the setups i30, i60 and ei60 at different time steps in Fig. 14: t = 10 P_out, t = 30 P_out and t = 50 P_out. The inner binary, the outer binary and the discs stayed in a coplanar configuration for the initially coplanar setups reƒand e05. In the setups where the companion is on an inclined orbit (namely i30, i60, and ei60) a warp forms steadily during the first orbits. The tilt profile Δ(r) shows a continuously but modestly increasing trend: the inner parts of the disc are misaligned by 1° to 4° with respect to the outermost parts after a few orbits of the companion. That tilt profile is maintained during the early moments of the simulation in all the simulations with an inclined companion. However the disc becomes planar again within the errorbars after a time t = 30 P_out. The companion bends the outer parts of the CBD at the beginning of the simulation, which creates a warp propagating in the bending wave regime (Papaloizou & Terquem 1995) to the innermost parts of the disc. Because of the loss of material in the disc and its truncation by the companion, the disc behaves more and more rigidly with time. In this way, the strongest warping occurs at the beginning of the simulation. After a few orbits of the companion once the disc truncated enough, the tidal torque of the companion is insufficient to warp the disc: the warp is damped and the disc becomes planar again Deng et al. (2021). At the end of the simulations, the disc can be considered planar within the errorbars and the amplitude of a potential warp in the disc would be only that of a few degrees. Measuring such a low warp our synthetic data is challenging and require a more in-depth analysis that we left for future studies. In the meantime, we favour a planar configuration for the discs in our simulations.

A warp in the CBD has been invoked to explain deviations from simple models (Long et al. 2021; Vides et al. 2023). Residuals from our Keplerian model are similar to these deviations: the arcs at the disc edges on the moment 0 and peak temperature residual maps (see Fig. 6 and Sect. 4.2) could be hints of a misalignment between the inner disc and the outer disc. According to the orbital parameters of the inner binary, the disc could be warped under the tidal torque of the binary if sufficiently misaligned (see Sect. 5.1). V892 Tau NE could also warp the disc, if on a misaligned orbit with respect to the disc. We found the CBD of V892 Tau to be a thin disc with an emission surface aspect ratio of z/r ~ 0.16 at 100 au, where most discs around young stars are found with z/r > 0.3 in the ¹²CO emission (Law et al. 2021). The disc of V892 Tau is also moderately extended. Thus we expect warping to be difficult in the disc. The high temperature in the disc up to ~110 K makes it prone to damp any potential warp (e.g. Rabago et al. 2023). In the end, even if the observed residuals could be consistent with a warped disc, the combined dataset used in this work can not robustly rule out nor confirm the presence of a warp in the disc.

Fig. 14 Tilt radial profiles of various simulations at t = 10 Pout (left), t = 30 Pout (middle), and t = 50 Pout (right). The orange, blue, and teal dots correspond to i30, i60, and ei60, respectively. The 1σ errors are indicated by the shaded areas.

4.7 The orbit of V892 Tau NE

4.7.1 The flyby scenario

The orbit of V892 Tau NE is currently unconstrained. This means that V892 Tau NE could in principle be totally unbound to V892 Tau or following a flyby parabolic/hyperbolic orbit. In the unbound case, V892 Tau would be considered as a binary system. In such a configuration, a coplanar disc with respect to the inner binary and does not precess (see Appendix F for details). It means that the inner binary alone cannot explain the observed disc tilt, which allows us to rule out the unbound scenario. Moreover it has been recently shown that V892 Tau and V892 Tau NE are likely related (Thomas et al. 2023). In the case of a flyby event, efficient truncation and the formation of prominent short-lived spirals are expected (Cuello et al. 2019, 2020, 2023; Smallwood et al. 2023). We can study the flyby scenario in our simulations by considering the state of the disc after the first orbit of the companion. Figure 8 shows column density views of the simulated discs after the completion of the first companion orbit. In all the simulations, dense spiral arms are triggered in the disc by the passage of the companion. The disc remains close to its original size, that is similar to the observed one, in all the simulations despite material being ejected. If the companion was to follow an inclined orbit with respect to the inner binary, Fig. 13 shows that the disc would be misaligned by 1–2 deg at the end of the first orbit considered here as the end of the flyby event. In our simulations, the companion approaches the disc closely with a closest distance no greater than 200 au. A more distant flyby would have produced fainter spirals and lower misalignments (Cuello et al. 2019). In any case, the subtle observed spiral patterns and the disc tilt of ~8 deg suggest that a flyby event is unlikely to have happened in V892 Tau.

4.7.2 Constraints on a bound orbit

Table 4 summarises the constraints placed on the V892 Tau NE orbit through our modelling of the V892 Tau disc. If bound, a medium eccentricity of the orbit of V892 Tau NE is needed to explain the observed disc extent while the orbit has also to be misaligned with respect to the inner binary to explain the slight misalignment between the disc and the inner binary. In such a configuration, the orbit of V892 Tau NE cannot be too eccentric (e > 0.5), as high eccentricities tend to align the disc with the inner binary (Fig. 13).

From the points of the parameter space we probed with our simulations, we suggest constraints on the parameter eccentricity and mutual inclination of V892 Tau NE of 0.2 < e < 0.5 and 30° < Δi < 60° respectively. If V892 Tau NE was to follow an eccentric inclined orbit, we expect the CBD to oscillate and precess in the light of the discs in the setups i30, i60, and ei60. According to our simulations, the CBD inclination oscillates with time and these oscillations are damped over a period of ~1000 yr. Then the tilt of the disc reaches a non zero final value of 3° to 8° which is in good agreement with the observed value of ~8°. In such a configuration, we expect the disc PA to oscillate with time on top of a linear precession. We discuss the non-coplanar configuration of the system and the resulting disc dynamics in the Sects. 5.1 and 5.2 below.

4.8 Stability of future planets

Considering the results of the previous Sect. 4.7, we consider the dynamics of planets in a system following the ei60 orbital configuration. Although not fully encompassing all disc features, the parameter set represented by ei60 delineates the boundary of the parameter space governing an inclined eccentric orbit of the companion. We anticipate that this particular configuration will induce the most pronounced perturbations on planets within the system, as alternative configurations would entail lower eccentricity and/or mutual inclination of the companion. Consequently, we present the following results as a worst-case scenario for the prospective planets in the system. The results for the other N-body setups can be found in Appendix G.

The stability of circumbinary planets in a system in the ei60 configuration depends on the planets semi-major axis first. Figure 15 shows the evolution of the test particles in the N-body simulation of ei60, which is described in Sect. 3.4. The test particles with a semi-major axis a ≥ 100 au were ejected from the system during the simulation. The innermost test particle at a = 20 au was also gravitationally ejected from the system by the inner binary during the first orbit of the outer companion. Particles in between a = 40 au and a = 80 au remained stable during the simulation. This can be explained by the competition between two dynamical effects. First due to the presence of the misaligned companion, the test particles undergo von Zeipel-Kozai-Lidov (vZKL hereafter) cycles (von Zeipel 1910; Kozai 1962; Lidov 1962) with a characteristic timescale of ~90 P_out at 100 au (Ceppi et al. 2023). This effect results in a rise and oscillations in eccentricity and inclination of the test particles leading to the ejection of the test particles closer to the external star. Second the inner binary induces a nodal precession of the test particles around its angular momentum vector with a timescale of ~98 P_out at 100 au (e.g. Farago & Laskar 2010; Lodato & Facchini 2013). Closer to the inner binary, the timescale of this precession is lower than the vZKL timescale, stabilising the particles in the inner half of the disc (a ≤ 100 au) (Verrier & Evans 2009). Still, the particles found stable at the end of the simulation were undergoing large oscillations in inclination (and eccentricity for the test particle at a = 80 au) that questions their future stability on longer timescales. Thus if the current gas disc is ensuring the stability of forming planets, those planets could become unstable once the gas disc dissipated. On longer timescales, it implies the migration of circumbinary forming planets towards stable regions at an intermediate distance from the inner binary between 40 au and 80 au. We discuss these results in more details in the Sect. 5.3.

5 Discussion

5.1 Final inclination of the disc

CBDs are expected to reach a stable coplanar or polar configuration depending on the properties of the inner binary and of the initial mutual inclination. If the initial misalignment and the eccentricity of the inner binary are low, the disc is expected to move into a coplanar configuration (Facchini et al. 2013). On the opposite, for binaries with higher eccentricity and a high mutual inclination the CBD is expected to evolve towards a polar configuration (Aly et al. 2015; Martin & Lubow 2017; Cuello & Giuppone 2019). Nonetheless CBDs have been found with a medium mutual inclination with their inner binary in a configuration that is not coplanar nor polar. For example the CBD around GG Tau A is inclined of about 25° with respected to the inner binary (Köhler 2011; Andrews et al. 2014). In SR24, the tilt of the CBD is around 35° with respect to the binary SR24N (Fernández-López et al. 2017). These non-zero tilt values could be explained by oscillations around the final configuration of the disc before reaching the equilibrium coplanar configuration. In hierarchical triple systems as V892 Tau, the picture is modified by the evolution of the stellar orbits that can occur on shorter timescales with respect to the disc lifetime. In particular in misaligned systems, the inclination of the orbits can change with time meaning that the intermediate tilt of CBDs could be due to the slow reaction of the disc to these variations (Ceppi et al. 2023). This means that the disc of V892 Tau could eventually be seen misaligned with the inner binary even if the orbit of the companion has a different inclination. Such a misalignment would be created over the timescale of the orbit evolution that we assume to be the precession timescale of the eccentricity vector. We find that the precession timescale of the inner (outer) binary is about 6000 Myr (24 Myr) for V892 Tau (Ceppi et al. 2023). Compared to a typical disc lifetime of 10 Myr (Ribas et al. 2015) and to our 0.2 Myr-long simulation runs, the orbital configurations tested in ref and e05 setups are not expected to become misaligned with respect to the disc.

The CBD of V892 Tau is found with a misalignment of approximately 8° or 113° with the inner binary of eccentricity 0.25 (L21; Vides et al. 2023). The two possible solutions derive from a lack of radial velocity information on the inner binary, making its true inclination difficult to constrain. If considered as a CBD in a triple system with the ei60 orbital configuration, the disc of V892 Tau could go to a stable polar configuration (Martin et al. 2022; Ceppi et al. 2023). The large semi major axis and the low mass of the companion do not hinder the polar alignment. However, polar alignment would require the disc to be initially misaligned by 66° or more (L21). In our simulations, the disc was initialised coplanar with the inner binary. Our simulations indicates that a misalignment of ~8° is created and sustained by the outer companion. The disc oscillates around that equilibrium position. Given the proximity of our best-fit tilt value, we prefer the nearly coplanar configuration for the disc in the end. This non-zero tilt value shows that intermediate configurations are stable for CBDs, independently of secular dynamics of the stars in the system. The V892 Tau disc is in a sweet-spot of the parameter space that would allow this kind of configuration. In order to achieve this, the semi-major axis ratio between the inner and outer binaries should be relatively large and the inclination of the companion with respect to the inner binary around a medium value of ~45°. Future studies to characterise CBDs inclination in triple system as a function of the companion inclination would help to better understand this problem.

Fig. 15 Evolution of the semi-major axis (top), eccentricity (middle) and mutual inclination with the inner binary (bottom) of the test particles in the ei60 N-body setup. Each colour represents a different test particle that started at a given semi-major axis. The semi-major axis of a particle falls to 0 once ejected from the system. The particle starting with a = 20 au was ejected from the system before the completion of the first orbit of the outer binary.

5.2 Oscillations and precession

Due to gravitational torques applied by misaligned stars, discs in multiple systems very commonly undergo precession in the sense that their PA is evolving with time. Linear precession can occur in binary systems in CBD and circumstellar discs (e.g. Papaloizou & Terquem 1995). Simulations of binary systems have already observed this phenomenon (e.g. HD 100453 Gonzalez et al. 2020; Nealon et al. 2020; GW Ori Smallwood et al. 2021).

We find that the V892 Tau CBD is linearly precessing but that its PA also oscillates with time (see Fig. 13). Assuming V892 Tau and V892 Tau NE as a binary system, the expected precession timescale of the disc is about 144 P_out for i30 and 251 P_out for i60 (following Eq. (6) from Gonzalez et al. 2020 derived from Bate et al. 2000 and Terquem 1998). This is well below the precession timescales found in Sect. 4.5. However these estimations are in line with precession timescales measured in simulations where the system is represented by only a binary. Those simulations are identical to i30 and i60 except that the inner binary is replaced by a single star of the same mass (see Appendix F). As seen in Fig. F.1, the evolution of the tilt and PA of the disc in those binary-system simulations is different from the one in the triple-system simulations. The discs in the i30 and i60 runs start to linearly precess at the same rate as in the binary-system simulations, but their precession is then slowed down by the inner binary. It indicates that the disc dynamics are set by the triplicity of the system, that can not be simplified as a binary system.

In binary systems, a misaligned circumstellar disc is expected to align to the binary orbital plane in about a precession timescale in the presence of hydrodynamical instabilities such as the warp developed at the beginning of our simulations (Bate et al. 2000). In the binary-system simulations, the disc starts to align on a timescale in line with that prediction. However in i30, i60 (and ei60) the disc oscillates and ends up misaligned with the orbital planes. The gravitational torques of both the inner and the outer binary try to align the disc in their respective planes that results in an equilibrium position in between those planes. It shows again that when considering the disc dynamics the misaligned hierarchical triple system V892 Tau can not be reduced to a binary system.

In the i60 and ei60 simulation setups, vZKL oscillations could be triggered in the CBD that would result in the disc eccentricity and inclination oscillating with time (Martin et al. 2014). As the inclination oscillation pattern are similar in i30 and i60, this strength of this effect is low in i60. However knowing that the vZKL timescale decreases with a higher eccentricity of the outer orbit (e.g. Kiseleva et al. 1998), vZKL cycles could explain the dynamical behaviour. The vZKL timescale for ei60 is about ~16 P_out at 100 au (Ceppi et al. 2023), which is comparable to the oscillation period in ei60 (see Fig. 13). Moreover the CBD of the ei60 SPH simulation shows clear eccentricity, showing the clear influence of the vZKL mechanism on this disc.

The disc is subject to the combined gravitational torques of the external companion and of the inner binary. While in the case of a binary system a linear precession is expected, in a triple system oscillations are also expected. This peculiar precession can be assessed by measuring the tilt Δ and the twist β. These oscillations, projected in the 2D sky plane, are traced by the disc inclination and PA. Previous studies have used complex variables to examine the dynamical evolution of discs (e.g. for eccentricity modes in the context of planet-disc interactions, Ragusa et al. 2018). A similar study of the complex variable Δ = |Δ| exp(jβ), where j is the complex number, could help to clarify the connection between the tilt and twist oscillations and the resulting evolution of misaligned CBDs in triple systems.

5.3 Resulting planetary architectures

In its current configuration, the misaligned disc of V892 Tau could result in the formation of misaligned planets with respect to the inner binary. Though those planets would have to migrate to the inner half of the disc to survive in the system (see Sect. 4.8). It is unclear to what extent this migration would impact the orbital elements of the planet. If the planet was to tilt relatively to the disc plane, this tilt would be damped and the planet would re-align (Burns 1976; Tanaka & Ward 2004). In the case of V892 Tau the resulting planets would stay close to the disc plane, inclined with respect to inner binary. The stability of inclined planets in hierarchical triple systems depends on the stellar orbital configuration. In triple systems with low semi-major axis ratio (<20) and eccentricity of the outer orbit, the outer bound of the stability zone is limited. If the mutual inclination between the companion star and the inner binary is not zero, the planets stability region extends vertically allowing for high inclination (Busetti et al. 2018). In that configuration, circular planets that form coplanar with the inner binary end up with their i and Ω parameters oscillating on timescales several order of magnitude longer than the orbital period (Busetti et al. 2018). The inclination between the planet and the inner binary also sets the stability regime. Planets are generally found stable close to coplanarity with the inner binary, regardless of the inclination of the outer star. In the ei60 configuration of V892 Tau, the stability region extends up to ~ 100 au from the inner binary and up to misalignments of the planet of ~30° (Gianuzzi et al., in prep.). Putting that result together with the results described in Sect. 4.8, we expect future planets in the V892 Tau system to form at intermediate distance from the centre of the system between 40 au to 80 au with mild inclination with respect to the inner binary.

The dust disc of V892 Tau is one of the most massive Class II disc (L21), potentially presenting an abundant pebble reservoir. The horseshoe asymmetry seen in the continuum could be explained by an over-density created by the inner binary, which could enhance planet formation at that particular location. Misaligned and warped discs around eccentric binaries are also subject to dust concentrations triggered by the difference in precession timescale between dust and gas, which results in dust traps and favours planet formation (Aly et al. 2021, 2023). More advanced dust+gas hydrodynamical simulations would allow for the presence of warp instabilities in the V892 Tau CBD to be assessed and this would help to understand planet formation in triple systems in general.

5.4 Circumbinary discs interactions with an external companion

In the system of V892 Tau, the inner binary dominates the mass budget and the low mass companion V892 Tau NE tidally interacts with the CBD. It induces spiral patterns within the disc, a truncation of the outer disc, a slight misalignment with the inner binary but no sustained warped geometry. In the following we explore how, given observations of a CBD, similar structures can be explained by ongoing interactions with an hypothetical companion. This companion may either follow a bound orbit around the inner binary or have been involved in a flyby event. In the absence of identifiable candidates in the vicinity of the CBD, Gaia astrometry remains a viable method to search for potential past external flyby occurrences.

Repeated or not, the passage of a star close to the CBD would dramatically reduce the size of the CBD by tidal truncation (e.g. Artymowicz & Lubow 1994; Miranda & Lai 2015), providing a potential explanation of small CBDs. The efficiency of truncation increases with a higher mass ratio between the inner binary and the outer companion, as well as with greater eccentricity or a smaller semi-major axis of the companion. During a close encounter, the companion may also accrete material from the CBD, leading to the formation of a new disc and associated streamers (e.g. VLA 1623-2417 Mercimek et al. 2023, FU Ori-like systems Borchert et al. 2022).

CBD-binary misalignments are common and expected in binary systems and can find their origin in interactions with surrounding stars (Czekala et al. 2019; Elsender et al. 2023). Indeed an external perturber would exert a gravitational torque on the CBD, if misaligned with respect to each other, aligning the CBD with the orbital plane of the perturber (Papaloizou & Terquem 1995). Should the misalignment between the CBD and the inner binary reach a critical angle, the CBD may transition towards a polar configuration relatively to the (eccentric) inner binary (Martin & Lubow 2017; Martin et al. 2022; Ceppi et al. 2023). The disc can also become warped in the process or break at a given radius depending on the inner binary parameters (e.g. Nixon et al. 2013; Lodato & Facchini 2013). All these effects are enhanced by large companion to inner binary mass ratios and high misalignment degrees to the perturber, or in the case of a retrograde flyby (Cuello et al. 2019).

The observation of two diametrically opposed spirals in the CBDs, as in AS 205 (Kurtovic et al. 2018), suggests the occurrence of a flyby event. Indeed, any close encounter would invariably induce spiral density waves within the disc through tidal effects (Rafikov 2002). Prograde encounters and close fly-bys tend to induce the formation of prominent spiral arms, whereas bound companions typically result in less conspicuous spiral arms (e.g. Cuello et al. 2019; Ménard et al. 2020).

However if none of the aforementioned structures are detected in a CBD, it remains plausible that the disc has undergone relaxation over time due to prior interactions with neighbouring stars. Lastly, some multiple systems show a lack of CBD (e.g. HT Lup; Kurtovic et al. 2018, Sz65&Sz66; Miley et al. 2024). This absence could be attributed to the considerable separation between the binary components precluding the formation of a CBD (Harris et al. 2012; Elsender et al. 2023). Moreover the hypothetical existence of such a disc suggests rapid dispersal, possibly stemming from intense interactions with a neighbouring star. These interactions could range from dramatic events like disc-penetrating flybys (Cuello et al. 2023) to the influence of additional radiation enhancing photoevap-oration (Dai et al. 2018). In any case, a detailed systematic study of CBD in triple systems exploring the wide range of available orbital configurations will help to properly understand the previously cited effects.

5.5 Limitations and caveats

The main limitation of our work is the low number of orbital configurations tested. Considering the total lack of orbital information on V892 Tau NE, we proposed orbital arrangements that could trigger the observed structures in the CBD. Admittedly our goal was not to thoroughly sample the available parameter space for the orbit of V892 Tau NE, which would have been too computationally demanding with SPH simulations. Instead we chose to model V892 Tau as a triple system with a CBD in order to put constraints on the companion’s orbit. Without any astrometric measurements available and considering the well-characterised inner binary orbit, disc forward modelling constitutes the most reliable method to infer the long-period orbit of V892 Tau NE.

We ensured an optimal quality of the data by combining previous datasets with new observations. Our disc model (surface density exponent, aspect ratio, viscosity, etc.) and the values used in the radiative calculations (dust to gas mass ratio, CO abundance, freeze-out temperature, etc.) are all typical values found in the literature for Class II protoplanetary discs. We used gas only simulations for computational reasons. Therefore the dynamics of large dust grains (marginally coupled to gas) is not fully consistent. At any rate, in this work we focus on an in-depth study of the gas morphology and kinematics, which is marginally affected by our assumption regarding the gas-dust coupling.

Differences between the data and our synthetic maps could come from our numerical setup. Due to a lack of numerical resolution, the inner system is not well reproduced in our simulations. Absorption of the parent cloud at the centre of the system also limits its precise modelling. For example the circumprimary disc can not be captured due to an accretion radius of the sink particles being much larger than its expected physical size. This could explain the flux difference between our images and the data, as this narrow disc could capture incoming radiation of the primary. Also the large accretion radii of the sinks translate into an inner cavity wider than the one observed.

Last, the expense of the presented hydrodynamical simulations prohibits us from evolving the system until its equilibrium state. However our simulations reached a state close to this equilibrium, as the size of the bulk of the disc does not significantly change in the last 20 orbits of the companion. Therefore our simulations are long enough to reasonably infer the long-term disc behaviour. For longer evolutionary times, we expect the companion to clear even more material in the outer regions, which should lead to smaller disc sizes–potentially in better agreement with the observations, especially for ei60.

6 Conclusion

In this paper, we presented new ALMA Band 6 observations of V892 Tau at a high spatial and spectral resolution. We confirmed the presence of structures in the disc, which we interpreted as hints of an interaction with the companion star V892 Tau NE. Using 3D SPH simulations, we modelled V892 Tau as a triple system with a CBD and produced synthetic ALMA maps comparable to the observations. We investigated the available orbits for V892 Tau NE on the basis of their ability to reproduce the data. Once the orbital configuration was constrained, we discussed the resulting disc dynamics and planetary architectures.

Our main findings are summarised as follows:

• 1.

  Our observations confirm the hints of an ongoing interaction between the circumbinary disc around V892 Tau and V892 Tau NE. The disc extent is consistent with tidal truncation models. Even if dominated by Keplerian rotation, the CBD harbours non-Keplerian structures that are seen in the channel maps and moment maps, such as spiral features. Finally, the data neither conclusively support nor refute the existence of a warp in the disc.

• 2.

  An eccentric inclined orbit of V892 Tau NE reproduces features compatible with the disc observations. A non-zero mutual inclination of the companion with respect to the inner binary best explains the slight tilt of the disc. Given our limited parameter space sampling, this tilt seems to be reproduced best by a mutual inclination of 60°. The eccentricity of the companion sets the disc extent, which seems to be reproduced best by eccentricity values around 0.5.

• 3.

  The misalignment between the inner binary and V892 Tau NE triggers disc damped oscillations in inclination and a non-linear precession. At equilibrium, the disc mid-plane is close to being coplanar with the inner binary, but on a slightly misaligned plane with respect to the orbital planes of the stars. Remarkably, the observed dynamical state of the disc cannot be explained using binary system models alone.

• 4.

  Given the likely orbital configuration of V892 Tau and assuming planets will form in the disc, these circumbinary planets would become unstable in the very inner parts and in the outer parts of the disc – as soon as the gas dissipates. Any surviving planet should have an intermediate semi-major axis between 40 au and 80 au and is expected to undergo significant oscillations in eccentricity and inclination due to the interactions with the external companion.

Based on our investigation, V892 Tau constitutes the archetype of a hierarchical triple system with a CBD, standing out as a unique test bed for planet formation in multiple stellar systems. However, to thoroughly infer the available formation channels of circumbinary planets in such systems, detailed studies on dust dynamics and growth are required. In addition, future near-IR observations would be key to further characterise the disc content and morphology, which would add complementary and independent constraints on the disc structure. By combining this information to gas-dust hydrodynamical simulations coupled to dust growth algorithms, we could in principle identify the possible regions of planet formation in V892 Tau. This would pave the way for a global understanding of planet formation and dynamics in the context of multiple systems.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union Horizon Europe programme (grant agreement No. 101042275, project Stellar-MADE). The Geryon2 cluster housed at the Centro de Astro-Ingeniería UC was used for the calculations performed in this paper. The BASAL PFB-06 CATA, Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon/Geryon2 cluster. Some of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. A.R. has been supported by the UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1 and by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (RISE DUST-BUSTERS project). R.N. acknowledges support from UKRI/EPSRC through a Stephen Hawking Fellowship (EP/T017287/1). Support for AI was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51532.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2021.1.01137.S., ADS/JAO.ALMA#2013.1.00498.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The data underlying this article will be shared on reasonable request to the corresponding author. The code PHANTOM used in this work is publicly available at https://github.com/danieljprice/phantom. The DISCMINER code used in this study is open source and can be found at https://github.com/andizq/discminer. This work makes use of SPLASH (Price 2007), NUMPY (Harris et al. 2020) and MATPLOTLIB (Hunter 2007). We also acknowledge the use of SARRACEN (Harris & Tricco 2023). A.A. would like to extend his sincere appreciation to Jean-François Gonzalez and Enrico Ragusa for their fruitful discussions throughout the course of this research. We would like to express our profound gratitude to the anonymous referee for their insightful comments and suggestions which have greatly contributed to the improvement of this work.

Appendix A Data combination

The data on which this paper is based come from a combination of data from two observational programs. A first dataset (OLD) has been acquired in the context of the ALMA programm 2013.1.00498.S in 2013, and analysed in details in L21. A second one has been acquired in end-2021 with longer baselines (LB) being covered during the ALMA programm 2021.1.01137.S. A third dataset has been acquired mid-2022 in the context of the same ALMA programm and covering shorter baselines (SB).

The data reduction and calibration are described in Section 2.1 and we add details about the self-calibration process of the data in the following. The three raw datasets were downloaded from the ALMA archive and calibrated using the available calibration scripts. At first, we reduced and self-calibrated each dataset independently. The line channels were flagged and the SPW averaged together to create continuum MS. Using the gaincal and applycal, self-calibration rounds in phase were conducted on the continuum data. The first calibration tables were created for combined SPWs, combined polarisation and a infinite solution interval to maximise the S/N. In the following rounds, the SPWs, polarisations were progressively decombined and the solution interval progressively turned down to the integration time. Between each rounds, the RMS was measured on images produced with tclean.

We ensured that the RMS was decreasing between the rounds without any decrease of the total flux and peak intensity.

Second, we created a common model for the three resulting continuum datasets using the tclean task. Before doing so, the observations were centred together with the phaseshift task using the centre of the cavity as a common reference. The different MS were also rescaled to a common flux reference using the DSHARP utilities ². From this model, additional self-calibration rounds in phase were performed on the individual MS. It ensured that the three datasets were properly centred together. The calibration were stopped before the RMS stopped decreasing between the rounds. The gaincal parameters were turned down to a solution interval of two times the integration time, combined SPWs and independent polarisations. The final continuum image was created using tclean from the three resulting MS of the previous process.

Third, ¹²CO channels were separated from the original MS. The continuum level was subtracted from each MS using the task uvcontsub. The resulting gas MS were scaled to a common flux reference using the same correction factors that the continuum MS. The calibration tables of the previous self-calibration rounds were applied to the corresponding gas MS. Finally, the ¹²CO datacube was created using the tclean task, that ensures a spectral smoothing of the MS to a common spectral resolution of 0.5 km s⁻¹.

Figure A.1 shows how the combination of the datasets affects the final continuum image and ¹²CO moment 0 images. The combination of the dataset used in L21 with our newer data enhances the final S/N compared to our data (SB +LB) alone. This S/N increase allows the continuum to be imaged with a more robust parameter controlling the Briggs weighting when applying the tclean task, and thus for more angular resolution than in previous studies.

Fig. A.1 Comparison of the ALMA band 6 continuum (top) and 12CO (2–1) moment 0 (bottom) images resulting from the different data-sets combination. The left image in each panel has been created from SB, LB and the data-set used in L21. The middle image in each panel results from the combination of SB and LB only. The images on the right are residual maps, which are uniformly positive, highlighting the benefits of the data combination.

Appendix B Channel maps

Fig. B.1 12CO channel maps of V892 Tau. The contours shows the [3, 7]σ emission levels. The dashed line represents corresponding negative emission levels. The velocity of the emission channels is indicated at the top left. The synthesised beam is plotted at the bottom left of each channel.

Appendix C Gas emission fitting with Discminer

C.1 Initial conditions and results

The initial values of the MCMC sampling were chosen in agreement with the values derived in L21when possible, namely for the orientation and kinematic parameters. The initial values of the parameters defining the emission surfaces were set to initialise a thin disc, in agreement with the observations. The rest of the initial values were put in the middle of the prior range or at typical values. The reference radius r₀ was taken as 100 au. Priors were defined to limit the sampling of the parameter space and the computational cost. The range of priors for all the parameters was set to a limit the disc to a reasonable extent. Following the need of a precise model, we initialised 256 walkers with a total of 40000 steps. The chains reached convergence well before the end of the sampling and the resulting auto-correlation length was 106 steps. The first 6300 steps are considered as a burn-in phase where the walkers explore the parameter space. After that, the walkers reached convergence and remained concentrated around the final value. The best-fit parameters are taken as the median of the posterior distributions. The errorbars are the 16 and 84 percentiles of the posterior distributions. Table C.1 summarises the results of the fit.

Fig. C.1 Same as Figure B.1 for DISCMINER best-fit channel maps.

C.2 Rotation curve

DISCMINER allows velocity radial profiles from the results of the fitting procedure to be computed. The rotation curve of the data can be compared to the best-fit model velocity profile. From the resulting residuals and and their projection in the skyplane, the velocity deviation to a Keplerian profile ∆v_ϕ, the velocity of radial flows v_r and the velocity of vertical flows v_ɀ can all be computed as a function of radius. In Figure C.2, we show these quantities from our best-fit model to the ¹²CO (2–1) emission of V892 Tau. The disc is dominated by Keplerian rotation, although slight offsets can be observed in the innermost and outermost parts of the disc. In the inner parts, this mainly translates into a rise of v_r, potentially triggered by interactions with the inner binary. The outer disc is characterised by inward radial flows, vertical motions of the gas and slight deviations from the keplerian rotation speed. These patterns are consistent with interactions with V892 Tau NE and support the hypothesis of a bound triple system.

Fig. C.2 Radial velocity profiles of the 12CO (2–1) emission of V892 Tau fitted by DISCMINER. Top : rotation curve. Bottom : Offset to Keplerian rotation (blue), radial (yellow) and vertical (black) velocities as a function of radius.

Appendix D Manual rotation of the disc before post-processing

Discs in multiple systems are expected to precess due to the gravitational torques applied by misaligned stars (Papaloizou & Terquem 1995). We observe that precession in our hydrodynam-ical models of the disc of V892 Tau : the disc is initialised at its observed inclination and position angle (PA), but these values evolve during the simulation because of the gravitational perturbations induced by the inner binary and the external companion. Values of inclination and PA of the disc at the end of each simulation can be found in Table D.1 (see Fig 8 for a visual representation). However, for a proper comparison with the observations, the simulated discs should have the same orientation as the observed disc. But, it would be too computationally expensive to simulate the disc during a whole precession timescale, at the end of which the disc would have come back to its original PA while the inclination would have slightly changed. We tackle this issue by manually rotating the simulated discs back to their observed values of inclination and PA before going into their post-processing into synthetic observations. This assumes that the disc properties other than inclination and PA, and the interactions with the stars are independent of the disc orientation. We are aware of the discrepancy between the observed orientation of the disc and the orientation at the end of the simulations. But, we argue that this discrepancy is reasonably small, which allows us to work in the context of the previous hypothesis. This way, rotating the disc to the observed orientation will allow for an adequate comparison of the observations and the synthetic data.

Appendix E Inclination sign of the disc

Usually, it is straightforward to determine the closest side of a disc when its upper and lower emission surfaces are both detected (e.g. IM Lup, Pinte et al. (2018)). Here and in previous studies of the system, the lower surface of the V892 Tau disc is not clearly detected, making the sign of the inclination uncertain. According to the fitted linewidths maps of the disc displayed on Figure 6 last row, the Keplerian model is able to reproduce the data well. Nonetheless, the linewidths seems larger in the NW side of the disc than in the SE side. When compared to the model, the NW side shows an excess in linewidth (see residual map in Fig. 6 last row right panel). Lines are expected to be seen broader in the closest side of the disc to the observer due to the greater contribution of the lower emission surface (e.g. MWC480, Izquierdo et al. (2023)). Applied to V892 Tau, this could indicate that the NW side of the disc is the closest side to the observer and that the disc inclination is negative. Trying to confirm that result, we ran an MCMC sampling to fit the disc observations with a model with a negative inclination sign with the same method as described in Sections 3.1 and C. In that model, the closest side to the observer is the NW side of the disc. The resulting fitted values were similar to the original model apart from the inclination sign. It lead to a disc model quite identical to the original one in the outer disc, as showed by the linewidth map of each model on Figure E.1 below. This could be explained by the disc being significantly flat and its lower emission surface not being detected, which doesn’t allow to conclude on the true orientation of the disc. Differences still can be seen between the two models in the central parts of the disc, for which our data are not reliable due to cloud absorption.

In general, the far side of the disc is observed with a larger brightness temperature than the near side (Law et al. 2023), which is in line with our synthetic observations resulting from the simulations where the closest side to the observer is assumed to be the SE side. According to our ¹²CO observations, the NW side is brighter in average than the SE side by approximately 13%. This agreement between the simulations and the observations is in favour of the SE side being the nearside of the CBD. Moreover, the temperature asymmetry plotted on Figure 7 and heading toward the NW cloud also trace the flared far side of the disc, since similar patterns have already been observed in discs with a known orientation (see RXJ 1615 in Wölfer et al. (2023), LkCa 15 and HD 34282 in Law et al. (2023)). Considering the azimuthal asymmetry in the continuum disc arises from warm dust emission, this emission may come from the visible inner rim of the disc far side (Ribas et al. 2024). Such a scenario would indicate that the disc near side is the SE side of the disc.

The disc thickness and faint lower surface make the conclusion on the disc true inclination difficult. Nonetheless, we present the previous results as tentative evidences for the SE side of the disc being the closest side to the observer. Future scattered light observations or high resolution polarisation maps could help to characterise the dust scattering properties which would bring additional constrains on the disc true orientation.

Fig. E.1 Comparison of the fitted linewidth models from DISCMINER. The original model is represented in the left panel, the model with a flipped inclination sign in the middle panel and the residuals between the two are plotted in the right panel. A mask defined by the 3σ contour of the 12CO (2–1) emission has been applied to the models.

Appendix F Binary-only simulations

We performed additional hydrodynamical SPH simulations of the V892 Tau system. In these simulations, we model the system as a binary system with a disc. The disc parameters are the same as described in Section 3.2. Then, we consider different setups where the companion star is ignored, or where the inner binary is merged into a single star with an equivalent mass. These setups are :

• the inner binary surrounded by the disc (IB).

• the outer binary initially misaligned of 30° with respect to the disc (OB30).

• the outer binary initially misaligned of 60° with respect to the disc (OB60).

The orbital parameters of the binary star in IB are the same as the orbital parameters of the inner binary in the triple system simulations. The orbital parameters of the binary star in OB30 are the same as the orbital parameters of the outer binary in i30. The orbital parameters of the binary star in OB60 are the same as the orbital parameters of the outer binary in i60. These simulations counted 10⁶ SPH particles and ran for a total time of t = 25 P_out, where P_out is the outer binary period in the triple system simulations.

The oscillations of the modelled CBD are triggered both by the inner binary and the outer binary. Yet, the contribution of each component to the dynamics of the disc is unclear : the simulated CBD is found in an intermediate plane between the inner and outer binary orbital planes, its global inclination oscillates without reaching a coplanar nor polar configuration and it precesses non-linearly. Simulations of the V892 Tau system modelled as a binary star could help to disentangle the dynamical effects of the inner and of the outer binaries on the disc.

Figure F.1 shows the inclination i and the PA Ω of the disc as a function of time in the simulations i30, i60, OB30, OB60, and IB.

In the inner binary only simulation, the CBD remains copla-nar with the inner binary and its PA remains at the initialised value. The PA of the disc follows a decreasing trend in all the other simulations. During the very first orbits, the disc has the same linear precession rate in the setups i30 and OB30, and i60 and OB60 respectively. Then, the PA of the disc starts to oscillates in the triple system simulations, while the precession stays linear in the outer binary systems OB30 and OB60. The precession slope in those setups is steeper than the average slope in the triple setups, but asymptotically reproduces the precession rate of the disc in triple systems during the first orbits of the companion. The precession rate is similar in the OB30 and the OB60 setups. We expect that kind of linear precession for binary systems (e.g. Papaloizou & Terquem (1995)). From dynamical considerations, we expect the precession rate of the disc to be t_prec,0 = 144 P_out in the setup OB30 and t_prec,0 = 251 P_out in OB60 (Terquem 1998; Bate et al. 2000). We measure a precession rate of t_prec = 189 ± 1 P_out and t_prec = 195 ± 1 P_out respectively, which is in rough agreement with that prediction.

We already saw that the disc inclination follows damped oscillations in the triple system setups i30 and i60. In the OB30 setup, the inclination of the disc drops with time, but that drops is not linear nonetheless. The OB30 disc could then align with the orbital plane of the binary, but the simulation stopped too early to clarify this statement. In the OB60 model, the inclination of the disc rises steadily in a linear way with a slope of 5.6 × 10⁻² ± 0.2 × 10⁻² P_out. That particular behaviour could be explained by the disc starting to align with the binary orbital plane. The timescale of this alignment is larger that in the OB30 case due to the higher mutual inclination of the disc to the orbital plane, which results in a smaller torque applied by the companion on the disc (Bate et al. 2000).

The three additional simulations IB, OB30 and OB60 show that the dynamical behaviour of the V892 Tau disc can not be reproduced by a binary system only. The triplicity of the system, which means both the inner binary and outer binary, is responsible for the inclination oscillations and the peculiar precession of the CBD.

Fig. F.1 Mean inclination (left) and mean position angle (right) of the disc as a function of time for the binary star setups and the triple systems i30 and i60. The top row shows the discs in IB, OB30, and i30 and the bottom row shows the discs in IB, OB60, and i60.

Appendix G Other N-Body simulations

Fig. G.1 Same as Figure 15 but for the ref orbital configuration.

Fig. G.2 Same as Figure 15 but for the e05 orbital configuration.

Fig. G.3 Same as Figure 15 but for the i30 orbital configuration.

Fig. G.4 Same as Figure 15 but for the i60 orbital configuration.

References

All Tables

All Figures

	Fig. 1 V892 Tau is located in a clustered sub-region of the Taurus molecular cloud. The background image from the Herschel/SPIRE survey (Griffin et al. 2010; André et al. 2010) is shown on a linear intensity scale. Magenta points show positions of pre-main-sequence stars in Taurus taken from Kenyon et al. (2008).
In the text

Fig. 2 ALMA band 6 continuum emission map (top) and depro-jected azimuthally averaged intensity profile (bottom) at 224 GHz of V892 Tau. In the top image, contour levels are at 10σ, 50σ, and 150σ. In the bottom plot Rcav, Rpeak, and R90,cort correspond to the radius of the cavity, of the peak of the ring, and of the disc, respectively. The shaded orange area represents to the standard deviation of each annulus divided by the square root of the beam numbers along the annulus at each radial bin. The shaded grey area represents the major axis of the synthesised beam.

In the text

Fig. 3 Projection of the initialised orbits of each tested setup in the sky-plane xy and in a plane perpendicular to it, namely xz. The coordinates are centred on the primary. The CBD is represented by the shaded grey area. The projected position of V892 Tau NE in the skyplane is represented by a star and by a black line in the xz plane. The orbit of the inner binary (IB) is plotted in blue and the other colours correspond to the tested orbits for V892 Tau NE.

In the text

Fig. 4 Channel maps of the 12CO (2–1) emission. The DISCMINER best-fit model is plotted in colours and the data is represented by the white contour lines at [3,7]σ. Negative emission at the same significance is indicated by dashed lines. The synthesised beam is plotted at the bottom left of each channel map, while the corresponding velocity is plotted at the top left. On the bottom left channel, an arrow indicates the direction to the companion star V892 Tau NE.

In the text

	Fig. 5 Comparison of channel maps at 3.5 km s−1 (left), 11 km s−1 (middle), and 12.5 km s−1 (right) of the observed emission (top row) and of the synthetic emission (following rows, in order: ref, e05, i30, i60, and ei60). The contours represent the [3,7]σ emission levels. The beam size is represented by the grey ellipses.
In the text

Fig. 6 12CO (2–1) moment maps of V892 Tau. The integrated intensity, the peak temperature, the peak velocity, and the line width correspond to the first, second, third, and last rows, respectively. The first column shows moment maps computed from the data emission channels, the second column from the model and the last column shows the residuals in the detected disc area computed from the two previous maps. A mask defined by the 3σ contour of the data has been applied to the model before comparison with the data.

In the text

Fig. 7 Sketch of the features identified in the CBD of V892 Tau from the residual maps shown in Fig. 6. The black line corresponds to the 3σ contour of the disc emission. The indicated high brightness residuals come from the moment 0 residuals. The spiral arms are identified on the velocity residual map, from which only the red part is plotted and the centre masked for clarity reasons. The temperature asymmetry coming from the peak temperature residual map is overlaid in dark brown. The arrow in the top left indicates the direction of V892 Tau NE.

In the text

Fig. 8 Rendered plots of the density spatial distribution integrated along the line of sight of the different simulated setups. The reƒ, e05, i30, i60, and ei60 setups are shown from left to right. The top row shows the system after one orbit of the outer companion. The bottom row displays the state of the system after approximately 50 orbits of the outer companion, corresponding to 0.2 Myr and 28 400 orbits of the inner binary. The white dots indicate the sink particles representing the stars.

In the text

Fig. 9 12CO (2–1) integrated intensity map of the observational data (top left) compared to the synthetic maps built from the simulations reƒ (top middle), e05 (top right), i30 (bottom left), i60 (bottom middle), and ei60 (bottom right). The discs were rotated back in the observed disc plane to allow for a proper comparison to the observations. The beam sizes are plotted in the bottom left corners.

In the text

	Fig. 10 Same as Fig. 9 but for the 12CO (2–1) peak temperature map.
In the text

	Fig. 11 Same as Fig. 9 but for the 12CO (2–1) emission velocity map.
In the text

	Fig. 12 Azimuthally averaged radial profile of the 12CO intensity for the observations and the simulations. The vertical dashed lines indicate the R90% radius for the continuum, the gas emission, and the synthetic observations. The shaded grey area corresponds to the major axis of the synthesised beam.
In the text

	Fig. 13 Tilt (top) and PA (bottom) evolution of the disc for the i30 (left), i60 (middle), and ei60 (right) setups. The data are represented by the coloured points while the fit is shown by the red curve. The tilt evolution was modelled by damped oscillations. The PA evolution was modelled by oscillations on the top of a linear precession. The best-fit parameters can be found in Table 3.
In the text

	Fig. 14 Tilt radial profiles of various simulations at t = 10 Pout (left), t = 30 Pout (middle), and t = 50 Pout (right). The orange, blue, and teal dots correspond to i30, i60, and ei60, respectively. The 1σ errors are indicated by the shaded areas.
In the text

Fig. 15 Evolution of the semi-major axis (top), eccentricity (middle) and mutual inclination with the inner binary (bottom) of the test particles in the ei60 N-body setup. Each colour represents a different test particle that started at a given semi-major axis. The semi-major axis of a particle falls to 0 once ejected from the system. The particle starting with a = 20 au was ejected from the system before the completion of the first orbit of the outer binary.

In the text

Fig. A.1 Comparison of the ALMA band 6 continuum (top) and 12CO (2–1) moment 0 (bottom) images resulting from the different data-sets combination. The left image in each panel has been created from SB, LB and the data-set used in L21. The middle image in each panel results from the combination of SB and LB only. The images on the right are residual maps, which are uniformly positive, highlighting the benefits of the data combination.

In the text

	Fig. B.1 12CO channel maps of V892 Tau. The contours shows the [3, 7]σ emission levels. The dashed line represents corresponding negative emission levels. The velocity of the emission channels is indicated at the top left. The synthesised beam is plotted at the bottom left of each channel.
In the text

	Fig. C.1 Same as Figure B.1 for DISCMINER best-fit channel maps.
In the text

	Fig. C.2 Radial velocity profiles of the 12CO (2–1) emission of V892 Tau fitted by DISCMINER. Top : rotation curve. Bottom : Offset to Keplerian rotation (blue), radial (yellow) and vertical (black) velocities as a function of radius.
In the text

	Fig. E.1 Comparison of the fitted linewidth models from DISCMINER. The original model is represented in the left panel, the model with a flipped inclination sign in the middle panel and the residuals between the two are plotted in the right panel. A mask defined by the 3σ contour of the 12CO (2–1) emission has been applied to the models.
In the text

	Fig. F.1 Mean inclination (left) and mean position angle (right) of the disc as a function of time for the binary star setups and the triple systems i30 and i60. The top row shows the discs in IB, OB30, and i30 and the bottom row shows the discs in IB, OB60, and i60.
In the text

	Fig. G.1 Same as Figure 15 but for the ref orbital configuration.
In the text

	Fig. G.2 Same as Figure 15 but for the e05 orbital configuration.
In the text

	Fig. G.3 Same as Figure 15 but for the i30 orbital configuration.
In the text

	Fig. G.4 Same as Figure 15 but for the i60 orbital configuration.
In the text