Context. Stars and planets form in regions of enhanced stellar density subjecting protoplanetary discs to gravitational perturbations from neighbouring stars. Observations in the Taurus starforming region have uncovered evidence of at least three recent stardisc encounters that have truncated discs HVDO Tau RW Aurigae and UX Tau raising questions about the frequency of such events. Aims. We aim to assess the probability of observing truncating stardisc encounters in Taurus. Methods. We generated a physically motivated dynamical model including binaries and a spatialkinematic substructure to follow the historical dynamical evolution of the Taurus starforming region. We used this model to track stardisc encounters and the resulting outer disc truncation over the lifetime of Taurus. Results. A quarter of discs are truncated below 30 au by dynamical encounters but this truncation mostly occurs in binaries over the course of a few orbital periods on a timescale 0.1 Myr. Nonetheless some truncating encounters still occur up to the present age of Taurus. Strongly truncating encounters ejecting 10 percent of the disc mass occur at a rate 10 MyrSUP1SUP sufficient to explain the encounter between HV and DO Tau 0.1 Myr ago. If encounters that eject only 1 percent of the disc mass are responsible for RW Aurigae and UX Tau then they are also expected with encounter rate SUBencSUB  100200 MyrSUP1SUP. However the observed sample of recent encounters is probably incomplete since these examples occurred in systems that are not consistent with a random drawing from the mass function. One more observed example would statistically imply additional physics such as replenishment of the outer disc material. Conclusions. The marginal consistency of the frequency of observed recent stardisc encounters with theoretical expectations underlines the value of future large surveys searching for external structures associated with recent encounters. The outcome of such a survey offers a highly constraining novel probe of protoplanetary disc physics.
['planets and satellites: formation', 'protoplanetary disks', 'binaries: general', 'circumstellar matter', 'stars: formation', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
Running with the bulls The frequency of stardisc encounters in the Taurus starforming region
{'The frequency of star-disc encounters in the Taurus star forming region': "Andrew J. Winter, 1, ⋆ Myriam Benisty, 1 , 2 Linling Shuai, 3 , 4 Gaspard Dûchene, 2 , 5 Nicolás Cuello, 2 Rossella Anania, 6 Corentin Cadiou, 7 Isabelle Joncour 2  \n- 1 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06300 Nice, France\n- 2 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France\n- 3 Astronomy Department, University of Michigan, Ann Arbor, MI 48109, USA\n- 4 Department of Astronomy, Xiamen University, 1 Zengcuoan West Road, Xiamen, Fujian 361005, China\n- 5 Astronomy Department, University of California Berkeley, Berkeley CA 94720-3411, USA\n- 6 Dipartimento di Fisica 'Aldo Pontremoli', Universita degli Studi di Milano, via Celoria 16, Milano, 20133, Italy\n- 7 Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden  \nReceived September 15, 1996; accepted March 16, 1997", 'ABSTRACT': 'Context. Stars and planets form in regions of enhanced stellar density, subjecting protoplanetary discs to gravitational perturbations from neighbouring stars. Observations in the Taurus star-forming have uncovered evidence of at least three recent, star-disc encounters that have truncated discs (HV / DO Tau, RW Aurigae, UX Tau), raising questions about the frequency of such events.  \nAims.  \nWe aim to assess the probability of observing truncating star-disc encounters in Taurus.  \nMethods. We generate a physically motivated dynamical model including binaries and spatial-kinematic substructure to follow the historical dynamical evolution and stellar encounters in the Taurus star forming region. We track the star-disc encounters and outer disc radius evolution over the lifetime of Taurus.  \nResults. A quarter of discs are truncated below 30 au by dynamical encounters, but this truncation mostly occurs in binaries over the course of a few orbital periods, on a time-scale ≲ 0 . 1 Myr. Nonetheless, some truncating encounters still occur up to the present age of Taurus. Strongly truncating encounters (ejecting ≳ 10 percent of the disc mass) occur at a rate ∼ 10 Myr -1 , su ffi cient to explain the encounter between HV and DO Tau ∼ 0 . 1 Myr ago. If encounters that eject only ∼ 1 percent of the disc mass are responsible for RW Aurigae and UX Tau, then they are also expected with encounter rate Γ enc ∼ 100 -200 Myr -1 . However, the observed sample of recent encounters is probably incomplete, since these examples occurred in systems that are not consistent with random drawing from the mass function. One more observed example would statistically imply additional physics, such as replenishment of the outer disc material.  \nConclusions. The marginal consistency of the frequency of observed recent star-disc encounters with theoretical expectations underlines the value of future large surveys searching for external structures associated with recent encounters. The outcome of such a survey o ff ers a highly constraining, novel probe of protoplanetary disc physics.  \nKey words. protoplanetary discs - planet formation - star forming regions', '1. Introduction': "Planet formation proceeds in a 'protoplanetary disc' of dust and gas over a time-scale of ∼ 3 Myr (e.g. Haisch et al. 2001). During this time, nascent planetary systems typically inhabit star forming regions with a local stellar density that far exceeds the average in the galactic neighbourhood (e.g. Lada & Lada 2003). Neighbouring stars in these regions can feedback on the planet formation process in a variety of ways. These may include: external irradiation driving thermal winds (Winter & Haworth 2022, and references therein), chemical enrichment (Bastian et al. 2013; Lichtenberg et al. 2016; Parker & Schoettler 2023), late-stage gas in-fall (Dullemond et al. 2019; Ku ff meier et al. 2020, 2021, 2023) and star-disc encounters (Cuello et al. 2019, 2020, 2023). While all of these processes are of great interest for understanding the diversity of the observed exoplanet population, in this work we focus on the latter phenomenon. 'Star-disc encounters' refer to gravitational perturbations expe-  \nenced by a protoplanetary disc during the close passage of a neighbouring star.  \nVarious phenomena such as stellar accretion outbursts (Pfalzner 2008; Forgan & Rice 2010; Vorobyov et al. 2021; Dong et al. 2022) and spiral arms in protoplanetary discs (e.g. De Rosa & Kalas 2019) and free-floating planets (e.g. Vorobyov et al. 2017) - possibly including enigmatic binary planet-mass objects in the Orion Nebula cluster (Pearson & McCaughrean 2023; Wang et al. 2023; Portegies Zwart & Hochart 2023) may all be feasibly attributed to stellar encounters. There remain alternative explanations; for example, spiral arms may be produced by gravitational instability (e.g. Douglas et al. 2013; Meru et al. 2017; Baehr & Zhu 2021) or companion stars / brown dwarfs / planets (e.g. Dong et al. 2015; Ren et al. 2020). A search for flyby candidates that may have generated spiral arms by Shuai et al. (2022) revealed no evidence that nearby stars recently had a su ffi ciently close approach. However, such an e ff ort is challenging due to proper motion uncertainties and incompleteness of reliable measurements; particularly if a perturber is  \nvery low mass or is itself a binary. To assess the role of star-disc encounters for protoplanetary disc evolution, we must estimate the rate of encounters in dynamically evolving star forming regions.  \nThe rate of stellar encounters, and their role for planet formation, is dependent on the local stellar density of the star forming region. The stellar mass density ρ ∗ of star forming regions can vary dramatically in the range 1 M ⊙ pc -3 ≲ ρ ∗ ≲ 10 6 M ⊙ pc -3 . The latter limit represents the most extreme stellar densities in the galaxy, such as in the cores of globular clusters that remain bound for a Hubble time (e.g. Krumholz et al. 2019). Intermediate densities ρ ∗ ∼ 10 3 M ⊙ pc -3 are typical in the cores of open clusters that may survive against galactic tides over 100 Myr or Gyr time-scales. However, possibly the most common environments in which stars and planets form are the low density regions that produce loose 'associations', which are globally unbound against galactic tides. The average density in such regions has been assumed to be too low to induce frequent tidal encounters throughout the disc lifetime (e.g. Winter et al. 2018c).  \nTaurus is an example of a (globally) low density region, which does however host convincing evidence of at least three recent close star-disc encounters. These cases are RW Aurigae (Cabrit et al. 2006; Dai et al. 2015; Rodriguez et al. 2018), HV and DO Tau (Howard et al. 2013; Winter et al. 2018a) and UX Tau (Zapata et al. 2020; Ménard et al. 2020). Each of these systems exhibits significant external structure which appears to be explained by recent flybys, su ffi ciently close such that they were capable of unbinding disc material. In an unstructured star forming region at the average density of Taurus, the probability of a one-o ff random encounter over the disc lifetime is vanishingly small. This puzzle may appear to be partially resolved for some cases if the system is bound, such that close encounters occur periodically over the binary orbit. However, in isolation this is not a su ffi cient explanation. While repeated encounters in a multiple system can drive spiral arms in the disc (e.g. Alaguero et al. 2024), the disc should also be rapidly truncated such that on subsequent close approaches the gravitational pertubation no longer produces large extended tidal tails (e.g. Ménard et al. 2020). It is therefore necessary to perturb star-disc systems in Taurus stochastically, and do so over its ∼ 1 -3 Myr lifetime. In this manuscript, we aim to answer the question: do we expect su ffi -cient star-disc encounters in Taurus to explain the frequency of observed recent encounters?  \nThe role of stellar encounters has been the focus of numerous studies focusing both on the evolution of protoplanetary discs (e.g. Clarke & Pringle 1993; Ostriker 1994; Pfalzner et al. 2005; Winter et al. 2018b) and mature planetary systems (e.g. Spurzem et al. 2009; Shara et al. 2016; Winter et al. 2022; Li et al. 2024). Studies using N-body simulations often perform parameter studies or implement scaling relations to model star forming regions. For example, one approach is to adopt the putative mass-radius (e.g. total stellar mass M c, half-mass radius R c) relationship for star forming regions (e.g. R c ∝ M 1 / 2 c - Adams et al. 2006). However, the nature of this relationship varies depending on the sample / definition of mass and radius (Pfalzner & Govind 2021), and exhibits significant scatter. Mature / massive star clusters also do not appear to follow this relationship (Krumholz et al. 2019), and particularly for low mass / density regions it is challenging to unambiguously define an individual star forming region. In addition, even within a well-defined 'individual' star forming region, internal structure has been shown to have a strong influence on the role of encounters (e.g. Craig & Krumholz 2013; Parker 2023). Assessing encounter rates therefore requires quantitatively matching present day position-velocity structure. Yet it  \nis not clear how well the widely adopted method for generating fractal initial conditions, sampling hierarchical boxes (Craig & Krumholz 2013), reproduces the observed spatial and kinematic structure in star forming regions.  \nStructure in giant molecular clouds is set by the turbulent fragmentation, from which the power-spectrum and Mach number determine the mass density distribution (e.g. VazquezSemadeni 1994; Padoan et al. 1997). An ideal approach to studying the role of dynamical encounters in young star forming regions is therefore to model the star formation process directly through hydrodynamic simulations. Such simulations have shown that encounters are common during early disc evolution, possibly determining the initial distribution of protoplanetary disc radii (Bate 2018). However, these experiments are computationally expensive, following star formation only for ∼ 0 . 1 Myr time-scales. This makes parameter studies or tailored modelling of individual star forming regions in this way impracticable.  \nHere we present a complementary approach to the above works, targeted at generating N-body initial conditions tailored to match young star forming regions. We achieve this by simulating a Taurus-like star forming region, with physically and empirically motivated initial conditions, including binaries. We draw stellar positions and velocities from an empirically constrained power spectrum, reflecting how turbulent energy is distributed across di ff erent scales. Our main aim is to develop a dynamical model to track the history of stellar encounters in the Taurus star forming regions. In doing so, we will assess whether the rate of disc-truncating star-disc encounters in Taurus is su ffi cient to produce the three known examples: HV / DO Tau, RW Aurigae and UX Tau. This goal requires ensuring that our dynamical model closely reproduces the present day spatial and kinematic structure in Taurus. We therefore review these structural properties in Section 2, which we then adopt for benchmarking our dynamical model. We discuss our approach for initialising initial conditions and dynamically evolving the model in Section 3, in which we also draw comparisons with the structure metrics introduced in Section 2. We discuss the rate of truncating encounters over time in our simulation in Section 4, quantifying the degree to which the observed examples of recent star-disc encounters are statistically expected. We summarise our conclusions in Section 5.", '2.1. Aim': 'An empirically motivated dynamical model for Taurus requires accurately quantifying kinematic substructure. If we had arbitrarily accurate and complete data for the 3D positions and velocities of all the stars in Taurus, then it would be trivial to use these data to generate N-body initial conditions to compute the future evolution of Taurus. However, as discussed in Section 2.2, while we have a fairly complete census of stars in terms of their projected spatial distribution, kinematic data is far more limited. Given the e ff ects of multiplicity and extinction, parallax measurements also typically have associated uncertainties that make them impractical for use in setting initial conditions. In addition, we are interested in quantifying the frequency of encounters in the past, while we do not have direct measurements of the early spatial-kinematic stellar configuration in Taurus. For these reasons, we require metrics that characterise kinematic substructure. These metrics both guide our choice of initial conditions and o ff er a benchmark comparison for our models.  \nIn the following, we first review the data for the stellar population in Taurus (Section 2.2). We then discuss how we quantify spatial structure in Section 2.3 and kinematic structure in Section 2.4.', '2.2. Data': 'We use the census of Taurus members by Luhman (2023), with astrometric data from Gaia DR3 (Gaia Collaboration et al. 2016, 2023; Babusiaux et al. 2023). This catalogue contains 532 members, with a high degree of completeness for spectral types earlier than M6-M7. When we consider proper motion di ff erences between neighbours (Section 2.4), we restrict the sample to the 271 that have a reliable astrometric solution by the canonical criteria that the Renormalised Unit Weight Error (RUWE 1 ) is smaller than 1 . 4 in Gaia DR3.  \nWecannot use the observed 3D positions directly to generate initial conditions for several reasons. Most obviously, the sample for which we have parallax is incomplete (417 / 532), and the typical uncertainties on the parallax correspond to a few pc. This means that small scale structure is not resolvable. We therefore generate our initial conditions parametrically, closely comparing with the observed structure in Taurus. We consider only the projected separations between stars when inferring the spatial and velocity structure of the region.', '2.3. Spatial substructure': "We aim to understand the role of close neighbours on disc evolution. A sensible metric to quantify structure is therefore the normalised pair separation function, which can be defined in two dimensions ˆ Σ pairs( ∆ R ). This is the averaged surface density of neighbours for any given star, normalised to unity when integrated over 2 π ∆ R d ∆ R . This surface density evolves over time, and we will therefore use this metric to ensure that we capture the time at which the dynamical state in our simulation is similar to that in Taurus.  \nTwo related but complementary metrics are the one- and twopoint correlation functions, Ψ ( ∆ R ) and ξ ( ∆ R ) respectively. These metrics have been applied by Joncour et al. (2017) and Joncour et al. (2018) to study the structure in Taurus. They are broadly defined as the excess of pairs with a given separation compared to a random, uniformly distributed population of stars in the same area. The one-point correlation Ψ is an excess of nearest pairs, while the two-point correlation represents the excess of all pairs. Of particular relevance for this work, Joncour et al. (2017) applied the one-point correlation function to demonstrate that approximately 40 percent of stars in Taurus are in ultra-wide binaries, with separations in the range 1 . 6 × 10 3 -5 × 10 4 au. We will use this inference to motivate our initial conditions, and the one- and two-point correlation functions to validate our model a posteriori . In particular, from the one-point correlation function Joncour et al. (2017) find a region of 'inhibition' (a smaller number of pairs than expected from random sampling) between ∼ 0 . 1 · and ∼ 0 . 5 · , while up to ∼ 0 . 2 · , Ψ ∝ ∆ R -1 . 5 . For the two-point correlation function, ξ > 1 out to ∆ R ≈ 2 · .", '2.4. Velocity substructure': 'Stellar velocities are inherited from the velocities of the material from which they form, and are thus dependent on the kinematic  \nFig. 1. Proper motion di ff erence as a function of angular separation on the plane of the sky for nearest neighbours in the Luhman (2023) sample for Taurus. Red data points are the observed sample, with uncertainties propagated from Gaia uncertainties. The black line shows the the velocity dispersion as a function of separation we adopt for our model, with a Keplerian component that dominates at small separations. The transition between the two power-laws is the break between the binary and individual systems ( ∼ 5 × 10 4 au, comparable to the galactic tidal radius for typical stars). The colour bar shows the corresponding MaxwellBoltzmann distribution, normalised for each separation. See text for details.  \n<!-- image -->  \nstructure of the parental molecular cloud. This velocity structure in a turbulent medium is driven by interacting waves that generate an energy cascade that is described by an energy spectrum E ( k ). Across a wide range of length scales λ , it can be approximated by a power-law (e.g. Elmegreen & Scalo 2004):  \nE ( k ) = P ( k ) k 2 ∼ σ 2 v , k k ∝ k -β ∝ λ β , (1)  \nwhere P ( k ) is the (three dimensional) power-spectrum, k is the wavenumber, σ v , k is the characteristic velocity for k .  \nA range of β values in the interstellar medium (ISM) have been inferred observationally. These range from the original estimate for the size-linewidth relation for molecular clouds yielding β ≈ 1 . 76 (Larson 1981) to larger β ≈ 2 . 3 (Münch 1958; Heyer & Brunt 2004, for example). For weakly compressible turbulence, the energy spectrum for density and velocity follow the same power-law. Qian et al. (2018) find β ≈ 2 on small scales ≲ 2 pc (suggestive of compressible turbulence) and β ≈ 5 / 3 in the range ∼ 5 -10 pc (suggestive of incompressible turbulence, see also Brunt 2010). Ha et al. (2022) estimated β ≈ 1 . 7 from the stellar population and β ≈ 1 . 8 from H α , both on ∼ 10 pc scales.  \nFor our purposes we adopt β = 2, appropriate for supersonic, rapidly cooling turbulence (Burgers 1948). The velocity dispersion relation we adopt is then:  \nσ v ( ∆ r ) = σ v , 0 ∆ r ∆ r 0 ! 0 . 5 , (2)  \nwhere we choose normalisation constants σ v , 0 = 0 . 6 km s -1 for ∆ r 0 = 1 pc. This yields a velocity dispersion similar to that empirically inferred for the youngest stars on galactic scales in the solar neighbourhood (e.g. Holmberg et al. 2009), while more importantly matching the dispersion in Taurus down to binary length scales. In Figure 1, we show the distribution of proper  \nmotion di ff erences for nearest neighbours as a function of their separation in Taurus. We also show the Maxwell-Boltzmann distribution:  \ngN D ( ∆ v | σ v ) = 2 ∆ v N D -1 Γ ( N D / 2) 1 4 σ 2 v ! N D / 2 exp -∆ v 2 4 σ 2 v ! (3)  \nin two dimensions ( N D = 2) for equation 2 with a Keplerian component appropriate for a star of mass m ∗ = 0 . 5 M ⊙ added in quadrature. This shows that the majority of neighbours follow the expected separation-velocity relation. Neighbours with a much larger relative velocity may have larger physical than projected separations. We will apply the size-velocity relation given by equation 2 to generate stellar velocities, and compare the appropriate Maxwell-Boltzmann distribution to the nearest neighbour relative velocities we generate in our model, both initially and after dynamical evolution.', '3.1. Overview': 'Our approach for simulating a Taurus-like star forming region is to use a physically and empirically motivated set of initial conditions, including binary systems. We then benchmark our model against the spatial-kinematic properties of the present day Taurus to ensure that the simulation reflects the observed dynamical state. With this model, we are able to extract encounters which we can convert into the protoplanetary disc radii evolution using analytic formulae fit to numerical experiments for truncating encounters.  \nIn the remainder of this section, we detail this process. First, we discuss our method for generating physically motivated initial conditions, starting with a gas density distribution in Section 3.3, which we convert to a stellar density distribution in Section 3.4. We then discuss implementing empirically motivated binaries in Section 3.6 and the appropriate velocity substructure in Section 3.8. We validate our dynamical model with respect to the observed dynamical state in Taurus in Section 3.9. Finally, we discuss our approach for extracting close stellar encounters from the simulation in Section 3.10 and the resultant disc evolution in Section 3.11.', '3.2. Nbody code': 'Throughout this work, we use NBODY6 ++ (Aarseth 2003; Spurzem 1999; Wang et al. 2015) to integrate the stars and binary systems under gravity for 3 Myr. Given our interest in the short term evolution of a low mass star forming region, we do not include stellar evolution. Nor do we include tidal binary circularisation, or an external potential. To capture the role of several neighbours, we insist on a short time-step factor for the irregular force polynomial 5 × 10 -3 . The parameter adjustment time in N-body units is 10 -4 in code time-units that are 3 . 22 Myr, at which the regularisation parameters are updated. The parameter and initial condition file for our fiducial model are available online.', '3.3. Lognormal gas density field': 'We aim to initially generate a physically motivated gas density distribution, from which to draw our stellar population. Numerical experiments have shown that the mass density ρ of isothermal turbulent flows is well approximated by a lognormal distribution (e.g. Vazquez-Semadeni 1994; Nordlund & Padoan 1999;', 'Gaussian field box selection': "Fig. 2. Illustration of how we select a box from which to generate a lognormal density field from a larger scale, low resolution lognormal density field. We illustrate the underlying density grid by drawing faint black points with density proportional to the local density field. Cyan points indicate our selection of the top 20 density maxima from the density field. The red star indicates the position of the adopted box centre. The dashed red lines indicate the boundary of the new zoom-in box.  \n<!-- image -->  \nOstriker et al. 1999). In order to generate a lognormal gas distribution, we must first generate a Gaussian density field over a 3D grid. Generating a density field that we can use to produce initial conditions for the N-body simulation is not trivial, since a Gaussian random field may exhibit peaks close to the edges of the grid. These peaks would be the site of high stellar density (or 'NESTS' - Joncour et al. 2018), and may be partially excluded by the grid. Indeed, the centre of the grid can represent a 'void' (underdensity) in the Gaussian field, which would undermine our goal to explore the interactions within and between stars in overdensities. Since we also need a grid that covers a large dynamical range (from the size scale of Taurus, down to the wide-binary scale), it is not practicable to simply choose a subset of the complete Gaussian random field with a very high resolution. It is therefore useful to apply an approach that ensures that we can centre the grid on an overdensity.  \nWhile we could attempt to do this by multiplying the density field by a centrally concentrated profile (such as a 3D Gaussian or Plummer profile), this would fundamentally change the nature of the density field and corresponding power spectrum. To avoid this we adopt a two-stage process, inspired by cosmological zoom-in simulations and applying the publicly available genet IC 2 code (Stopyra et al. 2020, 2021). This code generates a zoom-in or splice (Cadiou et al. 2021) self-consistently within a surrounding Gaussian random field by Fourier-space filtering. The net e ff ect is that we are able to generate an initial, coarse Gaussian field, and then zoom in on an overdensity that represents the centre of our model for Taurus.  \nOur approach for this is represented graphically in Figure 2. We first generate a Gaussian random field over a 256 3 grid with  \nside length 100 pc. We adopt a power spectrum index of -3. We do not quantitatively fit for this index, but select a posteriori between indices -2, -2 . 5, -3 and -4, finding that -3 reproduces a stellar density structure that matches the observed structure. Empirically, this index is comparable to that inferred in Taurus by Brunt (2010).  \nFrom the large scale density field, we then locate local maxima over 5 pc length scales using the ndimage.maximum\\_filter algorithm from S cipy (Virtanen et al. 2020). We choose the 20 greatest local maxima, and select the closest to the centre of the grid, ensuring the zoom-in grid will be within the range of the initial grid. We then use genet IC to generate a new Gaussian density field on a grid that is smaller by a factor three (box side length 33 pc), and with 1024 3 grid cells. This defines a field which has a high density centre, and a resolution of 0 . 032 pc, or ∼ 7000 au, down to length scales comparable to wide binary separations. Finally, we enforce a dispersion σ ln ρ = ln GLYPH<16> 1 + 0 . 25 ⟨ v 2 ⟩ / c 2 s GLYPH<17> = 5 on the logarithmic density field. This is approximately the expected dispersion given a sound speed c s = 0 . 2 km s -1 and the root mean square velocity ⟨ v 2 ⟩ 1 / 2 = √ 3 σ v ∼ 5 km s -1 for ∆ r ∼ 30 pc, according to equation 2.", '3.4. Stellar spatial distribution': 'We now wish to draw a stellar spatial distribution from our underlying gas density profile. To do this, we must consider which regions of the cloud are able to rapidly collapse to form stars. We consider the time-scale for this collapse from rest:  \nτ ff = 1 . 52 ρ 100 M ⊙ pc -3 ! -1 / 2 Myr . (4)  \nTo infer a free-fall time across our grid, we require an absolute density scale, for which we assume a total mass within the box of 10 4 M ⊙ , corresponding to the approximate total gas mass in the Taurus complex (Goldsmith et al. 2008). We then impose a constraint that grid cells must have τ ff < 1 Myr in order to host a star, comparable to the empirical age dispersion Luhman (2023). We assign probabilities to these grid cells p ∗ ∝ ρτ -1 ff ∝ ρ 3 / 2 - i.e. the probability of forming a star is proportional to the quantity of material multiplied by the rate at which that material is expected to collapse. For a selected grid cell, we then o ff set the location of the formed star such that it is placed randomly within the cell. We validate this choice of spatial distribution in our dynamical model a posteriori in Section 3.9.  \nFor the the total number of systems, we draw N sys that yields a total number of stars N ∗ that is broadly consistent with the findings of Luhman (2023). This is not trivial, both because there is not a clear detection limit for that sample and because we include binaries stochastically from the initial population of potential primaries, as well as the ultrawide binary fraction inferred by Joncour et al. (2017). Further, it is not straightforward to interpret what fraction of binaries are resolved / detected by Luhman (2023). We perform a two stage process as follows in Section 3.5 and 3.6.', '3.5. Ultrawide pairs': "First of all, we correct the total number of stars for a given fraction of ultrawide binaries F uwb = 2 / 3. We consider 400(1 -F uwb / 2) = 267 stars with masses > 0 . 08 M ⊙ (half of these will end up in ultrawide binaries). We then add the fraction of brown  \nFig. 3. Histogram of stellar masses in our dynamical model for Taurus. The solid black shows all stars including binary companions, while the red line is just the primaries and single stars. The black dashed line corresponds to the number that would be expected in each bin from the Kroupa (2001) IMF (equation 5). The bin sizes are similar to those adopted by Luhman (2004) in Figure 13, top panel.  \n<!-- image -->  \n∗  \n/circledot  \ndwarfs we expect below this limit from the initial mass function (IMF). For this, we use the following IMF (e.g. Kroupa 2001):  \nξ ( m ∗ ) ∝                      m -0 . 3 ∗ 0 . 01 M ⊙ ≤ m ∗ < 0 . 08 M ⊙ m -1 . 3 ∗ 0 . 08 M ⊙ ≤ m ∗ < 1 . 3 M ⊙ m -2 . 3 ∗ 0 . 5 M ⊙ ≤ m ∗ < 1 . 3 M ⊙ m -2 . 7 ∗ 1 . 3 M ⊙ ≤ m ∗ 0 otherwise , (5)  \nwith normalisation constants such that ξ is continuous and integrates to unity over all m ∗ . For each of these stars, we then draw an ultrawide companion with total probability F uwb. This yields 660 stars and sub-stellar objects in total, 38 percent of which have masses < 0 . 08 M ⊙ . We then draw the masses of all stars from the equation 5.  \nThe ultrawide pair fraction we adopt is somewhat higher than F uwb ≈ 0 . 55, inferred by Joncour et al. (2017); 186 of the 338 stars in their sample are in ultrawide pairs. However, some of these companions will separate during dynamic evolution. We also note that Joncour et al. (2017) found evidence that members of ultrawide pairs are ∼ 15 percent more likely to be themselves in a shorter period binary. However, given that this is a relatively minor correction to the binary fraction, and it is challenging to disentangle dynamical versus primordial origins for these statistics, we do not include this enhancement. While we do not forward model the initial fraction, we validate our choice by comparing the one-point correlation function in our model with that computed by Joncour et al. (2017) in Section 3.9.  \nThe separations for the ultrawide pairs is approximately loguniform between 1 . 6 × 10 3 au and 5 × 10 4 au (see Fig. 7 of Joncour et al. 2017), and we therefore draw semi-major axes of these companions similarly. We assume a uniform eccentricity distribution up to 0 . 9, and a random orientation (as described in Section 3.6).  \nWe define binaries separately from 'ultrawide pairs', described above and treated as single stars in terms of their mass-function. We consider the entire population of stars (including each member of the ultrawide pairs) to be potential primary stars in a binary system, and add companions to a subset.  \nThe binary population we include is empirically motivated, based on the findings of Moe & Di Stefano (2017). In brief, we use the same functional form for the probability of each primary having a companion in each dex of orbital period space, modified in a number of ways. Firstly, we exclude binaries with period P < 10 5 days ( ≲ 30 au); such binaries are close enough to be largely uninfluenced by encounters, and in the context of protoplanetary discs may host circumbinary discs. At longer periods, we remove the exponential taper that Moe & Di Stefano (2017) infers for P > 10 5 . 5 days, instead assuming a constant binary fraction per dex out to P = 10 7 . 7 days, approximately corresponding to the scale at which we transition to ultrawide binaries. The log-uniform distribution is consistent with the observed separations for the ultrawide pairs (see Fig. 7 of Joncour et al. 2017), and for wide binaries in the field (Lépine & Bongiorno 2007). For all binaries, we draw eccentricities uniformly from 0 -0 . 9 (e.g Abt 2006, and review by Duchêne & Kraus 2013).  \nOnce we have randomly drawn the periods for the companions of our initial population, we draw the mass-ratio and eccentricity. For the ultrawide binaries, we draw masses from the IMF, as described in Section 3.4. For the regular binaries, we draw mass-ratios q in the range 0 . 1 < q < 1 from a probability density function p ( q ). To construct p , we define two power-law regime with indices γ small for q < 0 . 3 and γ large for q ≥ 0 . 3. Following Moe & Di Stefano (2017), we also include an excess probability of the star having a twin, where a twin is defined as having a mass ratio 1 -∆ twin < q ≤ 1. This is defined in practice by computing the quantity:  \np ' twin =        F twin R 1 0 . 1 p ' ( ˜ q ) d ˜ q ∆ twin(1 -F twin) 1 -∆ twin < q < 1 0 otherwise , (6)  \nwhere  \np ' ( q ) = ( q γ small 0 . 1 ≤ q < 0 . 3 0 . 3 γ small 0 . 3 γ large q γ large q ≥ 0 . 3 . (7)  \nThen finally the probability density function for q is:  \np ( q ) = p ' + p ' twin R 1 0 . 1 p ' ( ˜ q ) + p ' twin ( ˜ q ) d ˜ q . (8)  \nSince we only consider P > 10 5 days and largely low mass stars, we have only two regimes for the power-law indices. For P < 10 6 days we have γ small = 0 . 4 and γ large = -0 . 4, and otherwise we have γ small = 0 . 5 and γ large = -1 . 1. We fix F twin = 0 . 1, for ∆ twin = 0 . 05. These values are consistent with observational constraints, although typical uncertainties are high for many of these values (El-Badry et al. 2019). When we generate our initial population, we also exclude any companions that are generated that have masses below our lower IMF limit (0 . 01 M ⊙ ), although in practice this only influences brown dwarf primaries. The overall binary fraction is ∼ 30 percent for binaries with orbital periods > 10 5 days, which is broadly consistent with the field population (e.g. Niu et al. 2021). We show the statistics of our binary and ultrawide pair distribution for the stellar population in Figure 4.  \nFinally, we generate the positions and velocities of the companion population by randomly sampling cos i (for inclination i ), the argument of periapsis, longitude of ascending node and true anomaly uniformly over the appropriate ranges. We then compute the appropriate position and velocity vector of the companion with respect to the primary. We thus produce the initial binary population that we evolve dynamically during our simulation.", '3.7. Mass function': 'We have already described our stellar mass drawing procedure for single stars and binaries in Sections 3.5 and 3.6. To validate our mass distribution including the binary population, we show the histogram of stellar masses in Figure 3. The bins used for this histogram are similar to those used by Luhman (2004, top panel of their Figure 13). The number of approximately solar mass stars found by Luhman (2004) is ∼ 40. This number is close to complete for Taurus, and is similar to our drawing from the IMF. More recently, Esplin & Luhman (2019) found a peak around m ∗ ∼ 0 . 15 and ∼ 25 stars with m ∗ ≳ 1 M ⊙ across the whole of Taurus, again comparable to our IMF draw. The mass function for all the stars, including binary companions, is shown as the black histogram in Figure 3. Despite including a di ff erent mass function for companions, the mass function is not greatly altered from the Kroupa (2001) IMF we initially assume for the primary population. We therefore conclude that we have drawn a similar stellar population to that of Taurus.', '3.8. Initial velocities': 'We must assign velocities to the primary stars we have generated, motivated by our discussion in Section 2.4. We tackle this by generating velocities following a Gaussian process, such that velocities of stars that are close to each other are more highly correlated than those of stars at large separations. If ∆ r < ∆ r max, we can write the corresponding kernel or covariance function (e.g. see equation 2.19 of Rasmussen et al. 2006):  \nk ( r , r \' ) = k ( ∆ r ) = σ 2 v , max -σ 2 v ( ∆ r ) , (9)  \nwhere σ v , max = σ v ( ∆ r max). By this definition, the covariance function is not well defined for large ∆ r > ∆ r max, where k becomes negative. However, we are free to choose an arbitrary large ∆ r max. In this case, we can also rewrite the covariance function:  \nk ( ∆ r ) ≈ σ 2 v , max GLYPH<16> 1 - D 2 r GLYPH<17> (10)  \nwhere we have defined:  \nD r = " 1 -exp -∆ r ∆ r max !# 0 . 5 . (11)  \nEquations 9 and 10 are equivalent for large ∆ r max / ∆ r → ∞ , but at large separations equation 10 defines a maximal dispersion between stellar velocities. We choose equation 10 because it yields well defined covariance for any ∆ r . In practice, we anyway choose ∆ r max = 100 pc so that our choice is not important. In order to assign velocities given this kernel function, we define the covariance matrix K = [ kij ] = [ k ( r i , r j )]. We then perform a Cholesky decomposition K = LL T , where L is a lower triangular matrix. We define a vector w α which is a vector with a length  \nFig. 4. Histogram showing the distribution of companions at a given semi-major axis (left) eccentricity (middle) and mass-ratio (right) for primary stars with masses m ∗ > 0 . 08 M ⊙ . The red histogram is for those pairs we define as ultrawide pairs, while the blue shows the binaries.  \n<!-- image -->  \nFig. 5. Data points show the proper motion di ff erence versus angular separation (assuming a distance of 140 pc) for nearest neighbours from the initial conditions in our model. The black line shows the one dimensional velocity dispersion as a function of separation and the colour bar shows the normalised Maxwell-Boltzmann distribution for each separation, as in Figure 1.  \n<!-- image -->  \ncorresponding to the number of stars for which we assign velocities. We draw each wj α ∼ N (0 , 1) independently for each spatial dimension α = 1 , 2 , 3. We then define the velocity components for the stellar population u α = Lw α .  \nIn Figure 5 we show the velocity di ff erence between di ff erent nearest neighbours in our model initial conditions, including binaries. It is clear from Figure 5 that our synthetic stellar population (scatter points) have relative velocities that follow a similar size-velocity relation, σ v ( ∆ r ) as found in Section 2.4 - cf. Figure 1.', '3.9. Model validation': "We validate our model by analysing the structure metrics discussed in Section 2. First, we consider the pairwise separation distribution, as shown in Figure 6. The pair surface density profile ˆ Σ pairs shows an excess at very small projected separations ( ∆ R ≲ 5 × 10 -4 pc), where the Luhman (2023) sample may be missing closer binaries. There is also a small excess around the  \nbinary transition at a few 10 -2 pc. This excess is quickly lost as the system evolves, and the spatial structure is in good agreement with the observed population at ∼ 1 Myr. We therefore adopt 1 Myr as the 'present day' in our simulation.  \nAt this time, the velocity structure (illustrated in Figure 7) remains similar to the velocity structure we inferred in Taurus in Section 2.4. When comparing Figure 7 to Figure 1, we note that there are some di ff erences in how they are constructed. For example, close binaries are complete in our model, but not for those in Taurus. Indeed, in Taurus the sample is restricted to only stars with Gaia proper motions, with all the biases that implies. Nonetheless, the correlation between velocity and separation remains broadly similar. The transition in nearest neighbour velocity di ff erence from 'binary' to 'field' coincides with the change in the power-law surface density profile in Figure 6 for ∆ r ∼ 3 × 10 -2 pc.  \nWe also show the one- and two-point correlation functions in Figure 8. These can be compared directly to Fig. 4 of Joncour et al. (2017). We find a similar functional form for both, again indicating that the structure in our model matches the physical structure of Taurus. We obtain a very similar region of 'inhibition' in which Ψ < 1 in the range of separations 0 . 1 -0 . 5 · . We find a similar power-law for Ψ where ∆ R < 0 . 2 · . We also recover ξ > 1 for ∆ R ≲ 1 · , close to the result of Joncour et al. (2017).  \nWe conclude that our model, with initial conditions based on physical and empirical arguments, reproduces the observed dynamical state of Taurus at 1 Myr. The age in our model of 1 Myr is somewhat younger than the 1 -3 Myr typically estimated for Taurus (e.g. Luhman 2023). If Taurus is typically older on average, this might suggest that the stellar distribution was initially more compact than the initial conditions in our simulation; although in this case we would also expect more rapid dynamical evolution. It is also possible that residual gas in the star forming region slows down the dispersal of structure. However, then it would be plausible that discs in the high density regions are replenished by accretion of this residual gas (e.g. Ku ff meier et al. 2020), blurring the lines of what ' t = 0' means for disc evolution. In truth, probably a mixture of these influences, as well as a finite period of star formation, somewhat influence the dynamical evolution of the region. Nonetheless, for the purpose of this work, we are satisfied that the current state of stars and discs in Taurus is well approximated by our simulation at 1 Myr. In this context, we note that a posteriori we find that the frequency of strong encounters does not change rapidly between 1 -3 Myr in our simulation (see Section 4.3).  \nFig. 6. The normalised pair surface density ˆ Σ pairs as a function of projected separation ∆ R in parsecs. The red line shows the observed surface density for Taurus from the Luhman (2023) census, which may be incomplete at the smallest separations. The coloured lines show the outcome from our model, at snapshot outputs indicated by the colour bar. The 1 Myr snapshot is also marked by black crosses, at which time the pair distribution function is in good agreement between observations and simulation.  \n<!-- image -->  \nFig. 7. As in Figure 5, except for our simulation after 1 Myr of evolution.  \n<!-- image -->", '3.10. Tracking encounters': "Due partly to our inclusion of binaries, interactions between stars are too complex to be easily identified with a simple criterion 'on the fly' during the simulation. We have therefore chosen to analyse encounter properties by post-processing high frequency outputs. Specifically, we have 9299 snapshots from our simulation over 3 Myr, corresponding to a time-step of 323 years. To accurately extract the correct encounter properties, we then perform the following analysis for each star i :  \n- 1. Identify the closest neighbour j ( t ) for each time-step, spatially separated by vector ∆ r i j and with velocity di ff erence ∆ v i j .\n- 2. If the nearest neighbour j to i has a bound companion k which is separated from j by vector ∆ r jk such that | ∆ r jk | <  \nFig. 8. One-point ( Ψ , blue points) and two-point ( ξ , grey points) correlation functions computed from our model at 1 Myr. The red line shows the Ψ ∝ ∆ R -1 . 5 relationship for ∆ R < 0 . 2 · and the two dashed black lines enclose the region of inhibition (where Ψ < 1), as inferred by Joncour et al. (2017, see their Figure 4).  \n<!-- image -->  \n- 0 . 1 | ∆ r i j | , we consider j and k to be a single star with the combined mass and momentum of j and k . In the following we will refer to star j as the binary barycentre.\n- 3. We search for a sign change in the vector ∆ r i j · ∆ v i j from negative (approaching) to positive (receding). We define snapshot l , at time tl , for which | ∆ r i j | is minimal between the two adjacent snapshots ( l , l + 1).\n- 4. We compute the analytic eccentricity e , closest approach distance r p, time of pericentre t p. If the predicted closest approach is not between tl and tl + 1, then this must be a nonhierarchical multiple interaction, occurring on a short timescale ( < 323 years). In this case, we assume that the closest approach is at tl (with closest approach distance given by the separation at this time), although in practice this is rare.  \nThis procedure ensures that, despite the finite temporal resolution of our simulations, we are able to resolve the majority of encounters that are relevant to disc truncation. We limit the number of encounters to one per time-step, thus we do not count every closest approach for close binaries which have an orbital period < 323 years. Although we could in principle include multiple encounters per time-step, these encounters anyway quickly truncate the disc on short time-scales. While our method may also not be accurate particularly for chaotic multiple interactions in cases where multiple interactions occur on time-scales less than 323 years, we show in Appendix A that our results are not influenced if we increase the output frequency.", '3.11. Disc truncation and initial radii': "To compute the post-encounter disc radius, we use the analytic functions inferred by Winter et al. (2018c). These functions were established by fitting a scale-free, angle-averaged expression to numerical test particle simulations, depending on the ratio of the closest approach distance r p to the outer disc radius R out, the eccentricity of the encounter e and the mass-ratio of the perturber q . Since these fitting functions were inferred for unbound encounters, we will adopt e = 1 for encounters with e < 1. We expect this to be a reasonable approximation. For example, Man-  \nara et al. (2019) fit an analytic functional form to the steady state truncation radius inferred from the numerical results of Artymowicz & Lubow (1994). For an equal mass system on a circular orbit, this estimate implies a truncation radius ∼ R out ≈ 0 . 33 r p, which is the same as the steady state truncation radius implied by our prescription (see Figure 4 of Winter et al. 2018c). In practice, the majority of encounters that strongly truncate the disc at late times also have e ∼ 1. We also consider all discs which experience encounters with r p < 10 au to be 'destroyed', making the assumption that the resultant compact disc may be rapidly accreted or photoevaporated (e.g. Clarke et al. 2001).  \nUnless otherwise stated, we initialise all discs with outer radii following (Andrews 2020):  \nR out , 0 = 250 m ∗ 1 M ⊙ ! 0 . 9 au . (12)  \nWhile this is empirically motivated by observed outer radii as inferred from CO emission lines, this relationship has a large scatter (which may in part be driven by dynamical truncation). The measurements themselves also come with numerous caveats originating from systematic uncertainties in outer radius definition, optical depth and chemistry (see e.g. Miotello et al. 2023). However, with these caveats we will here assume that the relationship is exact. We also do not consider viscous expansion (e.g. Lynden-Bell & Pringle 1974) or wind-driven contraction (e.g. Tabone et al. 2022). The discs therefore only evolve (shrink) as a result of star-disc encounters.", '3.12. Caveats for the physical model': 'There are two main caveats for the dynamical model we present here. The first is that, while we are interested in establishing the rate of strongly disc-truncating interactions at a given system age, this depends on there being a well-defined age of the system at large. More specifically, we implicitly assume that the spread in stellar ages is much smaller than the age of the system itself. This is almost certainly not the case, with the groups in Taurus having ages in the range ∼ 1 -3 Myr (Luhman 2023). However, our e ff orts are still valid in that a star-disc system will only undergo encounters with nearby neighbours, which we generally expect to have very similar ages. None of the existing examples of recent encounters have been suggested to be particularly young, but our conclusions can be reconsidered depending on the inferred ages of post-encounter systems.  \nThe other major caveat is that we do not include the selfgravity of the ISM. Without including a fully hydrodynamical model, it is not possible to include this potential in a physical way. Our choice to ignore the gas potential is in a sense a similar assumption to the one discussed above, specifically that local star formation is completed instantly, and the gas is instantaneously dispersed (e.g. by stellar feedback). If instead a large quantity of residual gas remained within the individual star forming regions that comprise the Taurus complex, then groups of stars may remain bound for longer. This may in turn increase the local encounter rate at later times. It is not possible to capture this gas potential without either full hydrodynamics (making matching exactly the structural properties of Taurus impossible), or by developing a scheme for integrating a complex potential that allows several components that co-move with stellar groups (beyond the scope of this work). However, we do find agreement between the dynamical state (i.e. correlation functions) in our model and observations, which would suggest the role of ISM self-gravity on the population dynamics as a whole cannot be dramatic. While  \nthis nonetheless remains a shortcoming of this work, we discuss in Section 4.5 that if high density gas alters our results this would have other interesting physical implications.', '4.1. Types of late-stage encounter': 'First, we qualitatively explore the nature of star-disc encounters in our simulation. For this we show some examples of encounters that occur at 1 ± 0 . 2 Myr in our simulation. The location of all these close encounters are shown in the central panel of Figure 9, from which we extract the following categories:  \n- a) Stable, long period eccentric binary: Possibly the most straightforward kind of encounter is a very long period binary which remains unperturbed by unbound stars, as shown in the top-left of Figure 9. This kind of encounter can still produce relatively large changes in outer disc radius if the orbital period is su ffi ciently long.\n- b) Chaotic multiple: Shown in the top-middle of Figure 9 is an example of a bound triple system that interacts to eventually eject one of the stars at ∼ 0 . 8 Myr, leaving a tighter binary that can further truncate the disc.\n- c) Perturbed eccentric binary: If a stable binary is perturbed by an encounter with an external (unbound) star, then this may also result in a tightening of the binary. A closer periastron distance can then further truncate the disc. We show an example of this in the top-right of Figure 9 (although in this case, the external star is in fact marginally bound).\n- d) Random unbound encounters: This is the best-studied type of encounter, the rate of which can be understood analytically at a given local stellar density and velocity dispersion. However, in our simulation, we do not find any examples of random encounters. These encounters are more common in denser regions, but in Taurus we find that the majority of stardisc interactions are mediated by wide binary companions.  \nWhile this categorisation illuminates the value of physically motivated structure and multiplicity within models for star forming regions, they are not rigid, and types of encounter may blur into each other. Our findings emphasise that unbound encounters and binary interactions cannot be studied separately during the dynamical evolution of a young star forming region.', '4.2. Disc evolution': "We consider how encounters shape the global disc radius distribution in Taurus. This is illustrated in Figure 10, where we show the cumulative distribution of disc outer radii over the course of the simulation. We find that the vast majority of discs are rapidly truncated below their initial radius by encounters, and ∼ 1 / 4 are truncated below 30 au. However, by the present time (1 Myr), dynamical encounters are not changing the overall distribution of disc radii significantly. This is because, while external perturbation does sculpt the outer disc for a large fraction of the population, disc truncation mostly occurs in stable binaries. Observations of discs in binary systems appear broadly consistent with theoretical expectations (Manara et al. 2019; Rota et al. 2022). We conclude that, for a region like Taurus, the role of encounters for the disc population as a whole is largely dominated by early interactions in binaries, rather than an ongoing process of disc truncation.  \nFig. 9. In the bottom panel, we show the spatial distribution of stars (mass > 0 . 08 M ⊙ ) at 1 Myr and contemporary disc-truncating encounters in our simulation. We show as coloured star or square symbols the location of all stars that underwent encounters yielding a fractional truncation -∆ R out / R out > 0 . 01 over the age range 1 ± 0 . 2 Myr. Encounters that yielded -∆ R out / R out > 0 . 1 are highlighted with red circles. In the panels along the top we show the separation between specific stars and their stellar neighbours over the first 1 . 5 Myr of the simulation. Each line represents the distance to a single stellar neighbour, linearly interpolated between snapshots, where neighbours that are one of the two closest stars at any given time-step are included (if they pass within 10 4 au). In each plot, we mark the location of a logged 'closest approach' as a cross (inferred analytically from the closest time-step - see Section 3.10).  \n<!-- image -->", '4.3. Truncating encounter frequency': 'We now turn to the primary motivation of this work, and ask the question: do we expect su ffi cient close encounters in Taurus to produce the examples of recent dynamical interaction? To answer this question we require two definitions:  \n- 1. Which encounters produce significant external structures (e.g. tidal tails)?\n- 2. For how long does this observable external structure persist?  \nTo answer these questions, we consider the rate of encounters that yield a fractional truncation -∆ R out / R out above some threshold, in the context of the observed flyby candidates in Taurus.', '4.3.1. HV and DO Tau': 'The huge extended dust bridge between HV and DO Tau (Howard et al. 2013) appears to be the result of a strongly truncating encounter, possibly occurring during the dynamical decay of a quadruple system (Winter et al. 2018a). While the model of Winter et al. (2018a) is probably not a unique scenario for producing the observed structure, we can make some quantitative arguments as to the requirements of such an encounter. Assuming some grain growth occurred within the disc, the mass of material expelled during the encounter is M ex ∼ 10 -4 M ⊙ (Winter et al. 2018a), which is ∼ 10 percent of the current mass M disc of  \nFig. 10. Disc outer radius R out evolution for the population of stars with m ∗ > 0 . 08 M ⊙ in our simulation. The colour bar shows the time at which the distribution is measured. The initial outer radius distribution is shown in red.  \n<!-- image -->  \nthe disc around HV Tau C (Stapelfeldt et al. 2003). If the surface density of the disc Σ ∝ R -1 , for R the cylindrical radius inside R out, this implies a fractional truncation -∆ R out / R out ∼ M ex / M disc ∼ 0 . 1. To reach the present day projected separation of ≳ 10 4 au, the encounter must have occurred ∼ 0 . 1 Myr ago. This would suggest a rate of ∼ 10 encounters per Myr for encounters that result in a fractional truncation -∆ R out / R out ≳ 0 . 1.', '4.3.2. RW Aurigae': 'For RW Aurigae, the best-fitting model explored by Dai et al. (2015) had an initial outer radius of 60 au and a final outer radius in the range ∼ 40 -57 au. This estimate is mostly based on the geometry of the spiral arm as inferred from their simplified radiative transfer. The 12 CO emission detected around RW Aurigae is at least partly optically thick, so it is not possible to reliably infer a total mass that has been ejected in the encounter. The conclusions of Dai et al. (2015) may change if the apparent morphology as seen in CO di ff ers when detailed chemistry, photodissociation or radiative transfer e ff ects are taken into account. If the physical arm is longer and wider than they assume, the authors show that a weaker encounter can be consistent with the observed structure. Dai et al. (2015) also infer a time since closest approach ∼ 600 years from their dynamical model. This conclusion may depend somewhat on the deprojected geometry of the spiral and the system orientation.', '4.3.3. UX Tau': 'In the case of UX Tau, Ménard et al. (2020) focused on two di ff erent flyby simulations, with outer radius R out = 60 au and 90 au, both with r p = 100 au, with a mass-ratio of the perturber q ∼ 0 . 08 -0 . 22. The former radius is close to the observed outer radius (post-encounter). Both simulations produce clear spiral arms and some extended structure, so a severely truncating encounter is not required. There are currently no estimates for the mass of the external structure with respect to the disc mass. From their simulations, Ménard et al. (2020) estimate ∼ 1000 years since closest approach, with some margin for uncertainty in the system geometry as in the case of RW Aurigae.', '4.3.4. Overall rates': "The considerations above lead us to conclude that a small fraction of the disc mass may be capable of producing some of the observed external structures. We therefore explore a range of fractional truncation thresholds -∆ R out / R out > 10 -3 , 10 -2 and 10 -1 . Choosing a threshold smaller than 10 -3 makes little di ff erence to the overall encounter rate. As discussed above, 'weak' encounters with -∆ R out / R out ∼ 10 -2 may be su ffi cient to reproduce the structures observed around RW Aurigae and UX Tau, while the stronger encounter threshold would correspond to systems like HV and DO Tau. While we cannot be confident of the exact -∆ R out / R out that result in RW Aurigae and UX Tau-like systems, our results will motivate future studies quantifying the mass in the external structure surrounding systems that have undergone recent encounters, given that M ex / M disc ∼ -∆ R out / R out.  \nThe overall encounter rate is summarised in Figure 11 for these thresholds. We have binned the encounters by the time at which they occur and by their fractional disc truncation. We also show the evolution of individual discs by connected lines between individual encounters, so that the evolution of the outer disc radii in (for example) binary systems is clear. As expected from Section 4.2, the majority of discs are truncated rapidly in binary systems during the first ∼ 0 . 1 Myr. However, there remain examples of individual encounters persisting throughout the course of the simulation. In order to explore the degree to which our results are stochastic, we run two additional versions of our experiment described in Appendix B.  \nQuantitatively, we can see that the rate of the strongest encounters -∆ R out / R out > 0 . 1 remains at ∼ 10 Myr -1 in the time range ∼ 1 -3 Myr (see also Appendix B). Thus HV and DO Taulike encounters, observable for ∼ 0 . 1 Myr, are likely (probability ∼ 50 percent) to be found somewhere in Taurus. We conclude that the occurrence of the HV and DO Tau encounter is expected in the context of our dynamical model.  \nThe more recent RW Aurigae- and UX Tau-like encounters require a considerably higher encounter rate. The rate of weaker encounters -∆ R out / R out > 10 -2 is Γ enc ∼ 100 Myr -1 , with approximately a factor order unity in stochastic variation (Appendix B). This rate is a factor ∼ 2 -3 larger for -∆ R out / R out > 10 -3 . We can write the probability of observing at least N encounters from a Poisson distribution:  \nP obs( N obs ≥ N ) = 1 -N -1 X i = 0 ( Γ enc τ obs) i exp( -Γ enc τ obs) i ! , (13)  \nwhere τ obs is the period for which the encounter is observable. Given that we are equally likely to observe the disc at any stage during this period of observability, this implies the average time of observation ⟨ t obs ⟩ = τ obs / 2. We therefore adopt a moderately generous τ obs = 2000 years, a factor ∼ 2 larger than the time since periastron for UX Tau and RW Aurigae. Taking a range of encounters Γ enc = 100 -200 Myr -1 for N obs ≥ 2 (RW Aurigae and UX Tau) yields P obs = 0 . 018 -0 . 062. Therefore, if encounters with -∆ R out / R out ≳ 10 -2 can produce these systems, then the tension with our model are not significant (or very marginally significant at ∼ 2 σ ). If a much more truncating encounter is required, then this tension may become significant. This marginal agreement underlines the importance of future studies quantifying the fraction of mass in the extended structure.  \nIn the absence of additional constraints, we conclude that the expected rate of encounters in Taurus is marginally su ffi cient to produce UX Tau and RW Aurigae without appealing to additional physics (see Section 4.5). However, this would not be the  \nFig. 11. The global rate of disc truncating encounters for discs around stars with mass m ∗ > 0 . 08 M ⊙ throughout the evolution of our dynamical model. The top panel shows the rate of all truncating encounters that decrease the disc radius by at least 0 . 1 percent (light gray), 1 percent (gray) and 10 percent (black) of the pre-encounter radius. The colour scale of the bottom panel shows the encounter rates binned into both the time of the encounter ( x -axis) and the degree of truncation ( y -axis). Each logged encounter is shown by a pink square, and subsequent encounters for the same disc (if any) are connected by faint pink lines.  \n<!-- image -->  \ncase if the sample of known recent, truncating star-disc encounters is incomplete. Indeed, there is some reason to suspect that this may be the case, as discussed below.", '4.4. Masses of stars undergoing close encounters': 'We can ask whether truncating encounters occur more or less frequently for high mass stars. We show the distribution of the masses of stars that undergo truncating encounters after 1 Myr in Figure 12. This mass distribution is indistinguishable to the overall mass function ( m ∗ > 0 . 08 M ⊙ ). By contrast, RW Aurigae A and B have masses ∼ 1 . 3 -1 . 4 M ⊙ and ∼ 0 . 7 -0 . 9 M ⊙ respectively (Ghez et al. 1997; Woitas et al. 2001). HV Tau C has a mass ∼ 0 . 5 -1 M ⊙ (Duchêne et al. 2010) and DO Tau ∼ 0 . 3 -0 . 7 M ⊙ (Beckwith et al. 1990; Hartigan et al. 1995). UX Tau A has a mass ∼ 1 . 3 M ⊙ and UX Tau C has mass ∼ 0 . 16 M ⊙ (Kraus & Hillenbrand 2009; Zapata et al. 2020). The origin of the extended structures are the discs around RW Aurigae A, HV Tau C, DO Tau and UX Tau A. The masses of these stars appear to be systematically greater than what would be obtained from randomly sampling from the IMF. Assuming nominal masses of 1 . 4 M ⊙ , 0 . 8 M ⊙ , 0 . 5 M ⊙ and 1 . 3 M ⊙ for RW Aurigae A, HV Tau C, DO Tau UX Tau A respectively yields KS test p -value  \n0 . 023 percent when comparing to the stars that undergo encounters in our model. It is quite possible that observations could be biased to detect evidence of encounters involving more massive stars with brighter discs. However, this would imply many more undetected external structures that are evidence of encounters involving low mass stars in Taurus. Following equation 13, if evidence for just one more encounter similar to RW Aurigae or UX Tau were uncovered in Taurus, this would being the 2 σ encounter rate to Γ enc ∼ 680 Myr -1 . Alternatively, this would yield P obs ≈ 10 -3 even for weak encounters with Γ enc ∼ 100 Myr -1 . Finding further examples would introduce significant tension between the frequency of such events and the encounter rates in our model (Section 4.3.4), necessitating additional physics. This highlights the importance of future unbiased surveys to search for evidence of recent star-disc encounters in Taurus.  \nWe caveat our findings with the fact that we do not invoke any primordial mass segregation (Zinnecker et al. 1993; Moeckel & Bate 2010; Plunkett et al. 2018). However, mass segregation is not evident in Taurus (Dib & Henning 2019).  \n/circledot  \nFig. 12. Cumulative distribution function for the masses m ∗ of stars with m ∗ > 0 . 08 M ⊙ in our simulation (black line), compared with those that underwent a truncating encounter -∆ R out / R out > 0 . 01 (red line). The distributions are not significantly di ff erent, with KS test probability p KS = 0 . 87. We also show mass estimates for the stars with discs that are responsible for the observed external structure. We do not show uncertainties in these estimates for clarity, but errors quoted are typically ∼ ± 0 . 3 M ⊙ . The distribution of masses of stars that have been inferred to have experienced recent truncating encounters are significantly different from those in our simulation, with p KS = 0 . 023.  \n<!-- image -->', '4.5. Mechanisms to increase late encounter rates': 'Given the apparently high probability that the observed sample of discs with external structure produced during a recent stardisc encounter is incomplete, we consider a number of possibilities that may enhance the frequency of encounters late-on during the dynamical evolution of Taurus.', '4.5.1. (Viscous) re-expansion': 'For discs that experience multiple encounters that are either unbound or in long period binaries, viscous angular momentum transport (e.g. Lynden-Bell & Pringle 1974) or in-fall of material from the ISM (e.g. Padoan et al. 2005; Manara et al. 2018; Ku ff meier et al. 2023; Gupta et al. 2023; Padoan et al. 2024; Winter et al. 2024) may replenish outer disc material such that subsequent encounters will have a stronger influence than under our assumption that disc radius is fixed. This replenishment should be correlated with stellar mass (as suggested by observed stellar accretion rates - Manara et al. 2017), so this would also preferentially enhance encounter rates for high mass stars, consistent with observations (Section 4.4). Conversely, if angular momentum is extracted by magnetohydrodynamic winds (e.g. Bai & Stone 2013) then this could exacerbate the encounter rate problem.', '4.5.2. Self-gravity of the interstellar medium': 'As discussed in Section 3.12, for practical reasons we have ignored ISM self-gravity in our calculations. If included, we would expect this self-gravity to increase the time for which groups of stars remain bound and therefore potentially enhance the frequency of truncating encounters at late times. This in-of-itself  \nFig. 13. Critical density at which the ISM is bound against turbulent pressure (black line, left hand side y -axis) and the resultant BondiHoyle-Lyttleton (BHL) accretion rate for a solar mass star (red dashed line, right hand side y -axis). Both are shown as a function of spatial scale λ .  \n<!-- image -->  \n/circledot  \nmay explain any shortcoming in the frequency of encounters in our model compared to observations. However, for this to keep stars bound together, gas must be at approximately the critical density to undergo gravitational collapse. On length scales λ relevant for the Taurus region ( σ v ≫ c s and λ ≪ h , the galactic scale height), this is e ff ectively the Jeans criterion: the gravitational potential must balance the turbulent energy. At such densities, as discussed above, a number of authors have suggested that disc replenishment by Bondi-Hoyle-Lyttleton (BHL) accretion may be substantial (Padoan et al. 2005; Throop & Bally 2008; Ku ff meier et al. 2023; Winter et al. 2024).  \nTo test whether we expect stellar dynamics-altering ISM self-gravity to also substantially replenish protoplanetary discs, we derive the critical density ρ c as a function of λ following Winter et al. (2024), shown as the black line in Figure 13. We then estimate the expected typical BHL accretion rate:  \n˙ M BHL ∼ 4 π G 2 m 2 ∗ ρ c σ 3 v (14)  \nfor a solar mass star m ∗ = 1 M ⊙ , shown as the red dashed line in Figure 13. Here σ v ∝ λ 0 . 5 as we assume when implementing kinematic substructure. For typical regions of size λ ∼ 1 pc (Schmalzl et al. 2010; Joncour et al. 2018), our calculations imply a BHL accretion rate ˙ M BHL ∼ 10 -8 M ⊙ yr -1 , comparable to observed stellar accretion rates (e.g. Manara et al. 2023). Therefore, if local dynamics is influenced by ISM self-gravity, we expect this also to replenish protoplanetary discs as discussed above in Section 4.5.1.  \nWe summarise that, while we cannot rule out ISM selfgravity as a possible driver of present day encounters in Taurus, this would have extremely interesting consequences for disc evolution. In particular, it would imply that replenishment from the ISM is indeed common.', '4.5.3. Tidal torques from the disc on the perturber': 'If the disc is su ffi ciently massive and the encounter su ffi ciently close, then it is possible that the tidal dissipation of the orbital energy of the perturber leads to capture and tightening of the bound orbit (Clarke & Pringle 1993; Ostriker 1994; Muñoz et al. 2015). In this case the subsequent encounters may be stronger as  \nthe periastron separation decreases. This is a promising way to generate delayed strong encounters in binary / multiple systems where multiple encounters have occurred, particularly given the evidence for several close passages in the RW Aurigae system (Rodriguez et al. 2018). Since the disc mass is a superlinear function of star mass (e.g. Ansdell et al. 2016), this mechanism could also be more e ff ective for perturbations to discs around massive stars. A key question however, is whether discs retain su ffi cient mass to generate orbital decay to the present day; this may be the subject of future hydrodynamic simulations for RW Aurigae and UX Tau.', '4.5.4. Embedded planets or brown dwarfs': 'If brown dwarfs or massive planets form and remain embedded in the protoplanetary discs, then it is possible that during the encounter they undergo large eccentricity perturbations (Heggie & Rasio 1996) or chaotic evolution, perhaps to be captured by the perturbing star (e.g. Fregeau et al. 2004). This would e ff ectively be a late-stage enhancement of the binary fraction, which may be expected to enhance the encounter rate. As for the previous two mechanisms, massive planets appear to be more common around high mass stars (Johnson et al. 2010). While an interesting possibility in the context of spiral arms in discs, given the presence of a stellar-mass perturber in the cases of RW Aurigae and UX Tau this does not appear to be a convincing explanation.', '4.5.5. Spatial variation of the binary fraction': 'In our simulations, we initiated binaries independently of their location. It is plausible that binary formation is more probable in regions of enhanced local density, where enhanced tides can in principle result in fragmentation into multiple systems (Horton et al. 2001). This would enhance the degree to which binaries typically interact with other stars, and may therefore lead to enhanced encounter rates. However, RW Aurigae at least does not appear to be located in a region of high stellar density (e.g. Pfalzner & Govind 2021).', '4.5.6. Prospects': 'While these explanations discussed in this section are not yet necessary based on existing constraints, two future developments may change this conclusion. First, the quantity of mass in the extended structure surrounding RW Aurigae and UX Tau could be a substantial fraction of the disc mass (that may imply | ∆ R out / R out | ≫ 10 -2 ). Second, any new discoveries of extended structure in Taurus originating from a star-disc encounter will considerably increase the required encounter rate to Γ enc ∼ 680 Myr -1 , much greater than the rate for even weak encounters in our model Γ enc ∼ 200 Myr -1 . Should an additional recent encounter be discovered, and given the arguments outlined above, possibly the most promising mechanism for enhancing the encounter rate is the (viscous) re-expansion or ISM replenishment of the disc. This hypothesis should be quantitatively investigated in the event of additional discoveries.', '5. Conclusions': "We have presented an N-body model for the Taurus star forming region which is physically and empirically motivated. Our initial stellar population is generated probabilistically from a turbulent gaseous medium, via a zoom-in on a Gaussian field us-  \ning the cosmology simulation tool genet IC (Stopyra et al. 2020). We applied the inferred size-velocity relation in Taurus to generate an initial velocity dispersion, combined with an empirically motivated binary population. Without any additional tuning, this yields a dynamical model with excellent agreement to the observed structure in Taurus by the metrics of the pair separation distribution and separation-velocity correlation. This has allowed us, for the first time, to accurately assess the global frequency of encounters in the Taurus star forming region.  \nIn this way we have shown that, like the bull that is its namesake, stars in the Taurus star forming region rarely settle for one close encounter. Instead, star-disc systems act as stellar matadors, often enduring several close approaches with a neighbouring star. The closest approach distances change stochastically over time as binaries are themselves perturbed by low velocity neighbours. High order multiples can also form and break up on time-scales much shorter than the lifetime of Taurus. As a result, strong encounters at the present day can occur in one of four ways. The most common are: during the evolution of a chaotic multiple system; when an eccentric binary is perturbed; or as a result of a close approach in a very long period eccentric binary. Random encounters between single stars are rare, with close encounters in Taurus being mostly mediated through a binary companion. These categories can be blurred, and in some sense to distinguish between binary and unbound encounters is a false dichotomy. Since binaries and nearby low velocity stars influence each other, inferring accurate encounter rates requires that the spatial and kinematic structure of a region is quantitatively taken into account.  \nOverall, ∼ 1 / 4 of discs are truncated below 30 au by dynamical encounters. However, the majority of these dynamical truncation events happen in the first few 0 . 1 Myr of the cluster evolution, over the course of a few binary periods. After this time, the role of encounters in sculpting the overall distribution of disc radii is largely finished in a low density star forming region such as Taurus.  \nNonetheless, individual strong encounters still occur over the region as a whole, and we consider whether the examples of HV and DO Tau (Howard et al. 2013; Winter et al. 2018a), RW Aurigae (Dai et al. 2015; Rodriguez et al. 2018) and UX Tau (Zapata et al. 2020; Ménard et al. 2020) should be observed in our model. We conclude that events resembling HV and DO Tau can occur at a rate of ∼ 10 Myr -1 and therefore, given that it has remained observable for ∼ 0 . 1 Myr, observing one such event at any given time is expected. If weak encounters that eject only ∼ 0 . 1 -1 percent of the disc mass are responsible for the external structure around RW Aurigae and UX Tau, then they are also expected given the encounter rate Γ enc ∼ 100 -200 Myr -1 in our model. However, we highlight that the systems that have been inferred to have experienced recent truncating encounters do not appear to be consistent with random drawing from the mass function, while those in our model are. This hints that the known sample of discs that have recently experienced truncating encounters is incomplete: just one more observational example would be in strong tension with our model, implying additional physics.  \nWe discuss a number of physical mechanisms that should be explored by future work that may yield enhancements in the rate of truncating encounters, such as the tightening of binaries due to star-disc interaction (e.g. Muñoz et al. 2015) or reexpansion / replenishment in the outer regions of protoplanetary discs. We also show that substantial replenishment via BHL accretion must be proceeding if ISM self-gravity, which we neglect  \nin this work, is significant for the dynamical evolution of high order multiples across the Taurus complex.  \nWe summarise that star-disc encounters are an important probe of disc physics. This work highlights the need for a systematic search for extended structure generated by star-disc encounters in Taurus.  \nAcknowledgements. We thank the anonymous referees for their useful comments that helped improve the clarity of this manuscript. This project has received funding from the European Research Council (ERC) under the European Union's Horizon research and innovation programme (grant agreement No. 101002188, project PROTOPLANETS, and grant agreement No. 101042275, project Stellar-MADE). AJW has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 101104656. LS was supported by NSFC grant Nos. 11890692, 12133008, and 12221003. LS thanks T Fang for support and acknowledges science research grants from NSFC, grant Nos. 11890692, 12133008, and 12221003. This work has made use of data from the European Space Agency (ESA) mission Gaia ( https://www.cosmos.esa.int/gaia ), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium ). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.", 'References': '- Aarseth, S. J. 2003, Gravitational N-Body Simulations Abt, H. A. 2006, ApJ, 651, 1151 Adams, F. C., Proszkow, E. M., Fatuzzo, M., & Myers, P. C. 2006, ApJ, 641, 504 Alaguero, A., Cuello, N., Ménard, F., et al. 2024, V892 Tau: A tidally perturbed circumbinary disc in a triple stellar system Andrews, S. M. 2020, ARA&A, 58, 483 Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46 Artymowicz, P. & Lubow, S. H. 1994, ApJ, 421, 651 Babusiaux, C., Fabricius, C., Khanna, S., et al. 2023, A&A, 674, A32 Baehr, H. & Zhu, Z. 2021, ApJ, 909, 135 Bai, X.-N. & Stone, J. M. 2013, ApJ, 769, 76 Bastian, N., Lamers, H. J. G. L. M., de Mink, S. E., et al. 2013, MNRAS, 436, 2398 Bate, M. R. 2018, MNRAS, 475, 5618 Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924 Brunt, C. M. 2010, A&A, 513, A67 Burgers, J. 1948, in Advances in Applied Mechanics, Vol. 1, A Mathematical Model Illustrating the Theory of Turbulence, ed. R. Von Mises & T. Von Kármán (Elsevier), 171-199 Cabrit, S., Pety, J., Pesenti, N., & Dougados, C. 2006, A&A, 452, 897 Cadiou, C., Pontzen, A., Peiris, H. V., & Lucie-Smith, L. 2021, MNRAS, 508, 1189 Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 Clarke, C. J. & Pringle, J. E. 1993, MNRAS, 261, 190 Craig, J. & Krumholz, M. R. 2013, ApJ, 769, 150 Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114 Cuello, N., Louvet, F., Mentiplay, D., et al. 2020, MNRAS, 491, 504 Cuello, N., Ménard, F., & Price, D. J. 2023, European Physical Journal Plus, 138, 11 Dai, F., Facchini, S., Clarke, C. J., & Haworth, T. J. 2015, MNRAS, 449, 1996 De Rosa, R. J. & Kalas, P. 2019, AJ, 157, 125 Dib, S. & Henning, T. 2019, A&A, 629, A135 Dong, R., Liu, H. B., Cuello, N., et al. 2022, Nature Astronomy, 6, 331 Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 Douglas, T. A., Caselli, P., Ilee, J. D., et al. 2013, MNRAS, 433, 2064 Duchêne, G. & Kraus, A. 2013, ARA&A, 51, 269 Duchêne, G., McCabe, C., Pinte, C., et al. 2010, ApJ, 712, 112 Dullemond, C. P., Kü ff meier, M., Goicovic, F., et al. 2019, A&A, 628, A20 El-Badry, K., Rix, H.-W., Tian, H., Duchêne, G., & Moe, M. 2019, MNRAS, 489, 5822 Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 211 Esplin, T. L. & Luhman, K. L. 2019, AJ, 158, 54 Forgan, D. & Rice, K. 2010, MNRAS, 402, 1349 Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 2004, MNRAS, 352, 1 Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1 Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1 Ghez, A. M., White, R. J., & Simon, M. 1997, ApJ, 490, 353 Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428 Gupta, A., Miotello, A., Manara, C. F., et al. 2023, A&A, 670, L8 Ha, T., Li, Y., Kounkel, M., et al. 2022, ApJ, 934, 7\n- Haisch, Jr., K. E., Lada, E. A., Piña, R. K., Telesco, C. M., & Lada, C. J. 2001, AJ, 121, 1512\n- Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736\n- Heggie, D. C. & Rasio, F. A. 1996, MNRAS, 282, 1064\n- Heyer, M. H. & Brunt, C. M. 2004, ApJ, 615, L45\n- Holmberg, J., Nordström, B., & Andersen, J. 2009, A&A, 501, 941\n- Horton, A. J., Bate, M. R., & Bonnell, I. A. 2001, MNRAS, 321, 585\n- Howard, C. D., Sandell, G., Vacca, W. D., et al. 2013, ApJ, 776, 21\n- Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905\n- Joncour, I., Duchêne, G., & Moraux, E. 2017, A&A, 599, A14\n- Joncour, I., Duchêne, G., Moraux, E., & Motte, F. 2018, A&A, 620, A27\n- Kraus, A. L. & Hillenbrand, L. A. 2009, ApJ, 704, 531\n- Kroupa, P. 2001, MNRAS, 322, 231\n- Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227 Ku ff meier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161\n- Ku ff meier, M., Goicovic, F. G., & Dullemond, C. P. 2020, A&A, 633, A3 Ku ff meier, M., Jensen, S. S., & Haugbølle, T. 2023, European Physical Journal Plus, 138, 272\n- Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57  \nLarson, R. B. 1981, MNRAS, 194, 809  \n- Lépine, S. & Bongiorno, B. 2007, AJ, 133, 889\n- Li, D., Mustill, A. J., Davies, M. B., & Gong, Y.-X. 2024, MNRAS, 527, 386\n- Lichtenberg, T., Parker, R. J., & Meyer, M. R. 2016, MNRAS, 462, 3979\n- Luhman, K. L. 2004, ApJ, 617, 1216\n- Luhman, K. L. 2023, AJ, 165, 37\n- Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603\n- Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 539\n- Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3\n- Manara, C. F., Tazzari, M., Long, F., et al. 2019, A&A, 628, A95\n- Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127\n- Ménard, F., Cuello, N., Ginski, C., et al. 2020, A&A, 639, L1\n- Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24\n- Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 501\n- Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15\n- Moeckel, N. & Bate, M. R. 2010, MNRAS, 404, 721\n- Muñoz, D. J., Kratter, K., Vogelsberger, M., Hernquist, L., & Springel, V. 2015, MNRAS, 446, 2010\n- Münch, G. 1958, Reviews of Modern Physics, 30, 1035\n- Niu, Z., Yuan, H., Wang, S., & Liu, J. 2021, ApJ, 922, 211\n- Nordlund, Å. K. & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco & A. Carraminana, 218\n- Ostriker, E. C. 1994, ApJ, 424, 292\n- Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 513, 259\n- Padoan, P., Jones, B. J. T., & Nordlund, Å. P. 1997, ApJ, 474, 730\n- Padoan, P., Kritsuk, A., Norman, M. L., & Nordlund, Å. 2005, ApJ, 622, L61\n- Padoan, P., Pan, L., Pelkonen, V.-M., Haugboelle, T., & Nordlund, A. 2024, arXiv e-prints, arXiv:2405.07334\n- Parker, R. J. 2023, MNRAS[ arXiv:2308.05790 ]\n- Parker, R. J. & Schoettler, C. 2023, ApJ, 952, L16\n- Pearson, S. G. & McCaughrean, M. J. 2023, arXiv e-prints, arXiv:2310.01231\n- Pfalzner, S. 2008, A&A, 492, 735\n- Pfalzner, S. & Govind, A. 2021, ApJ, 921, 90\n- Pfalzner, S., Vogel, P., Scharwächter, J., & Olczak, C. 2005, A&A, 437, 967\n- Plunkett, A. L., Fernández-López, M., Arce, H. G., et al. 2018, A&A, 615, A9\n- Portegies Zwart, S. & Hochart, E. 2023, arXiv e-prints, arXiv:2312.04645\n- Qian, L., Li, D., Gao, Y., Xu, H., & Pan, Z. 2018, ApJ, 864, 116\n- Rasmussen, C. E., Williams, C. K., et al. 2006, Gaussian processes for machine learning, Vol. 1 (Springer)\n- Ren, B., Dong, R., van Holstein, R. G., et al. 2020, ApJ, 898, L38\n- Rodriguez, J. E., Loomis, R., Cabrit, S., et al. 2018, ApJ, 859, 150\n- Rota, A. A., Manara, C. F., Miotello, A., et al. 2022, A&A, 662, A121\n- Schmalzl, M., Kainulainen, J., Quanz, S. P., et al. 2010, ApJ, 725, 1327\n- Shara, M. M., Hurley, J. R., & Mardling, R. A. 2016, ApJ, 816, 59\n- Shuai, L., Ren, B. B., Dong, R., et al. 2022, ApJS, 263, 31\n- Spurzem, R. 1999, Journal of Computational and Applied Mathematics, 109, 407  \nSpurzem, R., Giersz, M., Heggie, D. C., & Lin, D. N. C. 2009, ApJ, 697, 458  \nStapelfeldt, K. R., Ménard, F., Watson, A. M., et al. 2003, ApJ, 589, 410  \nStopyra, S., Pontzen, A., Peiris, H., Roth, N., & Rey, M. 2020, GenetIC: Initial  \n- conditions generator for cosmological simulations, Astrophysics Source Code Library, record ascl:2006.020\n- Stopyra, S., Pontzen, A., Peiris, H., Roth, N., & Rey, M. P. 2021, ApJS, 252, 28  \nArticle number, page 15 of 20  \nTabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290  \n- Throop, H. B. & Bally, J. 2008, AJ, 135, 2380\n- Vazquez-Semadeni, E. 1994, ApJ, 423, 681\n- Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261 Vorobyov, E. I., Elbakyan, V. G., Liu, H. B., & Takami, M. 2021, A&A, 647, A44\n- Vorobyov, E. I., Steinrueck, M. E., Elbakyan, V., & Guedel, M. 2017, A&A, 608, A107\n- Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070\n- Wang, Y., Perna, R., & Zhu, Z. 2023, arXiv e-prints, arXiv:2310.06016\n- Winter, A. J., Benisty, M., & Andrews, S. M. 2024, arXiv e-prints, arXiv:2405.08451\n- Winter, A. J., Booth, R. A., & Clarke, C. J. 2018a, MNRAS, 479, 5522\n- Winter, A. J., Clarke, C. J., Rosotti, G., & Booth, R. A. 2018b, MNRAS, 475, 2314\n- Winter, A. J., Clarke, C. J., Rosotti, G., & Giersz, M. 2022, MNRAS, 515, 2837\n- Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018c, MNRAS, 478, 2700\n- Winter, A. J. & Haworth, T. J. 2022, European Physical Journal Plus, 137, 1132 Woitas, J., Leinert, C., & Köhler, R. 2001, A&A, 376, 982\n- Zapata, L. A., Rodríguez, L. F., Fernández-López, M., et al. 2020, ApJ, 896, 132 Zinnecker, H., McCaughrean, M. J., & Wilking, B. A. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine, 429', 'Appendix A: Output frequency exploration': "Our approach to computing the encounter frequency has been to post-process high frequency outputs, following the method described in Section 3.10. For our fiducial model, we adopt an output frequency of one every 323 years. While this frequency is su ffi cient to capture even binary encounters if the semi-major axis a ≳ 50 au, it is possible that in some instances binary-single or binary-binary interactions result in miscalculation of the true encounter properties. Here we test if such events may change our results.  \nWhile it would be laborious (and challenging) to check every single star's encounter history in our sample for examples where our approach does not capture close encounters, we can more directly test our choice on the disc truncation history. We therefore decrease the output time-step by factor ten (one every 32 years), and repeat the same analysis on the rate of disc truncation as before. The results are summarised in Figure A.1, which can be compared directly to Figure 11. Statistically and qualitatively, our results are very similar. We therefore conclude that our approach for extracting encounters is not substantially altering our conclusions. However, we do notice some small di ff erences in the encounter history that correspond to a small number of individual encounter histories. These di ff erences are not important for our conclusions, but we investigate them as follows.  \nWe first identify an example of a system for which we observe di ff erences in the encounter history depending on the output frequency. We show one such encounter history in Figure A.2 for our fiducial model, which is a chaotic triple interaction. The close encounters identified by our algorithm are shown as crosses. Comparison with the high frequency output simulation (Figure A.3) shows that at early times, both encounter histories are identical. They diverge after ∼ 0 . 3 Myr, from which point the dynamics of the systems evolve chaotically to di ff erent endstates. This suggests that the di ff erences in the encounter histories inferred from Figure A.1 compared to Figure 11 are not due to di ff erences in our encounter-extraction algorithm, but di ff erences in the N-body integration itself. While in principle changing the frequency of outputs should not alter the integration, N body 6 ++ performs a number of accountancy operations at the time of output. It is possible that these operations slightly alter other numbers in the code that enter into numerical calculations via, for example, the adjustment of parameters. We demonstrate that altering the parameter adjustment time-step ( DTADJ ) can have a similar influence on the dynamical evolution of chaotic multiples in Figure A.4, where we reduce this time-step by five compared to our fiducial model.  \nWe do not here investigate what changes to the output timestep lead to an altered chaotic evolution of high order multiples when using the N body 6 ++ code. By the nature of chaotic interactions, such changes may be tiny (such as machine precision) and in this case no particular solution is obviously more accurate. This is particularly irrelevant physically, since for these cases other processes may also change the dynamical outcome. However, we are satisfied that the statistical distribution of star-disc encounters our model is not dependent on our choice of output time-step, and that our algorithm for extracting encounters is adequate for our purposes.", 'Appendix B: Stochastic encounter history': 'To ensure that we are not dominated in our quantitative conclusions by stochastic variations in the encounter rates, we perform two additional numerical integrations of a Taurus-like re-  \ngion. We do this by performing two resampling experiments for the same density field as we adopt for our fiducial model, then adding a new binary population. The outcome comparing the encounter rates across all the simulations is shown in Figure B.1. Figures comparable to Figure 11 are shown in Figures B.2 and B.3. While we are not able to repeat the experiment enough times to gain a full distribution at each time interval, we can estimate a factor ∼ 2 variation in the encounter rates is typical. We conclude that the rate of truncating encounters that we predict is only stochastic to within a factor of order unity.  \nFig. A.1. As in Figure 11, but for a factor ten higher output frequency.  \n<!-- image -->  \nFig. A.2. An example of a chaotic multiple interaction in our fiducial model. We show the separation of two neighbouring stars from a third star as a function of time. Crosses mark the location where the closest encounter is recorded by our post-processing analysis.  \n<!-- image -->  \nFig. A.3. As in Figure A.2, for the same stars but with a factor ten higher output frequency.  \n<!-- image -->  \nFig. A.4. As in Figure A.2, but with a factor five smaller time-step between the adjustment of simulation paramaters ( DTADJ parameter in N body 6 ++ ).  \n<!-- image -->  \nFig. B.1. Stochastic variation of the encounter rate between di ff erent simulations (shown by di ff erent line styles). The colour of the lines refers to the threshold truncation used to define the encounter (as in Figure 11).  \n<!-- image -->  \nFig. B.2. As in Figure 11 but for a second random drawing of the stellar and binary population.  \n<!-- image -->  \nFig. B.3. As in Figure B.2 but for a third random drawing of the stellar and binary population.  \n<!-- image -->'}