The observable effects of tidally induced warps in protostellar discs Free

Abstract

1 INTRODUCTION

Protostars have circumstellar discs from which they accrete matter. Once the local accretion rate on to the protostar falls below a certain value, given approximately by

1

where M is the stellar mass, R_* the stellar radius and L_* the stellar luminosity, then most of the energy radiated from the disc comes from re-radiated stellar emission, rather than its own internally generated accretion energy. Such discs are known as passive discs and most observed protostellar discs fall into this category. The nature of the spectrum emitted by the disc depends on how much of the stellar flux can be intercepted and at what radii this occurs. A flat, infinitesimally thin disc, has at large radii a temperature distribution of the form T_d∝R^−3/4 which gives rise to a spectral energy distribution (SED) of the form λF_λ∝λ^−4/3. Almost all protostellar SEDs are flatter than this, implying that real discs intercept substantially more of the stellar flux. For this reason, Kenyon & Hartmann (1987; see also the review by Dullemond et al. 2007) argued that protostellar discs are substantially flared, that is, that the disc scaleheights H are such that d(H/R)/dR > 0. Since then more sophisticated models of the flaring have been introduced (Chiang & Goldreich 1997; D'Alessio et al. 1999; Chiang et al. 2001; Dullemond & Dominik 2004a,b), together with truncated discs which display a heated inner rim (Dullemond, Dominik & Natta 2001; Natta et al. 2001; Muzerolle et al. 2003; Isella & Natta 2005). These models are generally successful in providing explanations of the observed SEDs in the near- and mid-infrared.

These models, however, all assume the disc to be cylindrically symmetric, and this seems to be a reasonable assumption in the inner regions (Monnier et al. 2006). Stars typically form within clusters, and so even single stars are likely to have undergone close encounters with neighbouring stars at a young age. Very close encounters are likely to be dissipative and lead to disc truncation and perhaps even capture and binary star formation (Clarke & Pringle 1991, 1993; Hall, Clarke & Pringle 1996; Boffin et al. 1998; Watkins et al. 1998a,b; Pfalzner et al. 2005; Moeckel & Bally 2006). More distant encounters occur more frequently, but can still provide the disc with significant warping. We note that there is observational evidence for such disc disturbances (Cabrit et al. 2006; Lin et al. 2006; Quillen 2006) and for the warping of the disc at large radius (Hughes et al. 2009).

Viscosity in the outer regions of protostellar discs is likely to be small, and typical figures of the dimensionless viscosity parameter α (Shakura & Sunyaev 1973) in the range 10⁻²–10⁻⁴ are thought likely, e.g. Terquem (2008). This is less than the disc thickness ratio H/R, and in this case such warps are propagated as waves (Papaloizou & Lin 1995). With such low viscosities, these waves can exist for a significant fraction of the protostellar lifetime (Section 2).

In this paper, we investigate the effects of these waves in a very simplified manner. We do not attempt to model perturbations to the structure or SEDs of discs with flares or other configurations. Rather we consider the effect on the shape and on the resultant SEDs emitted by a disc which is initially thin and flat. Although this is not likely to enable us to provide realistic fits to any observations, it does enable us to see at what wavelengths any complications deriving from such perturbations might manifest themselves, and to provide an estimate of their relative magnitude. We note that the effect of warping on the SED of a protostellar disc has been considered previously for the case of a disc in a misaligned binary system by Terquem & Bertout (1993, 1996), but that they do not consider the dynamic response of the disc, nor its subsequent evolution towards alignment (e.g. Bate et al. 2000).

In Section 2, we provide estimates of the periods and lifetimes of warp-like perturbations in a protostellar disc. In Section 3, we compute the time-dependent behaviour of warps which can be generated using a simple model of the flyby of a perturbing star. In Section 4, we compute the temperature distributions and resulting SEDs for these disc models. We discuss the implications of our results in Section 5.

2 WARP PROPAGATION AND LIFETIME ESTIMATES

We consider a protostellar disc in orbit around a star of mass M. For a Keplerian disc, the angular velocity of disc material is Ω= (GM/R³)^1/2. We take the disc semithickness to be H, so that at radius R

2

where c_s(R) is the local disc sound speed (Pringle 1981). Typical values of H/R are expected to be around 0.05–0.1 (Bell et al. 1997; Terquem 2008), independent of radius R. Typical estimates for the dimensionless viscosity parameter α for such discs are in the range α≈ 10⁻²–10⁻⁴ (Hartmann 2008). Thus, for protostellar discs, we expect that α < H/R, and therefore that warps propagate through them in a wave-like manner (Papaloizou & Lin 1995). The propagation speed for these waves is (Papaloizou & Lin 1995; Pringle 1999; Lubow, Ogilvie & Pringle 2002)

3

Suppose that a warp is induced in the outer disc regions by the fly-by of some other stellar object. If the relative velocity of the encounter is V, and the distance of closest approach is D then the characteristic frequency of the induced warp-like disturbance is

4

and the wavelength of the disturbance is

5

We expect the warp to propagate inwards to a radius R_crit at which R=λ(R). At that radius the propagating warp wave is reflected, and within that radius the disc is flat, but tilted by the warp. This is because the wavelength of the warp is greater than the size of the disc in the inner regions and therefore cannot communicate the warp. However, angular momentum is still communicated into the inner disc, so the disc tilts without warping. If H/R≈ const., independent of radius, then we find that at radius R_crit

6

If we make the further assumption that the fly-by is parabolic, so that to a first approximation

7

then we find that

8

Thus, we expect that while the inner disc regions R < R_crit will have a time-varying tilt, it is only the outer disc regions R > R_crit that will show a significant warp. As a numerical example, if we assume that the fly-by passes the disc at a distance of D= 250 au, and H/R= 0.1 then we have R_crit≈ 34 au.

The time taken for the warp to propagate to the disc centre and back is

9

where R_out is the radius of the outer disc edge. Equivalently, this may be written as

10

where P_out is the orbital time at the outer disc edge. Thus for an outer disc edge at 100 au, and

11

the warp crossing time-scale is around 10⁴ yr.

The decay time-scale for the warp is given approximately by (Lubow et al. 2002)

12

Thus for a disc with R_out= 100 au, we would expect the warp to last for a time of around 1.6 × 10⁶(α/10⁻⁴)⁻¹ yr.

3 WARP GENERATION – NUMERICAL RESULTS

3.1 The warp equations

We model the generation and propagation of the warp using the linearized equations derived by Lubow & Ogilvie (2000). These describe the warp in terms of the local unit tilt vector l(R, t) which is a function of radius and time. For a Keplerian disc these are

13

and

14

Here G is the internal torque which acts to realign tilted disc annuli. T is the external torque acting on the disc which in this case is caused by the stellar fly-by. In addition, Σ is the disc surface density, and is the vertically integrated pressure. We note that these equations do not make allowance for dispersive or non-linear effects in the wave propagation (Ogilvie 2006). For the waves, we discuss here dispersive effects may be of marginal significance, but only at the outer disc edge. However, for the large amplitudes required to substantially change the SED, non-linear effects may well play a role. However, to take these effects properly into account would require a full hydrodynamical model of the disc, which is beyond the scope of the current investigation.

To solve these equations numerically we adopt the approach of Lubow et al. (2002). The disc is taken to lie in the XY-plane, so that for small tilts we have to first order

15

where l_x, l_y≪ 1. The torques and lie in the XY-plane. We can then use W=l_x+i l_y, G=G_x+i G_y and T=T_x+i T_y as complex representations of l, G and T, respectively. We then replace the quantities W and G by the quantities A(R, t) and D(R, t) defined by

16

and

17

where the asterisk denotes a complex conjugate.

Equations (13) and (14) then become

18

and

19

We truncate the disc at R=R_out and do not allow for dissipation at the outer boundary. This is because the wavelength of the warp is long in comparison to any smoothing of the outer edge and so the warp will only see it as a sharp edge and reflect back inwards.

3.2 External torque due to fly-by

We now obtain an expression for the torque density T(R, t) exerted upon an annulus of the disc during a stellar flyby. The torque density exerted on the annulus at radius R with tilt vector by a mass of M₂ at position vector R_b is given to first order in R/R_b by (Lubow & Ogilvie 2000)

20

where b⁽¹⁾_3/2 is the Laplace coefficient. We approximate this in the form

21

We note that since we are using the approximation of only considering linear warp waves, it is only the T_x and T_y components which are relevant for the computations here.

3.3 The disc shape and its evolution

We start with a disc which has an outer radius of R_out= 100 au. We take the inner radius to be at R_in= 0.01 au, and note that this is much less than the expected value of R_crit (equation 11) within which we expect the disc to be tilted but not warped.

We take the surface density of the disc to be a power law of the form

22

(cf. Andrews et al. 2009; Isella, Carpenter & Sargent 2009) and take H/R= 0.1 = const. This then implies that the sound speed is a power law of the form:

23

We apply a small amount of damping by taking α= 10⁻⁴. We note that for these disc properties, we expect the warp amplitude to vary in such a manner as to keep the flux (F_⊥) of angular momentum associated with the warp constant, where

24

For the disc properties, we have chosen this corresponds to |W| ≈ const., so that the vertical local disc displacement due to the warp is roughly proportional to radius.

The disc evolution equations (18) and (19) are evolved numerically using a leapfrog scheme. Boundary conditions are chosen to ensure that the inner and outer edges of the disc are stress free.

To illustrate the possibilities, we consider two idealised flybys. We take the orbit of the perturber to lie along a straight line, constant velocity path. In practice, of course, the interaction is more dynamical than this (e.g. Cabrit et al. 2006; Moeckel & Bally 2006), however we can, to first order, obtain the correct warping structure for the disc. Thus, in both cases, we take the path to be of the form

25

Here, t₀ is some arbitrary time at which the perturber is at position vector ⁠.

3.3.1 Model 1

For Model 1, we take the flyby to lie along a straight line perpendicular to the disc of the form

26

Here, t₀ is the time of closest approach when the perturber is a distance of a= |a| = 300 au from the disc centre, and 200 au from the disc edge. As can be seen, the velocity is in the Z-direction. We take its magnitude V to correspond to the escape velocity from the point of nearest approach, so that

27

Because the path of the perturber lies solely in the XZ-plane, the disc annuli are tilted only about the Y-axis. Thus we can envisage the disc warp by considering a cut through the disc in the XZ-plane.

The results for Model 1a for which we take the mass of the perturber to be M₂= 1 M_⊙ are shown in Fig. 1. At time t= 0 the disc is assumed unperturbed and the perturber has a z-coordinate of z=− 1000 au, and thus t₀≈ 2000 yr. The outer edge of the disc is initially pulled downwards. Then as the perturber passes the disc (passing through the initial disc plane, at time t=t₀≈ 2000 yr) the outer edge of the disc starts to be pulled upwards. Thus, the disc is set oscillating with an amplitude at the outer edge of around 5 au. The period of the oscillation is

28

For the disc considered here this gives P≈ 2100 yr.

Figure 1

A sample of the XZ-profiles of the disc at various times for Model 1a, showing the induced warp and its propagation through the disc. For time t= 0, the disc is entirely flat in the XY-plane. Here we only show the positive X-axis. In the XZ-plane, the disc is antisymmetric around x= 0. The axes are in au. Note that the inner region of the disc (R≲ 30 au) is tilted but not warped.

Open in new tabDownload slide

For Model 1b we explore a more substantial warping of the disc. For this, we use a= |a| = 200 au and we take the mass of the perturber M₂= 10 M_⊙. We note that this changes the velocity of the flyby to V= 9.5 km s⁻¹. As can be seen in Fig. 2 this produces an amplitude for the warp wave of approximately 35 au.

Figure 2

An example of the XZ-profile of the disc for Model 1b, showing that in this case the warp is more substantial at the outer edge of the disc. The axes are in au.

Open in new tabDownload slide

3.3.2 Model 2

In Model 2, we take the path to lie in the direction ⁠, that is, at an angle to the disc plane. The perturber is initially at position

29

The perturber is then given a velocity of V= |V| = 2.7 km s⁻¹ which corresponds to the escape velocity at the point of closest approach R_b= (200, −100, −100) au. Note that this path passes through the XY-plane at the same point as the path in Model 1a. Because this path is not in the XZ-plane, the torque acting on disc annuli is now no longer simply in the Y-direction. Thus, the disc is not only tilted but also acquires a twist. For comparison with Model 1a, we take M₂= 1 M_⊙.

As the disc is twisted we can no longer view it via a cut through the XZ-plane. To illustrate this, in Fig. 3, we show the disc in three-dimensions, for Model 1b, and compare it to Fig. 4, which is a typical shape of the disc for Model 2. The disc shapes for Model 1a and 1b are the same but with differing amplitudes, so we use Model 1b here to illustrate the disc shape more easily.

Figure 3

We show here a representative disc shape in three-dimensions for Model 1b. The image was generated using the graphics package povray. The top image is the XZ profile of the 100 au disc, viewed and illuminated from (0, −200, 0). The second image is the same disc but viewed from (0, −200, 50), and illuminated from (0, 5, 0). Because the path of the perturber lies in the XZ-plane, each disc annulus is tilted about the Y-axis. Thus the disc is warped but not twisted.

Open in new tabDownload slide

Figure 4

We show here a representative disc shape in three-dimensions for Model 2. The image was generated using the graphics package povray. The top image is the XZ profile of the 100 au disc, viewed and illuminated from (0, −200, 0). The second image is the same disc but viewed from (0, −200, 50), and illuminated from (0, 5, 0). Because the path of the perturber is not perpendicular to the unperturbed disc, the individual disc annuli are tilted about different axes in the XY-plane. This means that the shape of the disc is twisted about the Z-axis. The central region is tilted in one plane and the outer rings are slowly twisted away from this plane, producing the shadowing effects.

Open in new tabDownload slide

4 DISC HEATING AND RESULTANT SEDs

4.1 Temperature profile

In order to keep things simple, we assume that the disc is thin, and therefore ignore any flaring. In this way, we ensure that any re-emission of radiation from the disc at large radii (long wavelengths) is caused predominantly by the warp. As we have noted, the inner disc (well within the radius R_crit≈ 30 au) remains flat, but with variable tilt. Thus for the inner disc we may take its surface effective temperature to be T_d(R) where

30

where T_* is the stellar effective temperature and R_* the stellar radius (Friedjung 1985; Kenyon & Hartmann 1987).

For the warped part of the disc, we write the local unit vector representing the disc tilt as

31

so that the angle β(R) represents the tilt of the disc normal, relative to the Z-axis, and the angle γ(R) represents the azimuth of the tilt. Since the warp starts only at a distance R≥R_crit≫R_* we are able to treat the illumination of the warped part of the disc by the central star as coming from a point source. In this case (cf. Pringle 1996) the disc surface temperature can be written as

32

This is valid only if self-shadowing of the disc can be ignored. In practice, this needs to be taken into account. Thus this temperature is applied only to those parts of the disc which are tilted towards the central star and which have a view of the central star which is unobscured by other (warped) parts of the disc. Parts of the disc which are obscured are assumed to have negligible temperature. This computation is carried out in a straightforward manner as in Pringle (1997).

Once the temperature distribution of the disc is found, we compute the emitted disc spectrum simply by assuming that each disc element emits locally as a blackbody. For the SED, we assume an observer placed on the Z-axis for simplicity. If the disc was to be viewed from an angle, there would be an apparent reduction in the magnitude of the flux emitted by the disc. If that angle is increased enough, then the central star and large portions of the disc would be obscured from the observer by the warp.

4.1.1 Model 1

In Fig. 5, we show a representative temperature distribution found for Model 1b, and in Fig. 6 the resultant SED (as viewed by a distant observer in the direction of the Z-axis) from the disc structure shown in Fig. 3. For comparison we also plot, in Fig. 6, an SED from Model 1a and the contribution to the SED from the central star. Because of the induced warp, the SED displays a secondary peak at around 100 μm. For shorter wavelengths, the radiation comes from the stellar surface and the flat illuminated part of the inner disc. The magnitude of the warp is larger in Model 1b, creating a more substantial flux from the warped disc at around 100 μm as opposed to the weaker flux from Model 1a. In Model 1, all disc annuli are tilted about the same axis, therefore outside the central region (taken to be at the reflection radius), only one half of the disc appears illuminated.

Figure 5

Temperature profile for the disc in Fig. 3, from Model 1b. The axes are in au, with the temperature coloured by the scale shown in Kelvin. The central star (with T_*= 5000 K) and innermost parts of the disc are omitted to give a higher resolution over the disc. Note that only one half of the disc appears illuminated to an external observer.

Open in new tabDownload slide

Figure 6

SED in terms of λF_λ. The solid line is for the disc in Fig. 3, from Model 1b. The dot–dashed line is a typical SED from Model 1a, and the dashed line represents the star which is taken to be a blackbody with T_*= 5000 K. The contribution from the star peaks at about 0.6 μm. The inner disc has a similar maximum temperature and a log wavelength dependence of the form λF_λ∝λ^−4/3. The contribution from the upturned part of the outer disc peaks at around 100 μm. As the warp is more substantial in Model 1b, the flux in this peak is greater than the flux in the corresponding peak for Model 1a. As is to be expected, both disc peaks have less magnitude than the primary peak, with Model 1b having a disc peak flux that is around 5 times less than the primary peak flux.

Open in new tabDownload slide

In general, the disc illumination can be more complicated than this because of self-shadowing, and we illustrate some possibilities in Fig. 7.

Figure 7

Some examples of the complicated illumination profiles that arise due to the self-shadowing as the warp propagates through the disc. Note that the characteristic property of an illuminated warp is that the disc surface brightness is lop-sided with respect to the central star. The images were generated by the graphics package povray, with the source of illumination on the central stars surface. The discs here have a radius of 100 au.

Open in new tabDownload slide

4.1.2 Model 2

As the flyby in Model 2 is misaligned with the original plane of the disc, the resultant disc is both tilted and twisted. This can be seen in Fig. 8 which shows a typical temperature profile for Model 2. This exhibits a spiral pattern in the illumination profile. The resultant SED, Fig. 9, is not changed greatly by this, which is to be expected as the disc intercepts a similar amount of stellar flux when both tilted and twisted as opposed to just being tilted. The contribution of the heated outer disc again peaks at around 100 μm.

Figure 8

Temperature profile for the disc in Fig. 4, from Model 2. The axes are in au, with the colour scale giving the temperature in Kelvin. The central star (with T_*= 5000 K) and innermost parts of the disc are omitted to give a higher resolution over the disc. Note that when the path of the perturber is not perpendicular to the unperturbed disc plane, the disc is twisted and so the pattern of illumination displays an apparent spiral pattern.

Open in new tabDownload slide

Figure 9

SED in terms of λF_λ for the disc in Fig. 4, from Model 2.

Open in new tabDownload slide

5 DISCUSSION

We have considered the response of a protostellar disc to a tidally induced warp, and the resultant change in the SED. We have used a simplified analysis in order to emphasize the relevant implications.

For the low viscosities currently envisaged for the outer disc regions, we note that a tidally induced warp is propagated as a wave, and that such a wave can be comparatively long lasting (in our case more than a Myr for α= 10⁻⁴). The warp is largest in the outer regions, and within some critical radius (in our case R_crit∼ 30 au) the disc remains flat but has a tilt which is variable on a time-scale of order the propagation time-scale of the warp (here a few thousand years). Note that any collimated outflow, or jet, originating from the inner disc regions might be expected therefore to vary direction on this time-scale.

In reality, the passive disc that we are considering here will be flared. Flaring of the disc contributes to the SED at wavelengths longer than a few μm. In this paper, in order to emphasize the effects of the warp, we have assumed the disc to be locally thin when computing the SED. In practice the two effects, that of the flare and that of the warp, must be combined.

From our computations, in agreement with Terquem & Bertout (1993, 1996), we conclude that the warp is only likely to affect the disc SED at wavelengths of around 100 μm and longer, but that the effects can be substantial. We also note that warp-induced variations in disc illumination have the effect of making the disc appear asymmetric with respect to the central star at radii R≳R_crit. Further, a twist in the disc, induced by a skewed flyby, can have the effect of introducing a spiral component into the disc illumination.

We thank Cathie Clarke for useful discussions on the topic. We also thank the referee for useful feedback.

REFERENCES