Evolution of Spin Direction of Accreting Magnetic Protostars and Spin-Orbit Misalignment in Exoplanetary Systems: II. Warped Discs

Magnetic interactions between a protostar and its accretion disc tend to induce warping in the disc and produce secular changes in the stellar spin direction, so that the spin axis may not always be perpendicular to the disc. This may help explain the recently observed spin-orbit misalignment in a number of exoplanetary systems. We study the dynamics of warped protoplanetary discs under the combined effects of magnetic warping/precession torques and internal stresses in the disc, including viscous damping of warps and propagation of bending waves. We show that when the outer disc axis is misaligned with the stellar spin axis, the disc evolves towards a warped steady-state on a timescale that depends on the disc viscosity or the bending wave propagation speed, but in all cases is much shorter than the timescale for the spin evolution (of order of a million years). Moreover, for the most likely physical parameters characterizing magnetic protostars, circumstellar discs and their interactions, the steady-state disc has a rather small warp, such that the whole disc lies approximately in a single plane determined by the outer disc boundary conditions, although more extreme parameters may give rise to larger disc warps. In agreement with our recent analysis (Lai et al. 2010) based on flat discs, we find that the back-reaction magnetic torques of the slightly warped disc on the star can either align the stellar spin axis with the disc axis or push it towards misalignment, depending on the parameters of the star-disc system. This implies that newly formed planetary systems may have a range of inclination angles between the stellar spin axis and the symmetry axis of the planetary orbits.


INTRODUCTION
In a recent paper [Lai et al. (2010), hereafter Paper I], we proposed a novel mechanism for producing misalignment between the spin axis of a protostar and the normal vector of its circumstellar disc.Our work was motivated by recent measurements of the sky-projected stellar obliquity using the Rossiter-McLaughlin effect in transiting exoplanetary systems, which showed that a large fraction of the systems containing hot Jupiters have misaligned stellar spin with respect to the planetary angular momentum axis [see Triaud Email: fvf2@cornell.eduet al. (2010); Winn et al. (2010) and references therein].Additional evidence for nonzero stellar obliquity came from the statistical analysis of the apparent rotational velocities (v sin i ) of planet-bearing stars (Schlaufman 2010).
The basic mechanism ("Magnetically driven misalignment") for producing spin -disc misalignment in accreting protostellar systems can be sumarized as follows (Paper I).The magnetic field of a protostar (with B > ∼ 10 3 G) penetrates the inner region of its accretion disc.These field lines link the star and the disc in a quasi-cyclic fashion (e.g., magnetic field inflation followed by reconnection; see Bouvier et al. (2007) and Alencar et al. (2010) for observational evidence).Differential rotation between the star and the disc not only leads to the usual magnetic braking torque on the disc, but also a warping torque which tends to push the normal axis of the inner disc away from the spin axis1 .Hydrodynamical stresses in the disc, on the other hand, tend to inhibit significant disc warping.The result is that, for a given disc orientation imposed at large radii (e.g., by the angular momentum of the accreting gas falling onto the disc), the back-reaction of the warping torque can push the stellar spin axis toward misalignment with respect to the disc normal vector.Planets formed in the disc will then have a misaligned orbital normal axis relative to the stellar spin axis, assuming that no evolution mechanism occurring after the dissipation of the disc forces the alignment of the system.
The process of planetary system formation can be roughly divided into two stages (Juric & Tremaine 2008).In the first stage, which lasts a few million years until the dissipation of the gaseous protoplanetary disc, planets are formed and undergo migration due to tidal interactions with the gaseous disc [Lin et al. (1996); see Papaloizou et al. (2007) for a review].The second stage, which lasts from when the disc has dissipated to the present, involves dynamical gravitational interactions between multiple planets, if they are produced in the first stage in a sufficiently closepacked configuration (Juric & Tremaine 2008;Chatterjee et al. 2008), and/or secular interactions with a distant planet or stellar companion (Eggleton & Kiseleva-Eggleton 2001;Wu & Murray 2003;Fabrycky & Tremaine 2007;Wu et al. 2007;Nagasawa et al. 2008).The eccentricity distribution of exoplanetary systems and the recent observational results on the spin -orbit misalignment suggest that the physical processes in the second stage play an important role in determining the properties of exoplanetary systems.Nevertheless, the importance of the first stage cannot be neglected as it sets the initial condition for the possible evolution in the second stage.Our result in Paper I shows that at the end of the first stage, the symmetry axis of the planetary orbit may be inclined with respect to the stellar spin axis.
At first sight, it may seem strange that the magnetic field effects can drive the stellar spin axis toward misalignment with respect to the disc symmetry axis, given that the spin angular momentum of the star ultimately comes from the disc and the disc contains a large reservoir of angular momentum.The key to understand this is to realize that when the gas reaches the magnetosphere boundary, its angular momentum is much smaller than in the outer disc (the specific angular momentum of the disc is j disc (r) = √ GM r for a Keplerian disc), and any magnetic torque, which in general can break the axisymmetry of the system, is of the same order of magnitude as the accretion torque on the star.
A key assumption adopted in Paper I for the calculation of the magnetic torque on the star from the disc is that the disc is flat.This is a nontrivial assumption.Indeed, the magnetic coupling between the star and the disc operates only in the innermost disc region (e.g., between the inner radius rin and rint ≈ 1.5rin), and this region has a much smaller moment of inertia than the star.Therefore, if there were no coupling between this inner disc region and the outer disc, the inner disc would be significantly warped on a timescale much shorter than the timescale for changing the stellar spin (Pfeiffer & Lai 2004).If there is any secular change in the stellar spin direction, the inner disc warp would then follow the varying spin axis.Clearly, in order to determine the long-term spin evolution of the star, it is important to understand the dynamics of the warped disc, taking into account the magnetic torques on the inner disc and the hydrodynamical coupling between different disc regions.This is the goal of our paper.
To be more specific, there is a hierarchy of timescales related to the combined evolution of the stellar spin and the disc warp: (i) The dynamical time t dyn associated with the spin frequency ωs, disc rotation frequency Ω and the beat frequency |ωs − Ω|.This is much shorter than the effects (steady-state disc warping and spin evolution) we study in this paper.
(ii) The warping/precession timescale of the inner disc [see Eq. ( 7)] where M , R , B are the mass, radius and surface (dipole) magnetic field of the protostar, respectively, θ is the inclination angle of the stellar dipole relative to the spin, Σ is the disc surface density, and ζ is a dimensionless magnetic twist parameter of order unity related to the strength of the azimuthal magnetic field generated by star-disc twist.
(iii) The disc warp evolution timescale t disc .This is the time for the disc to reach a steady-state under the combined effects of magnetic torques and internal fluid stresses (see Section 5).For high-viscosity discs, t disc is the viscous diffusion time for the disc warp [see Eq. ( 34)] and depends on the viscosity parameter α and the disc thickness δ = H/r: tvis ∼ (3000 yrs) α 0.1 For low-viscosity discs (α < ∼ δ), t disc is the propagation time of bending waves across the whole disc and depends on the sound speed.In general, t disc can be several orders of magnitude larger than tw.
(iv) The stellar spin evolution timescale.The magnetic misalignment torque on the star is of order µ 2 /r 3 in (µ being the magnetic dipole moment of the star), which is comparable to the fiducial accretion torque, given for Keplerian discs by N0 = Ṁ √ GM rin.Assuming the spin angular momentum Js = 0.2M R 2 ωs (the value for a Γ = 5/3 polytrope, representing a convective star), we find the spin evolution time In general tspin t disc .In this paper we will study the evolution of the disc warp on timescales ranging from tw to t disc , and the evolution of the stellar spin direction on timescales of order tspin.
It is important to note that we are not interested in disc warpings that vary on the dynamical timescale t dyn in this paper.In general, when the stellar dipole axis is inclined with respect to the spin axis, there will be periodic vertical forces at the rotation frequency of the star acting on the inner disc2 .These periodic forces will lead to the warping of the disc, particularly for low-viscosity discs in which bending waves propagate (Terquem & Papaloizou 2000;Lai & Zhang 2008).Indeed, there is observational evidence for such magnetically-warped discs.For example, the recurrent luminosity dips observed in the classical T Tauri star AA Tauri has been attributed to the periodic occultation of the central star by a warped inner disc (Bouvier et al. 2007).However, such dynamical disc warps average exactly to zero over a rotation period and have no effect on the secular evolution of the system.
The remainder of the paper is organized as follows.In Section 2, we summarize our analytical model of magnetopshere -disc interaction and derive the equation for the evolution of the stellar spin axis when the disc is warped.In Section 3 we present theoretical formalisms for determining the steady-state and time evolution of warped discs, for both high-viscosity regime (where warps propagate diffusively) and low-viscosity regime (where warps propagate as bending waves).An approximate analytical expression for the steady-state linear warp is also derived (see Section 3.2.2).In Section 4 we present numerical results for the steady-state disc warp profiles under various conditions and in Section 5 we study the time evolution of disc warps.We examine in Section 6 how the inner disc warp and the stellar spin evolution respond to variations of the outer disc, and discuss in Section 7 how this could, in principle, lead to anti-aligned planetary orbits, even for discs with initial angular momentum nearly aligned with the stellar spin.We conclude in Section 8 with a discussion of our results.

Magnetic Torques on the Disc
The interaction between a magnetic star and a disc is complex (see references in Paper I).However, the key physical effects of this interaction on the disc can be described robustly in a parametrized manner.The model used throughout this paper is detailed in Paper I. Here, we will limit ourselves to a brief summary of the magnetic torques acting on the disc.The stellar magnetic field disrupts the accretion disc at the magnetospheric boundary, where the magnetic and plasma stresses balance.For a dipolar magnetic field with magnetic moment µ, we have where η is a dimensionless constant somewhat less than unity (η ∼ 0.5 according to recent numerical simulations; see Long et al. 2005 3 ).We take rin to be the inner edge of the disc.Before being disrupted, the disc generally experiences nontrivial magnetic torques from the star (Lai 1999;Paper I).Consider a cylindrical coordinate system (r, φ, z), with the vertical axis Oz orthogonal to the plane of the disc.The magnetic torques are of two types: (i) A warping torque Nw which acts in a small interaction region rin < r < rint, where some of the stellar field lines are linked to the disc in a quasi-cyclic fashion (involving field inflation and reconnection).These field lines are twisted by the differential rotation between the star and the disc, generating a toroidal field ∆B φ = ∓ζB φ + ∆B φ differs above and below the disc plane, giving rise to a vertical force on the disc.While the mean force (averaging over the azimuthal direction) is zero, the uneven distribution of the force induces a net warping torque which tends to push the orientation of the disc angular momentum l away from the stellar spin axis ωs (see Paper I for a simple model for this effect, involving a metal plane in an external magnetic field).(ii) A precessional torque Np which arises from the screening of the azimuthal electric current induced in the highly conducting disc.This results in a difference in the radial component of the net magnetic field above and below the disc plane and therefore in a vertical force on the disc.The resulting torque tends to cause l to precess around ωs.In Paper I, we parametrized the two magnetic torques (per unit area) on the disc as where Σ(r) is the surface density, Ω(r) the rotation rate of the disc, and β(r) is the disc tilt angle (the angle between l(r) and the spin axis ωs).The warping rate and precession angular frequency at radius r are given by where θ is the angle between the magnetic dipole axis and the spin axis, and the dimensionless function D(r) is given by with H(r) the half-thickness of the disc.The function F (θ ) depends on the dielectric properties of the disc.We can write If the stellar vertical field is entirely screened out by the disc, the parameter f = 1; if only the time-varying component of that field is screened out, we get f = 0.In reality, f lies between 0 and 1.The magnetic torque formulae given above contain uncertain parameters (e.g., ζ, which parametrizes the amount of azimuthal twist of the magnetic field threading the disc); this is inevitable given the complicated nature of magnetic field -disc interactions.Also, while the expression for the warping torque [eq.( 5)] is formally valid for large disc warps, the expression for the precessional torque was derived under the assumption that the disc is locally flat [eq.( 8) is strictly valid only for a completely flat disc (Aly 1980)]; when this assumption breaks down (i.e., when |∂ l/∂ ln r| is large), we expect a similar torque expression to hold, but with modified numerical factors (e.g. the function D(r) in eq. ( 8) will be different).In the application discussed in the following sections, we find that the condition |∂ l/∂ ln r| < ∼ 1 is always satisfied.Thus we believe that our simple formulae capture the qualitative behavior of accretion discs subject to magnetic torques.
It is also worth noting that the expressions (5-6) for the torques only correspond to the zero-frequency component of the magnetic forces acting on the disc.The time varying components of these forces can also have significant effects.In particular, Lai & Zhang (2008) discussed how the components of the magnetic forces varying at the stellar spin frequency and at twice that frequency can excite bending waves in discs, while Terquem & Papaloizou (2000) showed that if the star has a dipole field misaligned with its rotation axis, magnetic effects create a steady-state warp in a frame corotating with the star.However, these "dynamical waves" average to zero over the stellar rotation period and do not affect the secular evolution of the stellar spin.In this paper, we concern ourselves only with long-term effects, effectively studying a disc profile averaged over multiple stellar rotations.

Spin Evolution of the Star
The effects of the magnetic torques on the evolution of the star -disc system are twofold.First, they will cause the orientation of the disc l(r) to deviate from a flat disc profile l(r) = lout = l(rout), set at the outer disc radius rout.These deviations will be studied in details for different disc parameters in Sections 3-5.Second, the back-reaction of the torques will change the orientation of the stellar spin axis on a longer timescale.The secular evolution of the stellar spin under the combined effects of matter accretion and stardisc interactions is explored in Paper I in the case of flat discs.Here we generalize the basic formulae derived in Paper I to warped discs.
In general, the spin angular momentum of the star, Js ωs, evolves according to the equation Here N l represents the torque component that is aligned with the inner disc axis l(rin) = lin.We parametrize N l by Equation ( 12) includes not only the accretion torque carried by the accreting gas onto the star, Ṁacc(GM rin) 1/2 l (where Ṁacc may be smaller than Ṁ , the disc accretion rate), but also the magnetic braking torque associated with the disc -star linkage, as well as any angular momentum carried away by the wind from the magnetosphere boundary (Shu et al. 1994;Romanova et al. 2009).All these effects are parametrized by the parameter λ < ∼ 1.In particular, if a wind carries away most of the angular momentum of the inner disc, we may get λ 1.The term N s = −|Ns| ωs represents a spindown torque carried by a wind/jet from the open field lines region of the star (e.g.Matt & Pudritz 2005).The terms N w and N p represent the back-reactions of the warping and precessional torques: Since both Nw and Np decrease rapidly with radius (as r −5 ), the integral can be carried out approximately, giving with where cos βin = ωs • lin.Note that both N 0np and N 0nw are of order µ 2 /r 3 in .For a fixed outer disc orientation lout, the inclination angle of the stellar spin relative to the outer disc, β = βout, evolves according to the equation Note that this does not depend on the specific form of N s.
For flat discs, equation (17) reduces to (Paper I) In the flat-disc approximation, the star -disc systems can thus be divided in two classes with very different long-term spin evolution (see Fig. 1).If ζ < λ, cos β always increases in time and the system will be driven towards the aligned state (β = 0).On the other hand, if ζ > λ, there are two "equilibrium" misalignment angles (defined by dβ /dt = 0): The smaller angle β+ corresponds to a stable equilibrium, while β− is unstable.Thus, the final state of the systems depends on the initial misalignment angle β (t = 0).If β (t = 0) < β−, the system will be driven towards a moderate misalignment β+ < 90 • ; otherwise it will evolve towards a completely anti-aligned configuration (β = 180 • ).From these results, we can see that, according to the flat-disc approximation, if λ 1 a misaligned configuration is strongly favored.
The probability distribution of the different cases for astrophysical systems will thus depend on the unknown value of the parameters of our model, as well as on β (t = 0)which depends on the formation history of the star -disc system and is quite uncertain (Bate et al. 2010).For example, for an isotropic distribution of lout on the unit sphere and ζ > λ, a fraction 0.5(1 − λ ζ ) of the systems would be anti-aligned, while the rest would tend towards a misalignment β+.The real distribution of disc inclination is certainly more complex, as λ and ζ will vary from system to system, and the distribution of the initial misalignment β (t = 0) is probably not isotropic.Additionally, the orientation of the outer disc might vary in time.For more details on the distribution of final inclination angles β , see Paper I (Sec.5), as well as Section 7 of this paper, which discusses a process to reach anti-alignment starting from β (t = 0) < β−.
In general, the magnetic torques induce disc warping so that l depends on r, and equation ( 17) must be used to determine the long-term spin evolution.Since the disc warp evolution timescale is much shorter than the stellar spin evolution timescale, the steady-state warp profile l(r) must be solved before equation ( 17) can be applied.In the following sections, we will show that for most (but not all) realistic choices of the free parameters in our model, using equation ( 18) instead of equation ( 17) does not significantly change our qualitative description of the long term behavior of the system.We will measure deviations from the flat-disc approximation through the parameter ξ defined by where the left-hand side is computed using the first line of equation ( 17).

DESCRIPTION OF WARPED DISCS: THEORY
As noted in Section 1, the evolution of the coupled star -disc system occurs over two different timescales.The first, t disc , characterizes the evolution of the disc under the magnetic torques and internal stresses towards a warped steady-state configuration, assuming that the spin of the star ωs is fixed.
The second, tspin, determines the evolution of ωs due to the combined effects of mass accretion and magnetic torques.
Since we expect t disc tspin, if the orientation of the outer disc is fixed we can consider that, at all times, the disc is in a steady-state leq(r; ωs).The evolution of the system is then described by a sequence of steady-state profiles leq(r, ωs(t)) where ωs(t) evolves according to equation ( 17) applied to l(r) = leq(r; ωs(t)).As discussed before, the disc itself will always show variations on shorter timescales (of the order of the stellar rotation period), which do not affect the secular evolution of the stellar spin and are averaged over in our description of the system.
Here we describe our method to calculate the evolution and steady-state of warped discs.
Systematic theoretical study on warped discs began with the work of Papaloizou & Pringle (1983) and Papaloizou & Lin (1995), who showed that there are two dynamical regimes for warp propagation in linear theory (for sufficiently small warps).For high viscosity Keplerian discs with α > ∼ δ ≡ H/r (where α is the Shakura-Sunyaev parameter so that the viscosity is ν = αH 2 Ω), the warp satisfies a diffusion-type equation with diffusion coefficient ν2 = ν/(2α 2 ).For low-viscosity discs, on the other hand, the warp satisfies a wave-like equation and propagates with speed ΩH/2.In the diffusive regime, the linear theory of Papaloizou & Pringle (1983) was generalized to large inclination angles by Pringle (1992) in the limit of small local variations of the disc inclination.A fully nonlinear theory was derived by Ogilvie (1999), with prescriptions for arbitrary variations of the inclination.The basic features of the theory were recently confirmed by the numerical simulations of Lodato & Price (2010).For low-viscosity Keplerian discs (α < ∼ H/r), the linearized equations for long wavelength bending waves were derived by Lubow and Ogilvie (2000) and Lubow et al. (2002), and a theory for non-linear bending waves was developed by Ogilvie (2006).
For protostellar discs, recent work by Terquem (2008) suggests that far away from the star the disc could have a very small viscosity parameter (α ∼ 10 −2 −10 −4 ), and would thus be described by the formalism of Lubow and Ogilvie (2000).However, close to the star (around a few stellar radii) where magnetic effects are most important and the disc warp can develop, the value of the effective viscosity is unknown.Thus in this paper, we will study both high-viscosity discs and low-viscosity discs.

Evolution Equations
For viscous discs satisfying α > ∼ H/r, we start from the equations derived by Ogilvie (1999).The main evolution equations for the disc are the conservation of mass and angular momentum where VR is the average radial velocity of the fluid at a given radius.The coefficients Q1,2,3 characterize the magnitude of the various viscous interactions, while depends on the vertical density profile of the disc.The term Nm = Nw + Np is the external magnetic torque per unit area.
In general, the viscous coefficients Q1,2,3 are functions of the viscosity parameter α, the warp amplitude ψ 2 ≡ |∂ l/∂ ln r| 2 , and the disc rotation law Ω.Their values can be obtained through numerical integration of a set of coupled ODEs (Ogilvie 1999).In the limit ψ 2 → 0, the viscous coefficients are given by equations [141][142][143] of Ogilvie (1999): For ψ 2 = 0 and Q3 = 0, this is equivalent to the formalism of Pringle (1992): the viscosities ν1 = ν and ν2 used by Pringle (1992), which correspond respectively to the shear viscosity usually associated with flat discs and the viscous torque working against the warping of the disc, are proportional to Q1 and Q2.The additional term Q3 was discovered by Ogilvie (1999), and contributes to the precession of a warped disc.For the disc configurations considered in this paper, the effects of finite ψ 2 are small -hence, our numerical results will be computed in the limit ψ 2 1.However, we do consider the effects of non-zero Q3.

Disc Model
For our numerical calculations, we consider Keplerian discs.
If we compare the projection of ( 23) along l for l = 0 with the standard flat-disc equation we see that (32) Using Q1[ψ 2 = 0] = −3α/2 and Ω = √ GM r −3 , we then have We also rescale the time coordinate by the viscous time evaluated on the inner edge of the disc: τ = t/tvis(rin), where and is the viscosity associated with the vertical shear in the disc.By projecting the evolution equation of the disc angular momentum onto directions parallel and orthogonal to l, we find where the new variables σ, S and ρ are defined by and Σ flat is the surface density of a flat disc Equations ( 36) and (37) form our model for the evolution of viscous discs interacting with a magnetic star.Note that as l is a unit vector, it only corresponds to two degrees of freedom in the system.Accordingly, equation (37) guarantees that ∂ l/∂τ is orthogonal to l.In practice, to avoid introducing a preferred direction in the system (as we want to allow arbitrary inclination angles for the disc), we evolve all 3 components of l, but normalize l at each timestep to the accumulation of numerical errors.

Steady-State Equations
From equations ( 36), ( 37) and (39), it is fairly easy to derive the equations defining the steady-state configuration of the disc.If we set ∂ l/∂τ = 0 and ∂σ/∂τ = 0 = S • l in (37), we obtain S = ( ωs. l) (42) Equation ( 39) projected onto l gives and projected in the plane orthogonal to l gives For Q3 = 0, we thus have a set of first order differential equations of the form U = F(U).Given appropriate boundary conditions at rin, it can easily be solved by numerical integration.For Q3 = 0, we can still perform numerical integration if we consider (44) as an implicit equation for l which has to be solved at each step of the integration algorithm.
In practice however, the boundary conditions are imposed partly at the inner edge rin and partly at the outer edge rout.Indeed, we consider the orientation of the outer disc to be fixed l(rout) = lout (45) and the mass accretion rate to be constant (the sign is chosen so that Ṁ > 0 for VR < 0).We also impose a zero-torque boundary condition at the inner edge l (rin) = 0 (47) and set the surface density there to for some freely specifiable scalar σin.Combining ( 46) with the zero-torque boundary condition and equations ( 28) and ( 41), we obtain a simple boundary condition on σ at rin: while (48) gives the value of σ[rin]: We thus have 4 boundary conditions at rin (on l , σ and σ ) and 2 at rout (on l).To solve the system numerically we use a shooting method starting at rin.Writing l = (cos β cos γ, cos β sin γ, sin β) and ωs = (0, 0, 1), we use a 2-D Newton-Raphson method to solve for the values of β[rin] and γ[rin] leading to a solution satisfying l(rout) = lout.The system of first-order ODEs which has to be solved at each iteration of the Newton-Raphson algorithm is treated using the 5th order StepperDopr5 method of Press et al. (2007), and the integration is performed under the constraint | l| = 1.

Evolution Equations
For discs with a viscosity parameter small compared to the thickness (α < ∼ δ = H/r), we can no longer use the evolution equations of Ogilvie (1999).In this case, disc warps propagate as bending waves.In the linear regime, the warp evolution equations were derived by Lubow and Ogilvie (2000): where cs = HΩz is the disc sound speed, Ωr and Ωz are the radial epicyclic frequency and the vertical oscillation frequency associated with circular orbits at a given radius from the star, G is the internal torque of the disc, and Σ = Σ flat is the surface density.These equations are only valid for α < ∼ δ, In the following, we shall use Ωr = Ωz = Ω, although we verified that small deviations from these equalities do not significantly modify our results.
Equations ( 51)-( 52) admit wave solutions.We define a Cartesian coordinate system so that lz 1 and | lx,y| 1, and the internal torque G acts in the xy-plane.Consider a local (WKB) wave with lx,y, G ∝ e ikr−iωt in a Keplerian disc with Nm = 0.For ω Ω, the dispersion relation of the wave is, neglecting the damping term αΩG in (52), with the eigenmodes satisfying lx,y = ( lx,y The + mode and − mode correspond to the outgoing and ingoing bending waves, respectively.A generic warp perturbation will not behave as pure eigenmodes.For numerical evolutions, it is convenient to define the variables Then, the evolution equations for the disc can be written as Here the dimensionless time τ = tδΩ(rin) and length ρ = r/rin are chosen so that the sound speed at the inner edge of the disc is cs(rin) = HΩz = 1.For the computation of the magnetic torque, we use the following approximations, accurate to first order in lx,y: The boundary conditions are particularly simple to implement for this choice of variables.At the outer edge of the disc, we require the ingoing mode to vanish At the inner edge, we impose the zero-torque boundary condition l = 0, G = 0, which can be written in terms of our evolution variables as In terms of the propagation of bending waves, this corresponds to the requirement that the waves be reflected at the inner edge of the disc.

Steady-State Warp
The steady-state profile of low-viscosity discs can be obtained by numerical integration of equation ( 56) or equations ( 51)-( 52), by setting ∂/∂τ = 0.In practice however, the steady-state profile of a low-viscosity disc is nearly always very well approximated by a flat disc profile.The amount of disc warping can then be evaluated analytically.Noting that Σr 3 c 2 s ∝ r 3/2 , equations ( 51)-( 52) can be combined to give Since Nm is falling rapidly with r (Nm ∼ r −5 ), and ∂ l/∂r = 0 at r = rin, we integrate the above equation from rin to r: Integrating from rout to rin, we then obtain Using Eqs. ( 5)-( 6), we have where tvis = r 2 in /ν2 is the viscous timescale for the warp, with ν2 = csH/(2α).Thus, the distortion of the disk can be seen as arising from the warping and precessional torques acting over the disc during a time of order the viscous time scale at the inner disc edge (where the magnetic torques are the strongest).Projecting Eq. ( 64) in the direction of the stellar spin axis ωs and using Eqs.( 5)-( 6), we have we see that as long as α 2 ζ η 7/2 , a condition satisfied for most parameters, the warp across the whole disc is small: For example, with η > ∼ 0.5, we find that for all discs |βin − βout| 1 if α < ∼ 0.15.This is almost certainly true for discs in which bending waves can propagate.
It is important to note that, although the approximate analytical expression of the global disc distortion derived above is based on low-viscosity discs, our result for |βin − βout| is also valid for higher-viscosity discs.Indeed, in the linear regime and for Keplerian discs, the steady-state equations are identical regardless of the viscosity regime considered.

STEADY-STATE PROFILE OF WARPED DISCS AND BACK-REACTION ON STELLAR SPINS
Using the numerical scheme presented in Section 3.1.3,we can now determine the time-averaged steady-state profile of the disc under the influence of the torques exerted by a magnetic star.The characteristics of the warped disc will of course vary with the choice of the free parameters included in our theoretical model.We begin our study by showing results for two standard discs, chosen so that they belong to the two classes of long term stellar spin evolution predicted in Section 2.2 when the accretion parameter defined in equation ( 12) is λ ≈ 0.5 (a typical value in the allowed range 0 λ 1).We then vary the disc parameters, and discuss their influence on the disc profile, and on the spin evolution.Finally, we check that, as predicted in section 3.2.2,low-viscosity discs, which follow the different evolution equations described in section 3.2.1 (valid for α < ∼ δ = H/R) have only negligible steady-state warps and are for all practical purposes well described by the flat-disc approximation.
Our base models are discs with viscosity α = 0.15 and thickness δ = 0.1.We fix the surface density at the inner boundary by setting σin = 1.0 so that it is equal to the surface density of a flat disc (41), choose the inclination angle of the outer disc βout = β = 10 • and the magnetic inclination angle θ = 30 • with respect to the spin ωs.The star is assumed to have mass M = M and radius R = 2R The strength of the magnetic field is chosen so that rin = 2.5R * (Corresponding to B ∼ 1kG for typical parameters, see Eq. [4]), and the action of the torque Nm is limited to the region rin r rint = 1.5rin.The accretion rate is Ṁ = 10 −8 M /yr, and we put the outer disc boundary at rout = 10 4 rin (corresponding to rout ≈ 250AU, a size typical of the observed protoplanetary discs).This disc has small values of ψ 2 everywhere, and accordingly we neglect the nonlinear terms in Qi (but we keep Q3 = 3/8).The other parameters are chosen to be ζ = 1, f = 0, and either η = 1 (so that the long-term evolution of the system aligns the spin axis ωs with the disc axis) or η = 0.5 (for which the flat-disc approximation predicts a long term misalignment toward β+ ≈ 45 • if the initial disc has β 135 • ).We should note that these parameters are purposefully chosen to test the limits of the flat disc approximation.Our choices M , R , B , δ and M are relatively standard values for protoplanetary discs around T-Tauri stars [see Bouvier et al. (2007b) and references therein], while there are no particular reasons to prefer any specific orientation of the magnetic dipole θ .But α = 0.15 is larger that recent estimates of the viscosity in the outer parts of the disc (Terquem 2008), and probably on the high end of what can be expected in the inner disc.However, we have shown that smaller values of α lead to smaller amplitudes of the steady-state warp (the warp amplitude is proportional to α 2 ).Thus, the flat disc approximation is more likely to be satisfied at low viscosities.
In order to analyze the radial variations of the disc warp profile, we define the tilt β[r] and the twist γ[r] by l with the convention that γ[rout] = 0. Some parameters of the system can be varied without modifying the dimensionless solution for the profile of the surface density σ(ρ) and the orientation of the disc l(ρ): modifications of Ṁ , M , R or rin (at constant η, rout/rin and rint/rin) will influence the values of the timescales tvis and tspin, but not Thus, the steady-state profile can be solved while keeping these parameters fixed without any loss of generality.The disc profile in physical units [Σ(r), l(r)] can easily be retrieved from the dimensionless solution [σ(ρ), l(ρ)].Additionally, the four parameters (η, θ , f, ζ) correspond to only two degrees of freedom in the model, through the quantities We will thus limit ourselves to variations of ζ and f .Varying the thickness δ of the disc has very similar effects: it changes the value of the function D(r) at small radii, effectively modifying the value of c2 close to rin.As the magnetic torques mostly affect the region close to the inner edge of the disc, the influence of δ is similar to that of F (θ ).Finally, we are also free to modify the boundary conditions used, and in particular the choices of rout and σin.Varying rout seems to have only negligible effects, as long as rout/rin is large enough for a steady-state solution to exist.Decreasing σin, on the other hand, leads to more significant changes in the warp profile.A small σin favors warping disc, so that a decrease of σin has an effect similar to increasing both c1 and c2.
Thus, the effects of varying various parameters of the system can be examined with our standard discs, by varying only two parameters, ζ (or c1) and f (or c2).In section 4.1, we present our results for our two standard discs, which are similar except for the value of the parameter η (changing η correspond to a rescaling of both the warping and the precessional torque).Then, in section 4.2, we study variations of the warping torque alone, by modifying the value of the parameter ζ characterizing the strength of the toroidal field in the disc.The influence of the precessional torque is studied in more details in section 4.3, through variations of the parameter f (related to the ability of the time-varying component of the vertical magnetic field to penetrate the disc).Finally, in section 4.4 we comment on the influence of the parameter Q3, which was usually neglected in previous studies of warped discs.

Standard disc results
The profile for the tilt and twist angles of our standard configurations (Fig. 2) show a relatively weak warping of the disc.For the disc with weaker magnetic interactions (η = 1), the difference in tilt between the inner and outer edges is about 0.17 • and the twist over the whole disc is 1.2 • , while for stronger interactions (η = 0.5) the disc is tilted by 1.8 • and twisted over 16 • .These warps are comparable in magnitude to what we could have predicted using the approximate equations ( 64) and (68).In particular, formula (68) applied to these two choices of parameters predicts tilts of 0.28 • and 3.2 • , respectively, with most of the difference between the approximate formula and the numerical results due to the cutoff applied to the magnetic torques at r = rint, neglected in the derivation of (68).
Note that we choose to vary the parameter η defined in equation ( 4) as it conveniently modifies the effective strength of both magnetic torques in our model.In practice, η is determined by the geometry of the accretion flow, while unknown physical parameters such as the dipole strength µ, its orientation θ , the surface density at the inner edge of the disc σin or the magnetic twist parameter ζ will vary from system to system.
Given the small warp, the evolution of the misalignment angle β between ωs and lout is well approximated by equation ( 17) with lin = lout: if we compute N from the steady-state profile l(r), we find that the parameter ξ in equation ( 21), which parametrizes deviations from the flat disc approximation (ξ = 1 for a flat disc) is ξ = 0.997 for η = 1 and ξ = 0.86 for η = 0.5, if we set the accretion parameter λ to 0 (we choose λ = 0 when computing ξ in order to measure directly differences in the effect of the backreaction magnetic torques between the flat-disc model and the warped disc steady-state, disentangled from the effect of angular momentum accretion).
However, even a small disc warp can significantly change the critical angles β± for which dβ /dt = 0.By varying β = β(rout), we can determine the values of β± numerically.For η = 0.5 and λ = 0.5, we find β+ = 32 • and β− = 148 • which are quite different from the prediction of the flat-disc approximation (β+ = 45 • , β− = 135 • ).This is due mainly due to the effect of the twist of the disc.In the flat-disc approximation, γ = 0 and the back-reaction due to the precession torque has no effect on the evolution of the stellar spin.But as long as γin < π, that back-reaction will tend to align the stellar spin and the disc orbital angular momentum, and this effect can be large enough to significantly shift the value of β± (see also subsection 4.3).The same effect can also modify the qualitative behavior of systems for which the predicted misalignment angle β+ is close to 0, in such a way that the stable configuration at β+ no longer exist.The orbital angular momentum of the disc would then be expected to align with the direction of the stellar spin.

Large warping torques
Realistic discs are expected to have parameters ζ ∼ 1 (characterizing the azimuthal magnetic twist) and η ∼ 0.5 (characterizing the inner disc radius; see Eq. [4]), but the exact values of those parameters are unknown (see Section 2).By increasing ζ we see that the approximation l(rin) = l(rout) can break down for high-viscosity discs.In Fig. 3, we show the variation of the steady-state disc profile when ζ is varied between 1 and 5 for η = 0.5 and our standard disc parameters.Clearly, there can be large differences between the orientation of the disc at its inner and outer edges when ζ > ∼ 1.In Fig. 4 influence of the precessional torque (which is independent of ζ) when the warping torque becomes small.The long-term evolution of the stellar spin direction will remain similar to the flat-disc predictions, with a stable configuration at some misalignment angle β+ = 0 for most values of the accretion parameter λ.But for ζ > ∼ 4.5, the twist is so large that the behavior is the opposite of what would be predicted by our approximate flat-disc formula: the back-reaction tends to align the disc and the spin of the star.This shows that at large ζ we must determine for each set of parameters the profiles β[r] and γ[r] in order to predict the long term evolution of the stellar spin.As an example, we construct sequences of steady-state disc configurations for a fixed ζ, varying the inclination angle of the outer disc β .In Figs.
5-7, we show the resulting d cos β /dt for ζ = 1, 3 and 5, and compare with the predictions of the flat-disc approximation.
For ζ = 1, 3, the general behavior is similar to what the flatdisc approximation predicts.As seen in subsection 4.1, the precessional torque will favor alignment of the stellar spin with the disc orbital angular momentum, so that the numerical results usually show that β+,Num β +,F lat -at least as long as F (θ ) is of order unity.For ζ = 5, however, significant differences become visible.At small inclination angles, the system will evolve towards β = 0, while at large inclinations, the system will evolve towards β ≈ 165 • .In the intermediate region 15 • < ∼ β < ∼ 135 • , the system will evolve towards β ≈ 50 • (for λ = 0.5).Finally, for larger ζ we are in a completely different regime: for some inclinations, two steady-state solutions exist.Clearly, to determine which of those steady-state solution is relevant requires numerical integration of the time evolution of the star-disc system.
Our current understanding of the effects of magnetic fields close to the inner edge of the accretion disc is not sufficient to determine with certainty the range of realistic values of the parameters ζ and η.However, their favored values lie in a region of parameter space where the flat disc approximation appears to hold relatively well (ζ ∼ 1, η ∼ 0.5).The large deviations from the flat disc model observed at high ζ are thus unlikely to be encountered in astrophysical systems, though they cannot be entirely ruled out.Thus, these results implies that the flat disc approximation is likely to be justified, with the caveat that it tends to overestimate the value of the misalignment angle β+.

Varying the precessional torque
The results presented in previous subsections were all obtained with f = 0, thus fixing the choice of the function F (θ * ) characterizing the magnetically driven disc precession rate.Using different values of f , even at low ζ it is possible to find discs which require numerical solutions to determine their warp profiles.For example, if we choose f = 1 instead of f = 0, the sign and magnitude of Ωp will change.
The twist of the disc becomes more important, so that even for ζ = 1, η = 0.5, there is a significant deviation from the behavior of the flat-disc configuration.Comparisons between the disc profiles for different f can be found in Fig. 8.The most important feature of these profiles is that, for the larger values of f , we have a large twist γ(rin).Hence, the precession term in equation ( 17) (proportional to np), which does not contribute to the evolution of β in the flatdisc approximation, now has a significant impact.For a twist γ such that sin(γ − γ[rout])F (θ ) > 0, the precession term directly contributes to the alignment of the outer disc axis with the stellar spin.This is always the case for disc twists |γin| 180 • , as a positive F (θ ) causes the inner disc to precess in the prograde direction, while F (θ ) 0 causes a retrograde precession.If the precessional torque becomes large enough compared to the warping torque (proportional to nw), the long-term evolution of the stellar spin direction will be modified.For our standard parameters and the choices of η = 0.5 and λ = 0.5, we find that discs with f > ∼ 0.5 will always lead to spin-disc alignment, contradicting the flat disc predictions (see Fig. 9).However, it is worth noting that some configurations with high f still allow for long term misalignments: for example, for f = 0.5, increasing the strength of the azimuthal B-field to ζ = 3 leads to a behavior very similar to what we found for f = 0, ζ = 3 (see Fig. 6), while decreasing the viscosity parameter to α = 0.015 (and choosing δ = 0.01) limits the twist of the disc, so that the flat-disc approximation remains valid.
The above results show that the precessional torque can in principle cause non-negligible deviations from the flatdisc model.Nevertheless, for the largest part of the favored parameter space (small α, or large α with small precessional torque), the flat-disc approximation is justfied.

Influence of Q3
As mentioned before, most previous works on warped discs have been done using the formalism of Pringle (1992), which corresponds to Q3 = 0 in the formalism of Ogilvie (1999).This is a good approximation, as long as the influence of the small precessional torque due to Q3 = 0 is negligible.For .Same as Fig. 5, except that we set λ = 0.5 and choose f = 0, 0.16, 0.5, 1.For f = 0.16, the disc twist is very small, and our numerical result matches the flat-disc approximation better than for f = 0.For larger values of f , the disc twist is large, and the disc will align with the stellar spin regardless of the initial value of β .
the system studied here a small change in the twist of the disc can affect whether a configuration will align over time, or be driven towards a stable misaligned steady-state.In Fig. 10, we show the difference in the disc tilt and twist our standard disc with η = 0.5, using both the formalism of Ogilvie (1999) and Pringle (1992).Differences in the warp of the disc of a few degrees are observed, though the warps are small in both cases.Because the precessional torque acting on the disc using Q3 = 0 is smaller, it will be less twisted.This leads to a behavior slightly closer to what the flat-disc approximation predicts.If we choose the accretion parameter λ = 0.5 (equation 12), then the flat-disc approximation predicts a stable misaligned configuration at β+ = 45 • .For warp discs, we find that the misalignment angle is significantly smaller, β+ = 32 • .The difference in β+ between profiles obtained using Q3 = 3/8 and Q3 = 0 is only 0.5 • , which is negligible at the level of accuracy our model can achieve.
Q3 can also influence the qualitative behavior of the steady-state solutions at high-ζ.For strongly warped discs, it is sometimes possible to have two solutions satisfying the steady-state equations.Choosing Q3 = 0 seems to limit the size of the region of parameter space where this happens.For example, for f = 0, η = 0.5, θ = 10 • and ζ = 5.5, two profiles are acceptable steady-state solutions if we chose Q3 = 0, while for Q3 = 3/8 the same parameters lead to a unique solution (see Fig 11).

Low-Viscosity Discs
In the linear regime, the equations determining the steadystate profile of the disc are identical for the α δ and α δ cases.In the previous subsections, we have seen that our approximate formulae for the amplitude of the warp, equations ( 64) and ( 68), give relatively good results for α ∼ δ = 0.1.We also confirmed numerically the α 2 dependence of the warp of the disc, shown in Figure 12.For smaller viscosities, α δ, we expect the warp to be even smaller, and the linear approximation more accurate.Hence, we can immediately deduce that the time-averaged warp of low-viscosity discs will be extremely small.For such discs, the flat-disc approximation will nearly always give accurate results for the secular evolution of the stellar spin.

TIME EVOLUTION OF DISC WARP TOWARD STEADY-STATE
Having established the steady-state of warped discs, we now study their time evolution starting from some generic initial conditions, when the symmetry axis of the outer disc is misaligned with the stellar spin.To this end, we evolve equations ( 36)-( 37 is necessary.Our numerical method is detailed in the Appendix.

High-Viscosity Discs
For viscous discs with α > ∼ δ = H/r, we expect the evolution of the system to occur over the timescale tvis(r) = r 2 /ν2 = 2α/(δ 2 Ω).Note that at the disc inner edge, tvis(rin) = (3α 2 ζ cos 2 θ /2η 3.5 )Γ −1 w (rin) [see Eq. ( 67)] is smaller than the warping timescale for typical parameters.In terms of the dimensionless time τ = t/tvis(rin), we expect the disc to reach the steady-state profile at radius r within a time of order τ ∼ (r/rin) 3/2 (assuming constant δ).To test this expectation, we evolve our standard disc model (see Section 4) for η = 0.5 and different locations of the outer radius (rout = 100rin and rout = 1000rin), as well as for a more viscous disc with α = 0.3.The disc is initialized in a flat configuration with l = lout and we observe its evolution towards the steady-state profile.In Figs.13-15, we plot the disc warp profiles at times τ = 10 3n/4 for n = 0, 1, ..., 4by which point the viscous forces should have brought the disc into its steady state up to radius r ∼ 10 n/2 rin.
In all cases, we see that the evolution occurs approximately on the expected timescales: the local distortion of the disc (i.e.∂ l/∂ ln r) up to radius r does not vary much past the viscous timescale at that radius.The orientation of the disc ( l), on the other hand, continues to change to accommodate the evolution of the disc at larger radii.Overall, the disc will reach its equilibrium profile within the viscous timescale tvis(rwarp), where rwarp is defined as the largest radius at which the warp |∂ l/∂ ln r| is significant.
For the two simulations with the outer disc boundary at rout = 100rin, we find that rwarp ∼ rout and the disc reaches its steady-state profile within τ ∼ 1000.At later times, the evolution of the profiles becomes negligible.For the larger disc (rout = 1000rin), the situation is slightly different.At τ = 1000, the disc has reached its steady-state distortion up to r = 100rin.The disc will still evolve up to τ ∼ 10 4.5 , but as the warp is very small for r > ∼ 100rin, the changes in the profile are minimal.As most discs studied in this paper show negligible warps for r > (10 2 − 10 3 )rin, we expect the steady-state to be reached within at most tvis(10 3 rin) ∼ 10 4.5 r 2 ν2 in ∼ 500 α 0.15 δ 0.1 −2 yrs, (72) regardless of the outer radius of the disc.As this is much smaller than the evolution timescale for the spin of the star, we are justified to consider only the steady-state configuration of the disc when attempting to determine the long-term evolution of the misalignment between the stellar spin and the orientation of the outer disc.

Low-Viscosity Discs
The evolution of low-viscosity discs (α < ∼ δ) is qualitatively different from high-viscosity discs.According to equation (56), perturbations around the steady-state propagates as bending waves, at roughly half the local sound speed.Thus, we expect the disc to settle to an equilibrium within the propagation timescale of these waves, In Figures 16-18, we show the evolution of the disc tilt profile β as the bending wave propagates across the disc, using our standard choice of parameters for the magnetic torques.We consider different discs: the first two use α = 0.01 and have their outer boundaries at rout = 100rin (Fig. 16) and rout = 1000rin (Fig. 17).The third has a higher viscosity α = 0.05, and rout = 1000rin (Fig. 18).All three simulations are started from a flat disc configuration, and show the same behavior: the magnetic torques perturb the inner disc, and the perturbation propagates outwards over the timescale twave.Again, this timescale is much less than the spin evolution timescale.
The location of the outer disc radius apparently does not have a significant influence on the final state of the system.The viscosity, on the other hand, affects the disc warp amplitude as predicted by equations ( 64) and ( 68): the warp is proportional to α 2 , with the amplitude |βout − βin| given by Eq. ( 68) to within a factor of two.

VARIATIONS OF THE OUTER DISC ORIENTATION
In the previous section, we have studied the time evolution of warped discs under the assumption that the orientation of the outer disc is fixed.However, a protoplanetary disc is formed inside the star forming core of a turbulent molecular cloud [e.g., McKee & Ostriker (2007)].Thus in general we expect the outer orientation of protoplanetary discs to have some variations in time.In this section, we study how the warped disc and particularly the inner disc orientation respond when the outer disc orientation varies by some finite amplitude (chosen to be 20 • ) over a period of time short compared to the evolution timescale of the disc, and how such variations affect the secular evolution of the stellar spin direction. .Time evolution of the disc tilt angle profile β for α = 0.15 and rout = 1000r in , when the outer disc orientation is changed from β(rout) = 10 • at t = t 0 − ∆t to β(rout) = 30 • at t = t 0 , with = 10 3 t vis (r in ).Time is in units of t vis (r in ).

High-Viscosity Discs
We first consider a viscous disc with α = 0.15 and rout = 1000rin.We choose to vary the outer disc orientation over ∆t = 1000tvis(rin) ∼ tvis(100rin).As in the case of the evolution towards the steady-state, the evolution of the disc occurs on the viscous timescale tvis (see Fig. 19).However, as significant changes now take place at the outer radius, the new steady-state configuration will be reached in a time of order the viscous timescale at the outer radius tvis(rout), whis is larger than tvis(rwarp) (see Section 5.1).Nevertheless, even though the steady-state is likely to be reached over a longer timescale than when the outer orientation is fixed, we still expect tvis(rout) to be significantly less than the evolution time for the stellar spin tspin.Thus, if the variation of the orientation of the outer disc occurs on a timescale shorter than tspin, the evolution of the stellar spin is well described by the approximation in which the disc is assumed to be in its steady-state configuration at all times, and adapting instantaneously to modifications of its orientation at the outer boundary.

Low-Viscosity Discs
The same type of evolution can also be studied for lowviscosity discs.If we choose the viscosity parameter α = 0.01, and change the orientation of the disc by 20 • over a timescale ∆t = twave(100rin), where twave(r) is defined by equation ( 73) with rout replaced by r, we obtain the evolution shown in Fig. 20.We see that a bending wave created at the outer boundary propagates inward, until it reaches the inner edge of the disc where it is reflected.The total time required for the disc to reach a new steady-state is thus twice the crossing time of the bending wave, ∼ 8/(3δΩ(rout)).For low-viscosity discs, the condition for the steady-state ap- proximation to be valid when the orientation of the outer disc is allowed to change over time is thus As tspin ≈ 1Myr [see Eq. (3)], this condition is easily satisfied.Also note that the evolution equations of bending waves adopted in our analysis are based on the flat-disc approximation.When the outer disc boundary evolves as fast as shown in Fig. 20, this approximation is no longer valid.Thus in practice, we should also require ∆t 8/[3δΩ(rout)].

APPLICATION: ANTI-ALIGNED EXOPLANETARY ORBITS
Our calculations in Sections 3-5 show that for the most likely physical parameters that characterize a magnetic star -disc system, the disc warp is small.Therefore the long-term evolution of the stellar spin is generally well-described by equation (18), as long as the orientation of the outer disc is kept constant.According to (18), three types of spin evolution trend are possible, depending on the parameters of the system and the initial conditions (Paper I).If ζ < λ, the stellar spin and the disc axis will always align (given enough time) regardless of their initial relative inclination.If ζ > λ, misalignment between the disc and the stellar spin will develop, evolving towards one of the two possible final states: either β = β+ < 90 • , or a perfectly anti-aligned configuration.
The second configuration can only be reached if the initial disc has a retrograde rotation with respect to the stellar spin, with β(t = 0) > 180 • − β+ = β−.In this case, to explain the observed expolanetary systems with retrograde orbits relative to the stellar spin (Triaud et al. 2010), we have to require that the disc rotates in a very different direction from the stellar rotation axis during the time of planet formation.
As discussed in paper I [called scenario (2) in Section 5 of Paper I], this is certainly possible if we consider the complex nature of star formation in molecular clouds and in star clusters [see also Bate et al. (2010)].
In Paper I, we describe another potential pathway to create retrograde exoplanetary systems [called scenario (1)] starting from prograde-rotating discs.If the disc axis and stellar spin axis are initially nearly (but not perfectly) aligned, and the magnetic torques are such that the aligned configuration is unstable, then the misalignment angle will tend towards β+, and no retrograde planets can be produced.However, this is only true if the orientation of the outer disc does not vary.If instead we assume that the outer disc experiences a change of its orientation ∆β > β− − β+ over a timescale ∆t sufficiently long that this change can propagate to the inner disc, but short enough that the stellar spin direction does not significantly evolve over ∆t, then the star -disc inclination can jump to β > β−, and continue to evolve towards anti-alignment.These conditions can be summarized as: where the disc warp evolution time t disc ∼ tvis(rout) if α > ∼ δ (high-viscosity disc) and t disc ∼ twave for α < ∼ δ (lowviscosity disc).As we have seen in Sections 6.1-6.2, the second and third conditions are fairly easy to satisfy, as t disc is at most of order 10 4 yrs for a viscous disc with rout ∼ 10 4 rin (and t disc would be significantly shorter for a smaller outer disc radius), while tspin ∼ 10 6 yrs for typical parameters [See Eq. ( 3)].The potential to satisfy the first condition, on the other hand, will depend on the fraction of the disc angular momentum which is accreted by the star (the parameter λ in equation 12) and the magnetic warp efficiency (the parameter ζ).If the star only accretes a small fraction of the angular momentum (λ 1), then the angles β± are both close to 90 • , and small variations of the outer disc are sufficient to allow the system to jump to the retrograde state and eventually evolve towards the anti-aligned configuration.

DISCUSSION
The main finding of our paper is that although magnetic interactions between a protostar and its disc have a strong tendency to induce warping in the inner disc region, internal stresses in the disc tend to suppress the warping under most circumstances.The result is that in steady-state, the whole protoplanetary disc approximately lies in a single plane, which is determined by the disc angular momentum at large radii (averaging out the dynamical warps which vary on timescales of order the stellar rotation periodsuch dynamical warps do not affect the secular evolution of the stellar spin).The reason for the small steady-state disc warp is that the effective viscosity acting to suppress disc warp, ν2 ν1/(2α 2 ), is much larger than the viscosity (ν1 = αHcs) responsible for angular momentum transfer within the disc (Papaloizou & Pringle 1983;Ogilvie 1999).In fact, our anaylsis of the steady-state magnetically driven disc warp shows that, in the linear regime, the disc inclination angle (relative to the stellar spin axis) varies from the outer disc to inner disc by the amount [see Eqs. ( 66 where tvis = r 2 /ν2 is the viscous time and Γw is the warping rate due to the magnetic torque.This result is valid regardless of whether the warp perturbations propagate diffusively (for α > ∼ H/r, high-viscosity discs) or as bending waves (for α < ∼ H/r, low-viscosity discs).Thus, for the preferred values of the parameters η ∼ 0.5, ζ ∼ 1, we find |βin − βout| 1 for α 0.3.Moreover, our analysis of the time evolution of warped discs shows that, starting from a generic initial condition, the steady-state can be reached quickly, on a timescale shorter than the characteristic timescale for the evolution of the stellar spin orientation.
Overall, our study of magnetically driven warped discs presented in this paper justifies the approximate analysis (based on the flat-disc approximation) of the long-term evolution of spin-disc misalignment presented in Paper I. Nevertheless, we note that even relatively small disc warps can modify the "equilibrium" spin -disc inclination angles β± (see Fig. 1) from the flat-disc values, thereby affecting the "attractors" of the long-term evolution of the spin -disc inclination angle.If we allow for more extreme parameters (but still reasonable by physical considerations) for the disc -star systems, much larger disc warps become possible and qualitatively different evolutionary trends for β may be produced (see . Taken together, the results of this paper and paper I demonstrate that at the end of the first stage of the planetary system formation (see Section 1), the inclination angle between the stellar spin and the angular momentum axis of the planetary orbit may have a wide range of values, including alignment and anti-alignment (see also section 7).Dynamical processes (e.g., planet-planet scatterings and Kozai interactions) in the second stage, if they exist, would further change the spin -orbit misalignment angle.More work is needed to determine the relative importance of the two stages in shaping the properties of planetary systems.Currently, the orbital eccentricity distribution of exoplanetary systems suggests that the second stage is important (e.g., Juric & Tremaine 2008).On the other hand, as noted in paper I, the 7 • misalignment between the ecliptic plane of the solar system and the sun's equatorial plane may be explained by the magnetically driven misalignment effect studied in this paper.Also, the recent discovery of Kepler-9 (Holman et al. 2010), a planetary system with two or three planets that lie in the same orbital plane, seems to suggest that at least some planetary systems are formed in a "quiet" manner without violent multi-body interactions.Obviously, measuring the stellar obliquity of such "quiet" systems would be most valuable.conditions.We encounter two types of boundary conditions: Dirichlet conditions of the type y = yBC are enforced by replacing the discretized version of (A1) by ỹ0,N−1 = yBC, while Neumann conditions of the type y = y BC are enforced by explicitly replacing y by y BC whenever necessary in (A1).If a second derivative is required to evaluate (A1) at the boundary, y−1 and yN are obtained using (yi+1 − yi−1) = (∆x)y i and the known value of y at the boundary.
, where ζ ∼ 1(Aly 1985;Lovelace et al. 1995) and the upper/lower sign refers to the value above/below the disc plane.Since the toroidal field from the stellar dipole B (µ) φ is the same on both sides of the disc plane, the net toroidal field B φ = B (µ)

Figure 1 .
Figure 1.Time derivative of the inclination angle of the stellar spin (in degrees) relative to the outer disc for a flat (unwarped) disc.The upper panel is for ζ/λ = 0.5 < 1, in which case the spin evolves towards alignment.The lower panel is for ζ/λ = √ 2 > 1, in which case the spin either evolves toward β + = 0 or toward β = 180 • , depending on the initial value of β .The quantity t spin is defined in equation (3).The arrows show the direction of the evolution of β in different regions of the parameter space.

Figure 3 .Figure 4 .
Figure 3. Upper panel: Disc tilt angle β for different choices of ζ values, all with η = 0.5.Lower panel: Twist angle γ for the same disc parameters.The outer edge of the disc is at rout = 10 4 r in .

Figure 5 .
Figure5.Secular evolution rate of the spin-disc inclination angle β for discs with ζ = 1.The time derivative of cos β is given for the flat-disc approximation (Flat) and for our numerical results for warped discs (Num), as well as for 3 different values of the accretion parameter λ = 0, 0.5, 1(see equation 12).The angle β will increase if d cos β /dt < 0.

Figure 7 .
Figure7.Same as Fig.5, except that we use ζ = 5.Note that when the disc is nearly aligned or nearly anti-aligned, the qualitative behavior of the solution is different from the predictions of the flat-disc approximation.

Figure 8 .
Figure 8. Upper panel: Disc tilt angle β for different choices of f .Lower panel: Twist angle γ for the same disc parameters.The outer edge of the disc is fixed at rout = 10 4 r in .

Figure 10 .
Figure 10.Upper panel: Disc tilt angle β for different choices of Q 3 .Lower panel: Twist angle γ for the same disc parameters.The outer edge of the disc is fixed at rout = 10 4 r in .

Figure 13 .
Figure13.Time evolution of the disc tilt angle profile for standard discs with α = 0.15 and rout = 100r in .Time is in units of t vis (r in ).

Figure 16 .
Figure 16.Evolution of the disc tilt angle profile for discs with α = 0.01 and rout = 100r in .The unit of time is (δΩ(r in )) −1 .
Figure19.Time evolution of the disc tilt angle profile β for α = 0.15 and rout = 1000r in , when the outer disc orientation is changed from β(rout) = 10 • at t = t 0 − ∆t to β(rout) = 30 • at t = t 0 , with = 10 3 t vis (r in ).Time is in units of t vis (r in ).

Figure 20 .
Figure 20.Same as Fig. 19 except for a disc with α = 0.01 and rout = 1000r in , and the outer disc orientation varies by 20 • over ∆t = twave(100r in ).