- Split View
-
Views
-
Cite
Cite
Fei Dai, Stefano Facchini, Cathie J. Clarke, Thomas J. Haworth, A tidal encounter caught in the act: modelling a star–disc fly-by in the young RW Aurigae system, Monthly Notices of the Royal Astronomical Society, Volume 449, Issue 2, 11 May 2015, Pages 1996–2009, https://doi.org/10.1093/mnras/stv403
- Share Icon Share
Abstract
RW Aurigae (RW Aur) is a binary star system with a long molecular arm trailing the primary star. Cabrit et al. noted the resemblance between this extended structure and the tidal arm stripped from the primary star in the simulations of star–disc encounters by Clarke & Pringle. In this paper, we use new hydrodynamical models and synthetic observations to fit many of the parameters of RW Aur. Using hydrodynamic models, we find that the morphological appearance of RW Aur can be indeed explained by a tidal encounter with the secondary star. We reproduce all the major morphological and kinematic features of the system. Using radiative transfer calculations, we find that synthetic CO and dust continuum observations of our hydrodynamic models agree well with observations. We reproduce all the main features of the line profiles, from emission fluxes to the optical depth of the different components of the system. The agreement between observations and simulations thus lends strong support to the hypothesis of a tidal encounter scenario. Finally, we propose a possible solution for the origin of the dimming of the primary star observed in 2010/2011 by Rodriguez et al.
1 INTRODUCTION
RW Aurigae (RW Aur) is a young stellar system comprising two roughly solar mass stars at a projected separation of 1.5 arcsec, i.e. 200 au at a distance d ≃ 140 pc (cf. Hipparcos catalogue; van Leeuwen 2007). The primary is a Classical T Tauri Star which is characterized by an unusually high ratio of accretion rate (2–10 × 10−7 M⊙ yr−1; Basri & Bertout 1989; Hartigan, Edwards & Ghandour 1995) to disc mass (∼10−4 M⊙ as derived from mm continuum dust emission; Osterloh & Beckwith 1995).
The most notable feature of the RW Aurigae system is the 600-au-long-arm structure trailing from the disc of RW Aur A which was imaged by Cabrit et al. (2006) in 12CO(J = 2–1) and 12CO(J = 1–0) with the Institut de Radioastronomie Millimétrique (IRAM) Plateau de Bure Interferometer. Cabrit et al. (2006) noted the resemblance between this structure and the tidal arm stripped from the primary star in the simulations of star–disc encounters by Clarke & Pringle (1993). Although these simulations were not designed to fit a particular observation, the superficial resemblance between the morphology of RW Aur and the prograde encounter of Clarke & Pringle (1993) was striking. It is however unclear without detailed modelling whether simulations can fit all the parameters of the system, particularly the gas kinematics and the available radial velocity and proper motion data for the two stars. Moreover, Bisikalo et al. (2012) presented evidence from proper motion data that the encounter might in fact be retrograde in which case the Clarke & Pringle (1993) models would predict a morphology unlike that observed in RW Aur. Bisikalo et al. (2012) have on this basis proposed an alternative model involving circumbinary disc accretion. In the detailed modelling that we present below, we shall attempt to fit all the available kinematic constraints through a star–disc encounter (i.e. we assume that the system is either highly eccentric or marginally unbound. For justification of this assumption, see Section 3.3), but note that, since the errors on the proper motion data in Bisikalo et al. (2012) are large, we shall place more emphasis in fitting the radial component of the stars’ relative space motions on the plane of the sky.
Interest in the RW Aurigae system has recently been revived by the photometric observations of Rodriguez et al. (2013). Their optical monitoring detected a 2 mag dimming event lasting 180 d which the authors intriguingly interpreted in terms of occultation of RW Aur A by a portion of the tidal arm. We shall investigate whether our best model for the encounter involves a geometry that would permit such an event.
In this paper, we use new hydrodynamical models and synthetic observations to fit many of the parameters of RW Aur. We aim to understand the geometry of the binary system, the mechanism/geometry responsible for its morphological appearance, other observed properties (such as molecular line profiles) and the recent dimming event found by Rodriguez et al. (2013). The rest of this paper proceeds as follows. In Section 2, we describe the range of observational constraints of the RW Aur system which we shall use to constrain our models. In Section 3, we outline the codes used in our investigation (the phantom1 smooth particle hydrodynamics code which is used for the hydrodynamical simulations and the torus2 Monte Carlo radiation transport code which is used to produce synthetic observations from the results of the hydrodynamic models) and also the initial conditions and parameter space surveyed. In Section 4, we present and discuss our best-fitting model solution and the degree to which various system parameters are constrained by the observations. In Section 5, we present the detailed line profiles and continuum emission map for our best-fitting model. Section 6 discusses a possible origin of the ‘dimming event’ in the light curve of RW Aur, and whether it could be associated with occultation by the tidal arm. Section 7 summarizes our conclusions.
2 OBSERVATIONS OF RW Aur
2.1 The disc around star A
The 12CO(J = 2–1) emission from close proximity to RW Aur A shows a perfect point symmetry between its redshifted and blueshifted lobes. Furthermore, the radial velocity gradient is perpendicular to the optical jet. Both of these features suggest an axisymmetric geometry, i.e. RW Aur A most likely has a rotating disc (Cabrit et al. 2006). The sky inclination of the disc was constrained using the ratio of the radial and tangential velocities of emission knots in its jet (López-Martín, Cabrit & Dougados 2003). Alternatively, the inclination was also constrained by the radial velocity of the jet itself. Combining these, the sky inclination was estimated to be in the range of 45–60° (Woitas et al. 2002).
The observed flux densities in 12CO(J = 2–1) and 12CO(J = 1–0) has a ratio of about 4 (see fig. 4 by Cabrit et al. 2006). The ratio of 4 is the square of the frequency ratio between these two lines. This suggests that the disc is optically thick in which case the Rayleigh–Jeans limit predicts that Bν = 2kTν2/c2. Beckwith & Sargent (1993, hereafter BS93) predicted that for an optically thick, Keplerian, vertically isothermal disc of radial temperature profile T ∝ R−q, the molecular line profile scales as ν3q − 5 towards the wings of the line profile. The disc line profile observed by Cabrit et al. (2006) suggests that q ∼ 0.73. The same parameter q was also independently estimated as ∼0.57 from the infrared spectral energy distribution (Osterloh & Beckwith 1995). The BS93 model also predicted that the peak flux density of disc line profiles occurs at the velocity of the disc's outer radius: Vpeak = sin (i)(GM/Rout)0.5. The observed peak for A's disc (∼3.5 km s−1) then constrains Rout to be in the range of 50–77 au, respectively, for an inclination of 45–60°. Cabrit et al. (2006) noted that stellar mass and outer disc radius are degenerated when considering the observed peak velocity. Therefore, they argued that the better approach is to compare the shape (rather than absolute peak velocity) of observed and theoretical line profiles in BS93; they claimed Rout to be 40–57 au.
Cabrit et al. (2006) estimated the temperature of the disc at its outer radius with two methods.
Measuring the peak flux density of the disc line profile and using the BS93 model which predicted that the peak flux density scales as |$F(\lambda ) \sim \lambda ^{-2} R_{\rm out} ^{5/2} T_{\rm out}^ {3/2} M_{\rm star} ^{-0.5} d^{-2}$|.
Using the brightness temperature since the disc is optically thick and in the Rayleigh–Jeans regime.
The two methods consistently gave a temperature of 60–107 K (the range corresponding to different sky inclination) at the outer radius of the disc.
Cabrit et al. (2006) further estimated the column density of the disc. The fact that the disc is optically thick in 12CO and undetected in 13CO placed lower and upper limits on the column density, respectively. These CO column density limits were converted to total column densities using standard 12CO and 13CO abundances in the local interstellar medium. Assuming that the disc surface density scales as 1/R, Cabrit et al. (2006) derived an estimate of the total mass of the disc of 10−5–10−4 M⊙. This is consistent with the mass estimate by Osterloh & Beckwith (1995, 3 × 10−4 M⊙), but not with the higher estimate by Andrews & Williams (2005, 4 × 10−3 M⊙).
2.2 The tidal arm
12CO molecular emission from the tidal arm is entirely redshifted (with the peak velocity at 3.2–3.8 km s−1 with respect to star A; Cabrit et al. 2006). The observed radial velocity is larger than the escape velocity at the corresponding distance, indicating that the motion of the tidal arm is unbound. The tidal arm is also optically thick as the observed 12CO(J = 2–1)/ 12CO(J = 1–0) flux density ratio is ∼4. The tidal arm extends about 4 arcsec (∼600) au at the distance of 140 pc.
2.3 The molecular complex around star B
12CO molecular emission around star B has a broad and asymmetric peak blueshifted by 3–6 km s−1 with respect to star A. The asymmetric peak cannot be attributed to an undisturbed rotating disc. Given that the tidal interaction hypothesis is favoured, the molecular complex around B most likely consists of material captured during the tidal encounter. The 12CO(J = 2–1)/ 12CO(J = 1–0) ratio was observed to be ∼5.7 suggesting that the molecular complex around B is marginally optically thick (Cabrit et al. 2006). Note, however, that the signal to noise of the 12CO(J = 1–0) emission line is very low, we therefore do not use this line as a constraint to discuss our model.
2.4 Proper motion of the system
Bisikalo et al. (2012) compiled the position angles of star B and the angular separations between star A and B over the past 70 years. Linear regression analysis showed that the projected AB separation was increasing at 0.002 arcsec yr−1. 70 years is short compared to the Keplerian orbital period (∼2000 yr) at the observed separation (∼1.5 arcsec, or 200 au for 140 pc away), thus justifying a linear approximation for the velocity. They similarly claimed that the position angle of star B was increasing at a mean rate of 0.02° yr−1, which would imply that the orbital motion of the pair was anticlockwise as projected on to the plane of the sky. This conclusion is however not firm given that the claimed increase is based on a fit involving data from the 1940s whose radial component was discrepant with the recently observed trend of increasing separation. Moreover, this fit to the tangential data is in the opposite sense to the apparent trend observed during the 1990s (although the uncertainties during this period were large). We therefore only constrain our model such that the sky-projected separation of star A and B is increasing at the current epoch. Future proper motion observations are needed to decide whether the motion of the two stars projected into the plane of the sky is clockwise or anticlockwise.
2.5 The dimming event
Rodriguez et al. (2013) reported a long (∼180 d) and deep (2 mag) dimming of RW Aur from 2010 September until 2011 March. Typical photometric variations of Classical T-Tauri Stars (CTTS) only last for days and weeks with a depth <1 mag. This is usually associated with circumstellar extinction or hotspots on the accretion disc that rotate into and out of the observer's view. Rodriguez et al. (2013) noted that there was no similar dimming event3 for the system in the past 50 years and proposed that the dimming may be due to the occultation of star A by the leading edge of the tidal arm. They further estimated that the occulting body was moving with a transverse velocity of 0.86–2.58 km s−1 (dividing the stellar radius of star A by the ingress time of the dimming). The size of the occulting body was estimated to be ∼0.09–0.27 au using the product between the duration of the dimming and the derived velocity. The large range in the estimates of transverse velocity and size are due to the large uncertainties on the ingress time-scale, approximately between 10 and 30 d.
More recently, Antipin et al. (2015) have again detected RW Aurigae in the ‘low state’ (i.e. 3–4 mag less than what observed by HST in the mid-90s; White & Ghez 2001).
2.6 Summary of observational constraints
Among the various observations, we identified the following key criteria to gauge how closely our simulations describe the observations. These features were more reliably constrained by observations; they are also easily quantifiable in simulations:
The size of the system: Cabrit et al. (2006) reported that the tidal arm extended about 4 arcsec (∼600 au at a distance of 140 pc); the projected AB separation was ∼1.5 arcsec (∼200 au). We aim to reproduce the size of the system as observed. However, we argue that it is more important to fit the relative ratio between the dimensions of various structures of the system. This is because the relative ratio is independent of the error in the distance measurement (a parallax of 7.08 ± 0.71 mas corresponding to a distance of 140 ± 14 au; Basri & Bertout 1989).
The relative orientation of the A disc and star B. We gauged the relative orientation using position angles (PA). For the disc, its position angle is defined as the angle between its radial velocity gradient and north. For star B, PA is defined as the angle between the vector joining star A and B and north. According to Cabrit et al. (2006), the PA of the disc is ∼40°; the PA for B is ∼105° (or 255° if using the convention of Bisikalo et al. 2012).
Sky inclination of the disc: the sky inclination of the disc was estimated to be in the range of 45–60° (Woitas et al. 2002; López-Martín et al. 2003).
Proper motion: as reported by Bisikalo et al. (2012), the angular separation between star A and star B was increasing at a rate of 0.002 arcsec yr−1.
The line profiles: the radial velocity of the peak, the peak flux density and the optical depth.
3 NUMERICAL METHOD
3.1 Hydrodynamics models
We model the stellar fly-by using the smooth particle hydrodynamics (SPH) code phantom (see e.g. Lodato & Price 2010; Price & Federrath 2010; Price 2012; Facchini, Lodato & Price 2013; Nixon, King & Price 2013). We assume that self-gravity in the disc is negligible and that the gravitational potential is dominated by the stars. We employ a locally isothermal equation of state and adopt standard artificial viscosity parameters (α = 0.1; β = 2; σ = 0.1) according to the presciption of Morris & Monaghan (1997).
3.2 Synthetic observations
We use the radiation transport and hydrodynamics code torus to compute synthetic observations in this paper. torus is predominantly a Monte Carlo radiation transport code used for calculating line and continuum transport (e.g. Harries 2000; Kurosawa, Harries & Symington 2006) but also features photoionization (Haworth, Harries & Acreman 2012), non-LTE molecular line transfer (Rundle et al. 2010) and hydrodynamics (Haworth & Harries 2012).
We map the result of the SPH simulation on to the torus grid using the technique described in Acreman, Harries & Rundle (2010). We do not explicitly derive the temperature using radiative transfer, but rather assume that the temperature distribution in the model is a function of cylindrical distance from the primary (see below; equation 1), as in Cabrit et al. (2006). For molecular line transfer calculations, we solve the equations of non-LTE statistical equilibrium and produce synthetic data cubes in the manner described in detail by Rundle et al. (2010).
Since we are comparing with the observations from Cabrit et al. (2006), we consider the 12CO(J = 2–1) and (J = 1–0) lines and assume the same parameters, i.e. the 12CO abundance relative to hydrogen is 8 × 10−5 and the system is at a distance of 140 pc from the observer. We assume a small micro turbulence value of 0.03 km s−1. The data cubes comprise 401 × 401 pixels and 36 velocity channels spanning −7.19–7.63 km s−1. This velocity range and number of channels is chosen to match the IRAM Plateau de Bure Interferometer, the instrument used for the observations in Cabrit et al. (2006) against which we compare our models. The beam size is chosen to match that of the IRAM Plateau de Bure Interferometer (2.42 arcsec × 1.55 arcsec and 0.89 arcsec × 0.58 arcsec for the J = 1–0 and J = 2–1 transitions, respectively). We convolve the resulting data cubes with a 2D Gaussian with a width corresponding to the semiminor axis of the two beams using aconvolve from ciao v 4.5 (Fruscione et al. 2006). Finally, the line profiles of the synthetic observations are obtained by integrating over ellipses that match the ones used by Cabrit et al. (2006).
We also compute a 1.3 mm continuum image using Monte Carlo radiative transfer, algorithmic details of which are given in e.g. Harries (2000) and Harries et al. (2004). We discuss the assumed dust distribution in Section 5.2.
3.3 Initial conditions
We model the encounter geometry shown in Fig. 1. in which a spherical polar system (with the z-axis along the normal to the plane of the secondary's orbit, and the x-axis pointing towards the apocentre) is used to define the coordinates (θ, ϕ) of the disc normal (where the separation vector of the two stars at pericentre is located at ϕ = π). We populate particles in a Keplerian disc around the primary according to a prescribed surface density law (see below) and with a Gaussian distribution normal to the disc plane which is the hydrostatic equilibrium solution set by the imposed local temperature (see below). For reasons of computational economy, we set an inner radius of 5.9 au and do not continue to integrate particles that stray within this radius.
We argue that the presence of a single well-defined tidal feature in the observations suggests that the disc of RW Aur A has been subject to either a one-off encounter or to an encounter that is well separated in time from any previous interactions. Thus, we consider it likely that the orbit is either of high eccentricity or unbound. It is unlikely that the encounter is highly hyperbolic since this would imply a relative velocity at infinity that exceeds typical stellar velocity dispersions in star-forming regions. We thus expect that the orbit – whether mildly bound or mildly unbound – is likely to be close to parabolic and hence model this situation here. This is consistent with the low relative radial velocity between the two stars reported by White & Hillenbrand (2014). The masses of RW Aur A and B are set to 1.4 and 0.9 M⊙, respectively (Ghez, White & Simon 1997; Woitas, Leinert & Köhler 2001).
The radial index, q, has been estimated to be in the range 0.5–0.73 by various authors (q = 0.57 and q = 0.73 Osterloh & Beckwith 1995; Cabrit et al. 2006, respectively). We adopt a temperature around 60 K at R = 57 au in line with the estimate of Cabrit et al. (2006) for an inclination of the disc i = 60°, and vary the value of q (which affects the pressure scalelength of the outer disc and has a modest effect on the morphology of the tidal tails produced).
The simulation is initiated when the stellar separation is ≥10 pericentre radii in order to ensure that the disc is not violently perturbed by the sudden introduction of the companion in its vicinity.
3.4 Constraining other parameters
After setting up our model with the plausible initial conditions discussed in the previous section, in principle one could constrain the various other orbital parameters by running the hydrodynamic simulation and synthetic observations in conjunction with a parameter study tool (e.g. Markov chain Monte Carlo). However, such a parameter study is unfeasible in our case: conventionally, 106 particles are needed for sufficient resolution in such a hydrodynamic simulation. In conjunction with the simulated observations, this would take several days of computation per model. Moreover, there are at least 10 remaining parameters in the model and so the time required for such an approach would be prohibitive.
Instead, we adopted the following strategy. We initially ran our hydrodynamic code with 104 test particles on a coarse grid in the 10D parameter space. We then focused on the region that is most consistent with various observational constraints. Parameters were then tuned until simulated values converged to observed values within respective uncertainties. Finally, we repeated the same model with 106 test particles and appropriate scalings (e.g. viscosity scales with number of particles) to reveal small-scale details of the system.
We remind the reader that although we refer to this final model as the best-fitting one, in reality this is only the best-fitting model among the few dozen models we tested. Given the uncertainties in the observations (Section 2) and the various shortcomings of our simulations, the statistical errors introduced from our sampling strategy are likely to be outweighed by other factors (see discussion in Section 5). We nevertheless argue that the good level of agreement between simulations and observations in the models performed is already sufficient to demonstrate that the tidal encounter hypothesis is indeed consistent with the existing observations.
4 RESULTS OF HYDRODYNAMICS MODELS
In this section, we first consider the parameters of the encounter that affect the gross properties of the flow (i.e. the morphology and kinematics of the tidal tail and disc). It is sufficient for such an analysis to be conducted with the hydrodynamic code alone in the limit of optically thin emission. We fix the stellar masses and eccentricity of the orbit for all simulations and for each simulation we adopt given values of θ and ϕ, of pericentre distance, disc outer radius, surface density and temperature normalization, and temperature power-law index. In Section 5, we consider whether a plausible initial disc model can then reproduce line fluxes and line profiles by using synthetic observations for a range of output times and viewing angles (constrained so as to give the correct projected separation of the two stars on the sky). The viewing angle is parametrized in terms of the direction of the line-of-sight vector identified by the polar coordinates θobs and ϕobs in the coordinate frame shown in Fig. 1.
4.1 Sensitivity of results to parameter choice
To determine which model is most consistent with observations, we vary each parameter around the initial estimates obtained from the pilot simulations. We then identify the system observables that are most dependent on each parameter in order to constrain that parameter.
The choice of periastron separation (70 au) was mainly decided by fitting observed length and width of the tidal tail both of which increase with pericentre separation. The length and width of the tidal tail is also dependent on time of observation. However this degeneracy can be broken, since the choice of periastron separation must also give a consistent radial velocity of material in the tidal tail which is similar to the escape velocity at periastron (see Fig. 8). Moreover, since the disruptiveness of the encounter depends on the ratio of the periastron separation to the outer disc radius, the observed scale of the residual disc constrains the initial disc radius (60 au).
Another important quantity that could be used to constrain the orbital parameters of the fly-by is the disc outer radius, after truncation. Since Clarke & Pringle (1993), other groups have worked on the angular momentum and energy transfer mechanisms occurring with a star–disc encounter (e.g. Kobayashi & Ida 2001; Muñoz et al. 2015). Recently, Breslau et al. (2014) focused on the sizes of protoplanetary discs after such encounters using numerical simulations, where the discs are described by test particles (hydro effects are neglected). They model the encounter via prograde coplanar parabolic orbits. They obtain the following prescription for the outer radius of the truncated disc: |$R_{\rm out} = 0.28 R_{\rm peri} m_{12}^{-0.32}$|, where m12 = M2/M1. However, note that they did not explore the effects of the misalignment between the disc and the plane of the parabola. Recent results by Lubow, Martin & Nixon (2015) analysed how the tidal truncation radius estimates are affected by the misalignment between a disc and an external binary (Artymowicz & Lubow 1994). Their results indicate that the external torque decreases quite rapidly with inclination, leading to a larger circumstellar disc. Similarly, we expect the same kind of behaviour for a fly-by interaction, even though in this case resonances do not come into play. Finally, the analytic prescription by Breslau et al. (2014) strongly depends on the definition of the outer radius. They define it as the point where the surface density profile presents a maximum in its radial gradient. However nothing prevents an observer from detecting regions of the disc outside this radius, where the material will be less bound and in more elliptical orbits (see again Breslau et al. 2014). To conclude, the outer radius is not the best parameter to be used to constrain a priori the orbit of the fly-by. We will check a posteriori that what we obtain is in reasonable agreement with the observations (see Section 5).
The mutual inclination of the disc and the binary orbit (θ) is the main determinant of the morphology of the system. While the azimuthal angle ϕ determines the relative position and separation between star A and star B, the polar angle θ determines the size and orientation of the tidal arm. The tidal arm is most prominent and most massive for a prograde coplanar tidal encounter. As θ increases, the tidal arm decreases in size and mass; it also becomes warped with respect to the original disc plane. If θ is larger than 90°, no extensive tidal arm is produced.
We therefore require a prograde encounter. We can obtain a clue about the likely orientation of the disc of A on the sky from the observation of its jet which is directed to the SE. Since the jet must be visible on the observer's side of the disc, it suggests that the NW side of the disc is nearer to the observer. The fact that the NE lobe of A's disc is blueshifted then suggests that the rotation of the disc is clockwise as seen from the earth. In the successful (rather close to coplanar prograde) model, we have found a viewing angle that well reproduces the observed blueshift of the material centred on B. However, the predicted proper motion in the plane of the sky is then clockwise in contrast to that claimed by Bisikalo et al. (2012). We have discussed above how it is desirable to better constrain the motion in the plane of the sky through future proper motion observations.
The power index q of the adopted radial temperature profile determines the rate at which temperature falls as a function of radius. A smaller value of q results in a hotter and thus thicker outer region of the undisrupted disc. Subsequently, a thicker outer disc produces a thicker tidal tail after the tidal interaction. Therefore, the value of q is constrained via the observed ratio between the width and length of the tidal arm.
4.2 The best-fitting model
Fig. 2 shows the gas morphology of the best-fitting model, the parameters of which are listed in Table 1 (compare Fig. 2 with Fig. 3, where we portray fig. 1 by Cabrit et al. 2006). The gas is rendered as individual SPH particles which are colour coded according to line-of-sight velocity. Fig. 4 shows the evolution of the fly-by with rendered plots. Fig. 5 explores the parameter space near the best-fitting model. It shows how observable features depend on variations of each orbital parameters.
Parameter varied . | Best-fitting value . | Effect . |
---|---|---|
Initial disc radius | 60 au | Post-encounter disc size and tidal tail morphology. |
Periastron separation | 70 au | Morphology and kinematics of tidal tail. |
Mutual inclination of disc and binary planes | θ = 18°, ϕ = 170° | θ determines the orientation of tidal tail and |
the captured material on star B; | ||
ϕ determines the velocity and positional angle of star B. | ||
Viewing angle | θobs = 76°, ϕobs = 114° | Morphology in plane of sky. |
Time of observation | 1.1 Pperi after pericentre | Kinematic signatures and ratio of AB separation to |
extent of tidal arm. | ||
Radial temperature profile | |$T = 230\,{\rm K} \left(\frac{R}{{\rm au}}\right)^{-0.5}$| | Power-law index determines aspect ratio of tidal tail. |
Initial disc mass | 0.00 167 M⊙ | Optical depth of different regions of the system. |
Parameter varied . | Best-fitting value . | Effect . |
---|---|---|
Initial disc radius | 60 au | Post-encounter disc size and tidal tail morphology. |
Periastron separation | 70 au | Morphology and kinematics of tidal tail. |
Mutual inclination of disc and binary planes | θ = 18°, ϕ = 170° | θ determines the orientation of tidal tail and |
the captured material on star B; | ||
ϕ determines the velocity and positional angle of star B. | ||
Viewing angle | θobs = 76°, ϕobs = 114° | Morphology in plane of sky. |
Time of observation | 1.1 Pperi after pericentre | Kinematic signatures and ratio of AB separation to |
extent of tidal arm. | ||
Radial temperature profile | |$T = 230\,{\rm K} \left(\frac{R}{{\rm au}}\right)^{-0.5}$| | Power-law index determines aspect ratio of tidal tail. |
Initial disc mass | 0.00 167 M⊙ | Optical depth of different regions of the system. |
Parameter varied . | Best-fitting value . | Effect . |
---|---|---|
Initial disc radius | 60 au | Post-encounter disc size and tidal tail morphology. |
Periastron separation | 70 au | Morphology and kinematics of tidal tail. |
Mutual inclination of disc and binary planes | θ = 18°, ϕ = 170° | θ determines the orientation of tidal tail and |
the captured material on star B; | ||
ϕ determines the velocity and positional angle of star B. | ||
Viewing angle | θobs = 76°, ϕobs = 114° | Morphology in plane of sky. |
Time of observation | 1.1 Pperi after pericentre | Kinematic signatures and ratio of AB separation to |
extent of tidal arm. | ||
Radial temperature profile | |$T = 230\,{\rm K} \left(\frac{R}{{\rm au}}\right)^{-0.5}$| | Power-law index determines aspect ratio of tidal tail. |
Initial disc mass | 0.00 167 M⊙ | Optical depth of different regions of the system. |
Parameter varied . | Best-fitting value . | Effect . |
---|---|---|
Initial disc radius | 60 au | Post-encounter disc size and tidal tail morphology. |
Periastron separation | 70 au | Morphology and kinematics of tidal tail. |
Mutual inclination of disc and binary planes | θ = 18°, ϕ = 170° | θ determines the orientation of tidal tail and |
the captured material on star B; | ||
ϕ determines the velocity and positional angle of star B. | ||
Viewing angle | θobs = 76°, ϕobs = 114° | Morphology in plane of sky. |
Time of observation | 1.1 Pperi after pericentre | Kinematic signatures and ratio of AB separation to |
extent of tidal arm. | ||
Radial temperature profile | |$T = 230\,{\rm K} \left(\frac{R}{{\rm au}}\right)^{-0.5}$| | Power-law index determines aspect ratio of tidal tail. |
Initial disc mass | 0.00 167 M⊙ | Optical depth of different regions of the system. |
This model successfully reproduced the following observed features of the system.
The north-east lobe of the disc is blueshifted with radial velocity of ∼−3 km s−1 with respect to Star A while the south-west lobe of the disc is redshifted with radial velocity of ∼3 km s−1.
The tidal arm has the observed spiral shape with the correct handedness. It is entirely redshifted with a radial velocity of ∼3 km s−1 which is larger than the escape velocity at the corresponding distance. The tidal arm is thus unbound, as observed.
The molecular complex on B is blueshifted with respect to star A. The shape of the captured material around B is similar to that observed with its major axis being tilted slightly to the east of north. This good agreement suggests that the gas observed around B may indeed be material captured from A's disc (even if B also originally possessed a small-scale disc of its own).
Comparing the observation and simulation more quantitatively, Table 1 lists the various quantifiable morphological features. Our simulations matches the observations in nearly all respects. The strikingly good agreement demonstrates that the system's morphology is most likely associated with a tidal interaction.
5 RESULTS OF SYNTHETIC OBSERVATIONS
5.1 Line profiles
We took snapshots from our hydrodynamic models at ∼570 yr (tobs ∼ 1.1Pperi, where Pperi is the Keplerian period at the pericentre) after pericentre and post processed them with torus to produce synthetic 12CO data cubes (c.f. Section 3.2). Fig. 6 shows the integrated 12CO (J = 2−1) emission map of RW Aur after correcting for beam dilution.
The synthetic observations allow us to constrain the temperature normalization and the initial mass of the disc Md, since these two quantities regulate the flux density normalization and the optical depth of the emission lines. In the case that the line is optically thick, the ratio 12CO(J = 2–1)/ 12CO(J = 1–0) should be ∼4 in Rayleigh–Jeans limit. We have therefore multiplied the 12CO(J = 1–0) by 4 for the regions of the tidal arm and the primary disc, to demonstrate whether these components of the system are optically thick. We do not consider the 12CO(J = 1–0) line for the gas around star B, since this line has been marginally detected by Cabrit et al. (2006), as discussed above. The radial velocities were corrected to heliocentric frame by adding the mean heliocentric velocity of star A (15.87 ± 0.55 km s−1; Cabrit et al. 2006).
First, we consider the impact of temperature normalization. By using the temperature profile as estimated by Cabrit et al. (2006, T(R) = 450 K × (R/au)−0.5), the most outstanding discrepancy is that the simulated flux densities are three to five times higher than the observed flux densities. We argue that the most likely cause is that the limitation of BS93 model led to Cabrit et al. (2006) overestimating the temperature of the disc. Although BS93 assumed local thermal equilibrium (LTE), this assumption may not be appropriate for RW Aur. torus is capable of treating both LTE and non-LTE, allowing us to test this. The comparison in Fig. 7 demonstrates that the resultant flux densities are indeed about two to three times lower in the LTE case.
This indicates that in order to explain the observed flux, Cabrit et al. (2006) would have overestimated the temperature by assuming an LTE model. Furthermore, if we use a lower temperature, more consistent with theoretical expectation from stellar irradiation and accretion, (T = 30 K at 40 au; Chiang & Goldreich 1997), the resultant theoretical emission is more compatible with the observed integrated flux densities. Note that lower temperatures are also indicated by Infrared SED fitting by Osterloh & Beckwith (1995, T = 25 K at 40 au).
We try to constrain the disc mass by requiring the tidal arm and the disc orbiting star A to be optically thick. This requires an initial disc mass of 1.67 × 10−3 M⊙ and a current disc mass of 1.05 × 10−3 M⊙. In order to estimate the latter, we have considered an outer disc radius of 57 au, as estimated by Cabrit et al. (2006). This estimate is higher than that estimated by both Cabrit et al. (2006) and Osterloh & Beckwith (1995) (10−4 M⊙ and 3 × 10−4 M⊙, respectively), but lower than the ones estimated by Andrews & Williams (2005, 4 × 10−3 M⊙). Note however that mass estimates from dust and CO isotopologues are always affected by very large uncertainties, and that our prescription for the temperature dependence is so simple that we do not expect to recover the observed mass with high precision.
After adopting the temperature normalization and initial disc mass as described above, the resultant line profiles show better agreement with the observations. Fig. 8 shows line profiles from subsets of these data cubes after applying convolution accounting for the IRAM beam size, for the best-fitting model. In particular, both the A disc and the tidal arm are optically thick as observed, and their peak flux is in broad agreement with the observations by Cabrit et al. (2006), together with the moderate asymmetry of the double-horned line profile of the disc. The disc has in fact gained a small eccentricity, as apparent from the right-bottom panel of Fig. 4. The 12CO(J = 2–1) line profile of the gas around star B recovers the main features of the observed profile, besides the peak flux, which is lower than the observed one by a factor of 2. This may result from incorrect assumptions about the temperature profile around star B (temperature has a simple dependence on the distance from star A), but may also simply reflect that star B also possessed a disc prior to the encounter (in contrast to what we assume here, where the final disc is composed of material entirely captured from the disc of star A).
Cabrit et al. (2006) obtained an indirect measure of the size of the disc around star A by using the position of the peaks of the typical double-horned profile of the emission lines (see Section 2.1). They observed a line profile with peaks at ∼3.5 km s−1 from the centre, whereas in our synthetic observation we obtain the two peaks at ∼2.5 km s−1 from the centre of the line. This seems to suggest that the disc in our simulations is larger than the one observed by a factor ∼2. However, if we look at the hydrodynamical simulation of our best model (see Fig. 9), or if we look at the synthetic continuum emission shown in Fig. 10, the disc size is in good agreement with the observations (∼50 au). After convolution, we are obtaining contamination from low-velocity gas in the inner part of the tidal arm and the elliptical outer region of the disc. The difference in our synthetic line profiles implies that we are overestimating emission from these outer regions. Detailed thermal modelling is however beyond the scope of this project; the agreement with the size derived from the thermal continuum reassures us that the dynamical effects of the encounter are well reproduced.
5.2 Dust emission
We used torus (c.f. Section 3.2) to produce a synthetic observation of the continuum emission at 1.3 mm of the same snapshot showed in Fig. 6. We consider standard assumptions on the dust properties. First of all we assume that spatial distribution of the dust follows exactly the gas distribution. This is likely to be the case in our framework, since the dynamics of the gas itself is mainly ballistic during the encounter. The dust component is dynamically driven by the same gravitational potential. At the moment we neglect secondary effects, that might be due to hydro terms in the dynamics of the gas. We assume a gas-to-dust ratio equal to 100, and a grain size distribution ranging between 50 Å and 1 cm. We use a typical grain size distribution, where dn/ds ∝ s−w, where n is the numerical density of dust particles, and s is the grain size. It is well known that in the interstellar medium w assumes a typical value w ∼ 3.5 (Mathis, Rumpl & Nordsieck 1977). More recently, interferometric observations have suggested a shallower distribution in the outer regions of protoplanetary discs, where w ∼ 3 (e.g. Ricci et al. 2010a,b, 2012; Testi et al. 2014). Following these recent results, for our synthetic observation we use w = 3.3. We assume that the dust is comprised of amorphous silicate grains (Draine 2003a,b), obtaining an opacity κν = 0.131 cm2 g−1 at 1.3 mm (assuming a gas-to-dust ratio equal to 100). Fig. 10 shows the emission we obtain in mJy beam−1, where the beam size is the one used by Cabrit et al. (2006) (0.87 arcsec × 0.56 arcsec, see their table 1). Contour levels are submultiples and multiples of the noise level σ = 0.47 mJy beam−1 in the same observation by Cabrit et al. (2006): (1/24, 1/6, 1, 3, 48)σ. From disc A, we obtain emission levels that are compatible with the one observed (given the large uncertainties on the opacity), and a disc size that is in agreement with the one estimated by Cabrit et al. (2006) from the line profiles (∼50 au from the unconvolved emission map). We underestimate the dust emission around star B, and this is again suggesting that our simple temperature treatment is underestimating the temperature of the material around the secondary star. By comparing the convolved emission map with the one shown by Cabrit et al. (2006) where the primary disc is located, it seems that we are overestimating the total emission by a very small factor ∼1.5 (compare the size of the 3σ contours), which is very small when compared to the uncertainties of the problem. Note that the normalization of the emission of our synthetic observation does depend on our assumptions, mostly the dust opacity (e.g. in discs is often assumed a lower opacity at 1.3 mm, as motivated by Beckwith et al. 1990), but also the uniform dust-to-gas ratio, the maximum grain size and the grain size distribution index. With the assumptions we have made, the emission from the tidal tail is below the noise level (∼0.3 mJy beam−1).
6 THE DIMMING EVENT
Rodriguez et al. (2013) argued that the leading edge of the tidal arm may wrap around the A disc and obscure the line of sight between Earth and RW Aur A, thereby producing the observed dimming of RW Aur. In our simulations, there is indeed a stream of gas linking star A and B after the tidal encounter (however, it is not part of the primary tidal arm as Rodriguez et al. 2013, suspected, see the bottom-left panel of Fig. 4). A portion of this stream of gas might provide the occulting body that caused the observed dimming event, as the integrated line-of-sight column density (∼10−2 g cm−2) exceeds the threshold column density (∼7 × 10−4 g cm−2) required to produce the dimming – assuming a flux decrease of the A star of 91 per cent (Rodriguez et al. 2013) and a dust opacity at optical wavelengths κ = 3000 cm2 g−1 (Men'shchikov & Henning 1997).
In order to verify this hypothesis, we have analysed the SPH simulation of the best-fitting model. We select all the foreground particles that are partially occulting star A at the present time by requiring h/D > 1, where h is the smoothing length associated with a single SPH single particle and D is the distance between the particle and star A in the plane of the sky. We then trace back their trajectories, in order to see whether the selected particles have a credible dynamical history. In Fig. 11, we plot all the SPH particles of the simulation with blue crosses, the two stars with yellow star symbols and the occulting particles with filled dots. The orange dots represent the particles that show a credible dynamical history. They start orbiting together within the circumstellar disc (see left-hand panels), and they then get incorporated in the well-resolved bridge linking star A and B. At the present time, they are located at the low-density periphery of the bridge (see right-hand panels). The occulting particles represented by the red dots are individually scattered particles that originate from the very poorly resolved outer regions of the initial disc and are not associated with any coherent fluid structure. Since the orange particles are very few (only four for the particular choice of viewing angle), we cannot determine the structure of the low-density margins of the bridge, and even less the size of the occulting body (0.09–0.27 au; as estimated by Rodriguez et al. 2013). An accurate determination of the density structure of the low-density periphery of this bridge would require an unfeasibly high resolution. However, the presence of particles in this region indeed indicates that possible ∼ au scale inhomogeneities at the edge of the stream linking A and B could cause the dimming event observed by Rodriguez et al. (2013). The velocity in the plane of the sky of these particles ranges between 0.97and1.09 km s−1 and is thus consistent with the estimates from Rodriguez et al. (2013, 0.86–2.58 km s−1).
Finally, we have computed the evolution with time of the average column density lying between star A and the observer. In Fig. 12, we show that in the past the star has been completely occulted (at optical wavelengths) by the density structure linking stars A and B. This bridge has then evolved, until the line of sight only grazes it nowadays (t = 0). In the future, the average column density stays well below the value required to observe a dimming of star A as observed by Rodriguez et al. (2013). Moreover, for t ≳ 0, the column density in the simulation starts to be dominated by particle noise. However, inhomogeneities in the bridge structure might entail new events in the next years.
7 SUMMARY AND CONCLUSIONS
We have used hydrodynamic models and synthetic observations to constrain the observed morphology, kinematics and line fluxes of the binary system RW Aur which was first proposed as a candidate star–disc encounter by Cabrit et al. (2006). We draw the following conclusions from this work.
With a suitable set of orbital parameters and viewing angle (see Table 2), we can reproduce almost all morphological and dynamical features of the system.
RW Aur A has a circumstellar disc with an outer radius of ∼50 au which is within the observed range of 40–57 au. The disc is optically thick in 12CO emission, as the 12CO(J = 2–1)/ 12CO(J = 1–0) ratio is ∼4.
RW Aur B is separated from star A by ∼200 au (∼1.5 arcsec for 140 pc). The projected AB separation is increasing at rate of ∼0.002 arcsecyr−1. The molecular complex on B is blueshifted with respect to star A with a broad distribution of radial velocity of 3–6 km s−1. The line profile from the star B is asymmetric and thus most likely due to material captured during tidal encounter instead of a circumstellar disc. The molecular emission is marginally optically thick with a 12CO(J = 2–1)/ 12CO(J = 1–0) > 4.
The expanding, spiral-shaped structure connected to RW Aur A disc is most likely a tidal arm produced in a fly-by of star B. The tidal arm extends about 600 au (∼4 arcsec for 140 pc). The arm is entirely redshifted with a radial velocity ∼3 km s−1 with respect to star A. This is faster than the escape velocity at corresponding distance, indicating the arm is indeed gravitational unbound.
Observed features . | Observed values . | Best-fitting . |
---|---|---|
. | . | simulation values . |
Length of tidal arm | 600 au (4 arcsec)a | 650 au |
AB separation on sky | 200 au (1.5 arcsec)a | 200 au |
PA of disc (°) | 40a | 47.2 |
PA of AB (°) | 105a | 105.9 |
Inclination of disc (°) | 45–60b | 64.0 |
Rate of change of | 0.002 arcsecyr−1c | 0.003 arcsec yr−1 |
Observed features . | Observed values . | Best-fitting . |
---|---|---|
. | . | simulation values . |
Length of tidal arm | 600 au (4 arcsec)a | 650 au |
AB separation on sky | 200 au (1.5 arcsec)a | 200 au |
PA of disc (°) | 40a | 47.2 |
PA of AB (°) | 105a | 105.9 |
Inclination of disc (°) | 45–60b | 64.0 |
Rate of change of | 0.002 arcsecyr−1c | 0.003 arcsec yr−1 |
Observed features . | Observed values . | Best-fitting . |
---|---|---|
. | . | simulation values . |
Length of tidal arm | 600 au (4 arcsec)a | 650 au |
AB separation on sky | 200 au (1.5 arcsec)a | 200 au |
PA of disc (°) | 40a | 47.2 |
PA of AB (°) | 105a | 105.9 |
Inclination of disc (°) | 45–60b | 64.0 |
Rate of change of | 0.002 arcsecyr−1c | 0.003 arcsec yr−1 |
Observed features . | Observed values . | Best-fitting . |
---|---|---|
. | . | simulation values . |
Length of tidal arm | 600 au (4 arcsec)a | 650 au |
AB separation on sky | 200 au (1.5 arcsec)a | 200 au |
PA of disc (°) | 40a | 47.2 |
PA of AB (°) | 105a | 105.9 |
Inclination of disc (°) | 45–60b | 64.0 |
Rate of change of | 0.002 arcsecyr−1c | 0.003 arcsec yr−1 |
The synthetic molecular line data cubes obtained from a purely radial temperature estimated by Cabrit et al. (2006) have flux densities three to five times higher than observations. Nonetheless, by lowering the temperature to typical disc temperature for CTTS, we can reduce the flux densities to the observed values.
The agreement between observations and simulations thus lends strong support to the hypothesis that the morphology of the RW Aur system is consistent with a tidal encounter scenario. Although the observational evidence for tidal encounters on Galactic scales is well established, this is the only system that is a good candidate for such encounters on the scale of protoplanetary discs (note that a tidal encounter scenario has also been proposed as a possible model to explain the recent ALMA observations of AS 205 S; see Salyk et al. 2014).
Finally, we note that our simulations provide support for the hypothesis of Rodriguez et al. (2013) regarding the origin of observed dimming of star A in 2010/2011. We find that it is plausible that this was caused by inhomogeneities in the outskirts of a ‘bridge’ structure linking star A and B which may graze the line of sight at current epochs. Our simulations do not have the resolution to trace the structure of the low-density periphery of the optically thick bridge, nor to model the ∼ au scale inhomogeneities that would be required to explain such an event. Nevertheless our simulations predict good agreement between the velocity in the plane of the sky of material in this region with that required by the observed ingress time of the dimming event. This provides good support for the occultation hypothesis. It suggests that – depending on the filling factor of such inhomogeneities – further dimming events should not be ruled out (as confirmed by the observations by Antipin et al. 2015) and that these would be expected to share similar kinematic properties.
We thank the anonymous referee, whose comments helped improving this manuscript. We thank Daniel Price for giving us the permission to use the phantom code, Joey Rodriguez and Keivan Stassun for a careful reading of the manuscript, and Sylvie Cabrit and Tim Harries for useful discussions. The synthetic observation calculations in this paper were performed using the Darwin Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E–infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E–Infrastructure. SF is funded by an STFC/Isaac Newton Trust studentship. Fig. 4 was produced using splash (Price 2007), a visualization tool for SPH data. This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG.
Note that the long-term light curve does show a number individual photometric data points at the ‘low’ level recorded in the 2011 event but a well-resolved eclipse light curve has not been detected hitherto.