Issue 
A&A
Volume 531, July 2011



Article Number  A163  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201016314  
Published online  07 July 2011 
An investigation of magnetic field distortions in accretion discs around neutron stars
II. Analysis of the toroidal field component
^{1}
Key Laboratory of Solar Activity, National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, PR China
email: luca.naso@gmail.com
^{2}
SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy
^{3}
Department of Physics (Astrophysics), University of Oxford, Keble Road, Oxford OX1 3RH, UK
Received: 14 December 2010
Accepted: 4 April 2011
Millisecond pulsars are believed to be old pulsars spun up by a surrounding accretion disc. Magnetic fields are thought to play a leading role in this, both by determining the location of the inner edge of the disc and by exerting an additional torque on the star (as a result of the interaction between the stellar magnetic field and the disc plasma motion, which creates a toroidal component B_{φ}). In some wellknown analytic models, developed in the 1980s, the B_{φ} profile was taken to be proportional to the relative angular velocity between the disc plasma and the neutron star, multiplied by a vertical dipolar field. The present work stands in the line of improving those models, suggesting a new profile for B. In a previous paper, we discussed the poloidal component of the magnetic field and here we consider the toroidal component, again making the kinematic approximation and looking for steady solutions of the induction equation for axisymmetric models. The poloidal magnetic field is not assumed to be dipolar and the poloidal velocity field is not taken to be zero everywhere. We also do not use the thin disc approximation to simplify the induction equation but instead solve it numerically in full 2D. The profile obtained in the earlier analytic models is shown to have very limited validity and a more general semianalytic solution is proposed.
Key words: accretion, accretion disks / magnetic fields / magnetohydrodynamics (MHD) / turbulence / methods: numerical / Xrays: binaries
© ESO, 2011
1. Introduction
In the present work we study the deformation caused in a neutron star’s magnetic field because of the interaction with the matter in a surrounding accretion disc. A basic description for this kind of system was given by Ghosh and Lamb in 1979, with the model subsequently being improved by Wang (1987) and Campbell (1987), who suggested an analytic expression for the toroidal component of the field. This expression has been widely used since then, on account of it being both simple and physically plausible.
In these analytic models, the authors made the kinematic approximation, looking for an axisymmetric stationary solution of the induction equation with a given unchanging structure for the fluid in the disc. They further took the disc to be thin, the poloidal component of the magnetic field to be exactly dipolar, and the velocity field to have zero poloidal component, with its azimuthal component being Keplerian inside the disc^{1} and corotating at the top of the disc. Using cylindrical coordinates (ϖ, φ, z), they found (1)where γ_{a} is the amplification factor, Ω_{K} and Ω_{s} are the Keplerian and stellar angular velocities, respectively, and τ_{d} is the dissipation time scale. The amplification factor γ_{a} was taken to be a constant not much greater than 1 (it depends on the steepness of the transition – in the zdirection – between Keplerian motion inside the disc and corotation with the star at the top of the disc). The precise profile of τ_{d} depends on what is the dominant mechanism for dissipating the magnetic field. Wang (1995) considered three different cases, with τ_{d} being dominated by the Alfven velocity, turbulent diffusion and magnetic reconnection, respectively.
Equation (1) was then used for calculating the net magnetic torque exerted on the neutron star. The vertically averaged torque can be written as (2)where h is the semithickness of the disc. Regions of the disc inward of the corotation point therefore give positive contributions to the torque (because B_{φ} > 0), while the remainder of the disc gives negative contributions (because B_{φ} < 0). The total magnetic torque is obtained by integrating the local values from the inner edge to the outer boundary, and it can be either positive or negative depending on the location of the inner edge of the disc with respect to the corotation point.
The aim of the present paper is to develop a semianalytic model that can improve on those of Wang and Campbell, while remaining simple enough to be useful for people discussing the behaviour of astrophysical sources, giving a conceptual picture to go alongside results from large numerical calculations where the full set of the MHD equations is solved.
Our approach, for the time being, is to continue to retain axisymmetry and the kinematic approximation but to calculate a consistent steadystate solution for the magnetic field, relaxing the assumptions on the poloidal components of the magnetic and velocity fields and using a 2D model without any vertical averaging of the Taylor expansion of the induction equation. In the main region of the outer disc (see Fig. 1) we use a simple Keplerian velocity profile, but this is something that will be improved on later. In a previous paper (Naso & Miller 2010, hereafter Paper I) we analysed the distortion of the poloidal component of the magnetic field using a similar approach, and found that deviations away from a dipole field can be quite significant. Here we focus on the toroidal component and use the results of the previous model to solve the φ component of the induction equation. We find that in general B_{φ} follows a profile different from that of the analytic models, i.e. Eq. (1), and reduces to that only in a very particular case.
Following this introduction, in Sect. 2 we briefly describe our model, which is the same as that of Paper I; in Sect. 3 we recall the equations used (obtained from the induction equation), give expressions for the velocity and diffusivity profiles and outline our solution method (details of tests made on the code are given in an Appendix); in Sect. 4 we present our numerical results; in Sect. 5 we comment on these, comparing them with those of the earlier analytic models, and develop our new suggestion for the B_{φ} profile. Section 6 contains conclusions.
2. Model
In this study, we are considering disc accretion by a neutron star having a dipolar magnetic field. The model is the same as that of Paper I. For a detailed description of it, see Sect. 2 of that paper; here we recall the main points.
Fig. 1
Schematic representation of our model (not to scale). We use r_{in} = 10 r_{g}, r_{tr} ~ 22 r_{g} and r_{lc} ~ 115 r_{g}. The opening angles are 8° for the disc alone and 10° for the disc plus corona. The outer disc extends much further out than the main region shown here: the grid continues until r_{out} = 380 r_{g}. 

Open with DEXTER 
We are assuming that the star is rotating about its magnetic axis, and that this axis is perpendicular to the plane of the disc; also, we assume that the fluid flow is steady and has axial symmetry everywhere. We use the kinematic approximation and do not consider any dynamo action, but turbulent diffusivity is included. The velocity field is not constrained to be purely azimuthal but is allowed to have nonzero components also in the other directions. We use spherical coordinates (r, θ, φ), with the origin being at the centre of the neutron star. Boundary conditions are imposed at the inner and outer radial edges of the disc (r_{in} is at the Alfven radius, and r_{out} is at 38 r_{in}), on top of the corona (taken as being a layer above and below the disc) and on the equatorial plane. Having the inner edge of the disc at the Alfven radius justifies the kinematic approximation to some extent, since in this configuration the magnetic pressure is smaller than the gas pressure within the region that we are considering, and so the effects of the plasma on the magnetic field should be larger than the magnetic feedback on the plasma flow. The ratio h/r is taken to be constant, with the opening angle being 8° for the disc (measuring from the equatorial plane to the top of the disc), and 10° for the disc plus corona.
The inner disc region (r < r_{tr} ~ 2 r_{in}) and the corona are modelled with a larger value of η than the other parts. In these regions, the kinematic approximation does not provide a good description of the system for two different reasons: in the corona this is because of the low density of the plasma (which therefore tends to follow the magnetic field behaviour rather than being followed by it); in the inner region, it is because the magnetic field intensity is still quite large – although the magnetic pressure is smaller than the gas pressure, it is not yet negligible. Using a larger value for η in these regions makes the magnetic field less sensitive to the plasma motion; a somewhat similar approach was used by Kueker et al. (2000). We recall that the present knowledge of the turbulent magnetic diffusivity is quite poor and it is not a simple task to obtain a reliable expression for the η profile.
As regards the velocity field: for v_{r} we use the expression given for the “middle region” of αdiscs by Shakura & Sunyaev (1973). For Ω we take Keplerian rotation in the disc and corotation at the top of the corona and at the inner edge of the disc, giving a maximum for Ω between r_{in} and r_{tr}. These different parts are smoothly connected using error functions. Regarding v_{θ}: we put it to zero in the disc but near to the boundaries we are forced to have a nonzero value in order to be consistent with the dipolar boundary conditions (as shown in Sect. 3.2 of Paper I) and so we use a nonzero profile in the corona. In this way we are including in the model an outflow from the surface of the corona, and this is in agreement with recent hydrodynamic simulations of accretion flows (Jiao & Wu 2011).
Summarising, we divide the surroundings of the central object into four parts (see Fig. 1, which is repeated from Paper I): (1) the inner disc, extending from r_{in} out to a transition radius r_{tr} ~ 2 r_{in} (where the diffusivity changes); (2) the outer disc, extending from r_{tr} to an outer radius r_{out} = 38 r_{in}; (3) a corona, above and below the disc; and (4) everything else, which we take here to be vacuum. As a unit for radial distances, we use the Schwarzschild radius r_{g}. Within the outer disc, we focus on what we call the main region, extending from r_{tr} out to the light cylinder at r_{lc} ~ 11 r_{in}.
3. Equations
In the kinematic approximation, one assumes that the velocity field remains fixed as specified, and the interaction between the magnetic field and the plasma is then described by the induction equation alone. In the presence of turbulence, it is more convenient to write this equation for mean fields rather than for the actual fields (which contain fluctuating parts as well).
The time dependence of the mean field is given by (3)where η_{Ohm} = c^{2}/4πσ is the Ohmic diffusivity and ℰ is the turbulent electromotive force. A common procedure is to expand ℰ in terms of the mean field and its derivatives and within the first order smoothing approximation one has ℰ = α_{T}B − η_{T}∇ × B, where the α_{T}B term generates the socalled αeffect. As in Paper I, we are neglecting this effect here and the induction equation then reduces to (4)where η = η_{Ohm} + η_{T} and is ~ η_{T}, because the turbulent diffusivity is much stronger than the Ohmic one.
We note that the effects of a dynamo action on the disc structure have recently been studied by Tessema & Torkelsson (2010a,b), who estimated the toroidal magnetic field generated by the dynamo to be about an order of magnitude larger than the B_{φ} calculated according to the early models. Here we show that the profile of the toroidal field can be very different from the one suggested by those models, if the poloidal component is not forced to be a dipole, but is instead calculated selfconsistently.
As described in Paper I, our strategy consists of writing Eq. (4) in spherical coordinates, applying the axisymmetry and stationarity assumptions (i.e. putting ∂_{φ} [... ] = ∂_{t} [... ] = 0) and then solving the final equations with the velocity field and magnetic turbulence profile given by the model. The three components of the induction equation are The first two equations contain only poloidal quantities and have been solved in Paper I (making use of the magnetic stream function). Here we focus on the third equation and solve it using the results for B_{r} and B_{θ} from the previous analysis.
We rewrite Eq. (7) in the following dimensionless way: (8)where x = r/r_{g} (not to be confused with the Cartesian coordinate x used for the plots), B_{φ} is the toroidal field in units of a reference field (as described in Sect. 3.4 below) and the dimensionless a coefficients are All of these coefficients have direct analytic expressions, except for the last one (a_{0}) which contains B_{r} and B_{θ}, whose values are taken from the previous numerical calculations. Note that the magnetic field components appearing in Eq. (13) are dimensionless.
3.1. Poloidal velocity and diffusivity
The poloidal components of the velocity field and the diffusivity have already been discussed in Paper I. Here we use the same profiles as before; for the sake of clarity, we give the expressions again below.
For v_{r}, we use (14)where m is the stellar mass in solar mass units, ṁ is the accretion flux in units of the critical Eddington rate and α is the standard ShakuraSunyaev viscosity coefficient. Using typical values (α = 0.1, ṁ = 0.03 and m = 1.4) this gives (15)For v_{θ}, we use (16)where the transition in v_{θ} between the two regions is made using(17)with (18)where θ_{D} is at the upper surface of the disc (i.e. 82°) and d = 10^{3} radians (i.e. 0.057°).
For the diffusivity, we use (19)where η_{0} and η_{c} are the values in the main disc region and the corona (see Fig. 1), for which we use the values 10^{10} cm^{2} s^{1} and 10^{12} cm^{2} s^{1} respectively (we also use other values to study the impact on the results of varying η). For η_{θ}(θ) and η_{r}(r) we use joining functions of the form (20)with x_{c} = θ_{D} and r_{tr} for η_{θ} and η_{r}, respectively, and with d being the width of the transition in the error function, for which we use d = 5 r_{g} in the radial direction and d = 2° in the θ direction^{2}.
3.2. Azimuthal velocity
We are modelling the main disc region as rotating with Keplerian velocity, with the corona being taken as a transition layer where the velocity goes from Keplerian to corotation. In the inner region, we join the Keplerian flow to corotation at the radial inner edge, again using an error function. We recall that in Paper I we showed that a strictly dipolar field, without distortions, can in principle be a stationary solution of the induction equation provided that the velocity field fulfils the two conditions: (21)and (22)From Eq. (22) one sees that corotation, which is obtained by choosing γ = 0, is consistent with having dipolar conditions (which is what we are using here).
For the magnetic field intensities and neutron star spin rates which we are using as standards (B ~ 10^{8} G and P ~ 10 ms), the corotation point is outward of the inner edge of the disc (which is the standard condition required for being in the accretion regime) and therefore Ω should reach a maximum at some location and then decrease again, as one moves inwards. Summarising, we use the following profile: (23)where θ_{C} is the upper surface of the corona, Δθ is the angular resolution, Ω_{s} is the stellar spin rate, Ω_{K} is the Keplerian angular velocity and the two smooth connections are made using the error functions given in Eqs. (25) and (27) below.
As regards the smooth joins, in the θ direction we write (24)where: (25)with d = 10^{2} radians (i.e. 0.57°), and we then have a modification in the radial direction giving (26)where (27)with d = 2 r_{g}. See Fig. 2 for a contour plot of Ω(r,θ) (made in terms of Cartesian coordinates x and z).
Fig. 2
Contour plot of Ω(r,θ). This is the same for all of the configurations. The straight white line indicates the boundary between the disc and the corona. The corotation point is at r = 18.8 r_{g}. 

Open with DEXTER 
3.3. Boundary conditions for B_{φ}
In our model we take the region outside of the disc and its associated corona to be vacuum, and suppose that there is no toroidal component of the magnetic field there (i.e. B_{φ} = 0). We impose the same condition also at the equatorial plane because B_{φ} has to be antisymmetric across it.
When B_{φ} = 0, Eq. (7) reduces to a_{0} = 0, i.e. (28)at the top of the corona, at the inner edge of the disc and on the equatorial plane. We need our choice of the boundary conditions (i.e. purely dipolar field and corotation) to represent a solution of this equation. In Paper I, we showed that corotation is consistent with a pure dipolar field, as already mentioned in the previous subsection (see Eq. (22)). In addition we note here that, on the equatorial plane, Eq. (28) is trivially satisfied because B_{r} = 0 and both v_{φ} and B_{θ} are symmetric with respect to this plane.
3.4. Solution method
In order to get to Eq. (8) we first expanded out all of the derivatives present in Eq. (7) and then put the result into a dimensionless form using x = r/r_{g} as the radial coordinate and measuring B_{φ} in units of , which is a reference value for the magnetic field, taken to be , where B_{0} is the magnitude of the dipolar field at the stellar equator and r_{0} is the dimensionless neutronstar radius. As in Paper I, we use canonical values for the mass and radius of the neutron star, 1.4 M_{⊙} and 10 km respectively, and take B_{0} = 3 × 10^{8} G, as typical for a millisecond pulsar.
Equation (8) is a nonhomogeneous elliptic partial differential equation for B_{φ} and we have solved it using the GaussSeidel relaxation method which uses a finitedifference technique, approximating the operators by discretizing the functions over a grid. At any given iteration step, the values of B_{φ} at the various grid points are written in terms of values at the previous step, or at the present step in the case of locations where it has already been updated. Details of the numerical scheme are given in the Appendix of Paper I, where we used the same method to solve the elliptic partial differential equation for the magnetic stream function.
Before applying the numerical scheme to the actual problem that we want to solve, we performed a series of tests on the code, which are described in detail in the Appendix. We used several configurations, with many different numbers of grid points, profiles of the turbulent diffusivity, initial estimates for the toroidal field, locations for the outer radial boundary of the grid and values for the iteration time step. In this way we have checked the code stability and convergence, have optimised the iteration procedure and have determined where best to place the outer radial boundary of the grid (so that the outer boundary conditions do not significantly influence the solution in our region of interest).
All of the results presented in this paper have been obtained using a grid spacing of Δr = 0.74 r_{g} and Δθ = 0.125°. The final maximum residual was of the order of 10^{13} and saturation was reached after about 4 × 10^{6} iterations, the iteration time step being 4.46 × 10^{4} (about 95% of the critical one, beyond which the method gives a divergent solution). As a first estimate for the profile of the magnetic field we have used a Gaussian. The numerical domain was x ∈ [10,380] , θ ∈ [80,90] (in degrees), which we covered with a homogeneous grid of 501 × 81 points. The profile of the poloidal field (which is present in the expression for the coefficient a_{0}) was calculated by running the code used in Paper I with the same profile of η and v and the same resolutions as used in this analysis. However, since for the poloidal calculation there is a stronger dependence on the outer radial boundary conditions, we have placed r_{out} at 750 r_{g} and used a 1001 × 81 grid. In general the number of radial points used in the poloidal calculation () and in the toroidal one () are related through the following condition (which comes from equating the spatial resolutions): (29)For the poloidal calculation the final maximum residual was of the order of 10^{9} and was reached after about 7 × 10^{7} iterations. The iteration time step was the same as for the toroidal calculation.
4. Results
4.1. Reference configuration
As mentioned in Sect. 3, we have used a slightly different configuration from that considered in Paper I (the transition in the η profile is wider and the resolution in the angular direction is increased). The poloidal field for the new configuration, resulting from solving Eqs. (5) and (6), is shown in Fig. 3. If we additionally assume a profile for the angular velocity, we can then solve Eq. (8). The result, for the profile given by Eq. (23), is shown in Fig. 4 as a contour plot of B_{φ}, for all of the region of interest (i.e. for r ∈ [10,115] r_{g} and including the corona). Contours for positive values of B_{φ}, representing toroidal field lines rotating in the same direction as the neutron star, are shown with black solid lines; while those for negative values of B_{φ} (rotating in the opposite direction) are shown with black dashed lines. Triple dotteddashed white lines show where B_{φ} = 0.
Fig. 3
Poloidal field lines from the numerical solution (solid) and those for a dipole (dotted). In this configuration, v_{0} = 10^{5} cm s^{1} while the diffusivity η_{0} = 10^{10} cm^{2} s^{1} with a transition width in θ direction of 2° (3.5 × 10^{2} radians). Straight lines indicate the top surface of the corona (solid) and the boundary between disc and corona (dashed). 

Open with DEXTER 
Fig. 4
Contour plot of the toroidal field in Cartesian coordinates, in the region r = [10−115] r_{g} and θ = 80°−90°. The black solid lines are lines of positive B_{φ}, while black dashed lines are used for negative B_{φ}. B_{φ} = 0 is shown with triple dotteddashed white lines. The straight white line indicates the boundary between disc and corona. 

Open with DEXTER 
Fig. 5
The same as in Fig. 4 but showing an expanded view of the region r = [10−35] r_{g}. The corotation point is at r = 18.8 r_{g}. 

Open with DEXTER 
The toroidal field shows a quite structured profile, with two (positive) maxima and two (negative) minima. Surprisingly, the global maximum is located outward of the corotation point and is positive, in contrast to the standard picture according to which the sign of B_{φ} is the same as that of the relative angular velocity between the disc matter and the star. There is also a quite striking vertical structure, the global minimum being located just above and to the left of the global maximum. As in the analytic models, B_{φ} becomes zero near to the corotation point (r_{cor} = 18.8 r_{g}), but here one sees that for almost every value of r > r_{cor} there is a value of θ where the field is zero, so that there are places with zero B_{φ} throughout all of the main disc region.
In Fig. 5 we show the B_{φ} contour plot in the region r ∈ [10,35] r_{g}, so that the structure of the magnetic field in this part can be seen more clearly. In Table 1 we give the coordinates of the minima and maxima and the magnitude of the toroidal field at those locations, measured in units of the stellar field strength B_{0} (which is taken to be 3 × 10^{8} G). All of these main features are located in the region r ∈ [10,55] r_{g}; the remainder of the disc can be divided into two zones: one where B_{φ} > 0 (in Figs. 4 and 5 this is below the long white triple dotteddashed curve that crosses the equatorial plane at x = 23 and x = 110), and the other where, instead, it is negative. The minimum latitude of the region with positive B_{φ} is about 84°.
Fig. 6
Toroidal field strength plotted against r at different values of θ, in the region r = [10−50] r_{g} (it is equal to zero on the equatorial plane). Negative values mean that B_{φ} is pointing backwards with respect to the disc rotation. The thick solid line is the shell average of B_{φ}. 

Open with DEXTER 
Locations of extrema of the toroidal field.
Radial profiles of B_{φ} are shown in Fig. 6 for several values of the latitude and for the shell average. We show the profiles for θ = 87.1° (where the global maximum is), θ = 83.6° (where the local minimum is, and which is very close also to the local maximum and the global minimum) and for an intermediate value θ = 85°. The curves all pass through zero at the corotation point (r = 18.8 r_{g}), at least to within the accuracy of the calculation. The large positive B_{φ} peak at about 30 r_{g} is progressively reduced as one moves from the mid plane of the disc towards the corona, and eventually becomes a negative local maximum. One can calculate the shell average of the radial profile, i.e. an average of B_{φ} over θ for each r, and this is also shown in Fig. 6 as a thick solid line. This average reproduces quite well the general behaviour of the toroidal field and shows all of its main features.
The θ dependence of the sign of B_{φ} is a key result, that can have quite dramatic consequences for the calculation of the magnetic torque exerted on the neutron star. In fact, it is usually assumed that the torque depends only on r, being positive inside the corotation radius and negative outside it (see Eqs. (1) and (2)), whereas we now see regions of positive B_{φ} even outward of the corotation point. Therefore we need to rethink the discussion of which regions in the disc tend to spin the star up or down. An appropriate approach for calculating the torque requires integration in both directions: we plan to perform this analysis in a future work.
Finally, we note here that the magnitude of the toroidal component of the magnetic field is typically larger than that of the poloidal component. For example, the maximum of the shell average of the poloidal component is ~ 3.6 × 10^{3}, while the maximum of ⟨ B_{φ} ⟩ _{θ} is ~ 2.3 × 10^{1}. In terms of energy conservation, one has to bear in mind that any change in the magnetic energy has to be compensated by a corresponding opposite change in the plasma energy, while in the present model we are taking the flow pattern to be fixed. It is important for the backreaction to be consistently taken into account and this will be done in subsequent work. Also, the distortions of the toroidal field will be limited by magnetic reconnection.
4.2. Configurations with larger η_{0}
We have already mentioned in Sect. 2 that there are big uncertainties about how to model the turbulent magnetic diffusivity within the disc and the surrounding corona. This quantity is often discussed in terms of the turbulent magnetic Prandtl number, P_{m} ≡ ν/η, which links it with the turbulent viscosity.
Within the kinematic approximation, one does not solve selfconsistently for the velocity field and so it is necessary assume some profile for it. As outlined in Sect. 3.1, in the present calculations we are using the velocity prescription given by the Shakura & Sunyaev thin disc model, which also embodies a particular connection between the viscosity and other disc quantities. This gives (following e.g. Szuszkiewicz & Miller 2001) ν = α h c_{s}, with (assuming vertical hydrostatic equilibrium), where c_{s} is the sound speed, p is the pressure and ρ is the density. Using the isothermal sound speed, the Keplerian angular velocity profile and the parameters given in Sect. 3.1, one gets ν ~ 10^{13} cm^{2} s^{1} which, in turn, gives P_{m} ~ 10^{3} in the disc and ~10 in the corona of our reference model. These values may seem rather high, but one should be cautious about using them because there are several further factors which need to be taken into account.
Note, first, that calculating the viscosity with the α disc model certainly gives an overestimate for P_{m}, because one is neglecting the contribution of the magnetic field in the equation of motion. Also, there is a degeneracy in the profile of the velocity field, in the sense that one can obtain the same velocity profile, and hence the same results for our numerical calculations, using different combinations of α and ṁ: an increase/decrease by a factor of λ in the accretion rate, together with a decrease/increase by a factor of in α (and hence ν), causes no change in the velocity profile. These considerations can easily bring down the true values of P_{m} for our calculations well below the approximations quoted above. In any case, there is currently no general consensus about the correct value for P_{m}: we note that rather high values have recently been found in some numerical MRI simulations (e.g. see Takahashi & Masada 2010; Romanova et al. 2011).
We chose our values for η bearing in mind which values would give significant field distortions in the disc and therefore be most interesting. However, it is clearly important to investigate the effects of varying these numbers, and so we also made some calculations using larger values of η (smaller P_{m}). We show here results for η_{0} = 4 × 10^{10} cm^{2} s^{1} and η_{0} = 10^{11} cm^{2} s^{1} with η_{c} being, as usual, two orders of magnitude larger (in the reference configuration we used η_{0} = 10^{10} cm^{2} s^{1}).
Results for both configurations are presented in Figs. 7–9, where we show the poloidal field lines, the contour plot of the toroidal field component and the radial profile of its shell average, respectively^{3}.
Fig. 7
Poloidal magnetic field lines. Comparison between two configurations having values of η_{0} larger than that of the reference configuration. 

Open with DEXTER 
Fig. 8
Contour plots of B_{φ}. Comparison between two configurations having values of η_{0} larger than that of the reference configuration. 

Open with DEXTER 
Increasing η_{0} (i.e. decreasing the magnetic distortion function D_{m}, our generalisation of the magnetic Reynolds number introduced in Paper I) makes the field diffuse more efficiently and therefore the solution gets progressively nearer to being a dipole field (see Fig. 7). Reducing the poloidal field distortions in turn changes the toroidal component, and the modifications are quite substantial. The structure with four extrema gradually turns into one with only two (see Fig. 8), where the sign of B_{φ} is positive for radii smaller than the corotation radius and negative for larger ones. This transition is clearly seen when plotting the shell average of the toroidal field (see Fig. 9). The behaviour of the sign of the shell average of B_{φ} is then the same as that for B_{φ} in the early analytic models, although the details of the profiles have significant differences (see Fig. 13 below).
Fig. 9
Shell average of the toroidal component. Comparison between two configurations having turbulent magnetic diffusivity, η_{0}, larger than that of the reference configuration. The values of η_{0} are 4 × 10^{10} cm^{2} s^{1} for the solid curve and 10^{11} cm^{2} s^{1} for the dashed curve. 

Open with DEXTER 
Finally we note that the magnitude of the toroidal component decreases with increasing η_{0}, the maximum of the shell average being ~14 × 10^{3} B_{0} for the configuration with η_{0} = 4 × 10^{10} cm^{2} s^{1} and ~6 × 10^{3} B_{0} for the one with η_{0} = 10^{11} cm^{2} s^{1} (compare the shell averages in Fig. 9 with each other and with the one in Fig. 6), while the maximum for the poloidal component has remained at ~3.5 × 10^{3} B_{0} (this is because the maximum for the poloidal componant occurs at the inner edge, where the field depends more on the boundary conditions than on the value of η).
4.3. Configurations with constant η
Although the η profile that we have used so far is the one that we consider to be the most appropriate when studying accretion within the kinematic approximation (for the reasons given in Sect. 2), we have considered also configurations where the diffusivity is constant through all of the disc and the corona, in order to have a clear understanding of how sensitive the results are to this quantity. Cases with η = 10^{10} cm^{2} s^{1} and 4 × 10^{10} cm^{2} s^{1} proved unsatisfactory because of having very abrupt changes away from the dipole field at the top of the corona (this is exactly the behaviour that we have tried to avoid by using a larger value of η in the corona in our reference configuration). Results for η = 10^{11} cm^{2} s^{1} and η = 10^{12} cm^{2} s^{1} are shown in Figs. 10–12.
Fig. 10
Poloidal magnetic field lines. Comparison between two configurations with constant η. 

Open with DEXTER 
Fig. 11
Shell average of the toroidal component for η = 10^{11} cm^{2} s^{1}. For η = 10^{12} cm^{2} s^{1}, the curve has exactly the same shape, but the values are about 100 times smaller (in absolute value). 

Open with DEXTER 
Fig. 12
Contour plots of B_{φ}. Comparison between two configurations with constant η. 

Open with DEXTER 
For configurations with constant η, the analysis is made simpler since there are no regions with diffusivity gradients. As stated in the previous subsection, using a larger value of η reduces the deviations away from the dipolar field. However, in contrast with the previous case, the distortions are now more uniform throughout the disc (compare Figs. 7 and 10), because η is constant and the magnetic distortion function is monotonically decreasing with r (while previously ∂_{r}D_{m} had a peak). We note that deviations away from a pure dipole are very small when using η = 10^{12} cm^{2} s^{1}, and for η comparable to or larger than this, the poloidal component of the field can be well approximated by a dipole. As regards the toroidal component, it has only two extrema which are both located just beneath the surface of the disc (see Figs. 11 and 12). These two extrema have been observed in all of the configurations which we have studied; they are also present in the models of Wang (1987) and Campbell (1987, 1992) and therefore seem to be robust features.
5. Discussion
5.1. Comparison with analytic models
In the analytic models of Wang (1987) and Campbell (1987), who we will refer to now as W&C, the toroidal component of the magnetic field is written as being proportional to the relative angular velocity between the disc and the star multiplied by the vertical field, which is taken to be a pure dipole (see Eq. (1)). For models where Ω_{disc} is Keplerian and the inner edge of the disc is inward of the corotation point, B_{φ} has a positive global maximum at the inner edge of the disc, is zero at the corotation point, reaches a global negative minimum somewhere outward of this and then tends towards zero at very large r. If Ω deviates from Keplerian in the inner part of the disc (reaching a maximum and then decreasing again as one moves inwards), then the location of the maximum of B_{φ} is not in general at the inner edge, but depends on the precise profile of Ω in this inner region.
The above description is only partially in agreement with the results of our present two dimensional calculations. They share the feature of having a zero of B_{φ} at the corotation point (or extremely close to it, see Fig. 6). As regards the predicted global maximum of B_{φ} in the inner part of the disc, all of our calculations show a positive maximum close to where Ω has a maximum (compare Figs. 5, 8 and 12 with Fig. 2). However, when η is not constant, magnetic field lines accumulate and a second maximum can appear outward of the corotation point, and this can also become a (positive) global maximum depending on the value of η_{0}. Finally, regarding the minimum: as in the W&C models, we always find a global negative minimum before the outer edge of the disc; however, when η is not constant a second minimum can appear, thus producing a structure with four extrema: two maxima and two minima (see Fig. 4). We note here that even in the configurations with only two extrema, as in W&C, the profile of the toroidal component is still different from that predicted by those models. In particular, the location of the minimum and the magnitude of the field at both extrema can be very different.
The main features of the toroidal field component, as given by our present numerical calculations and by the analytic models, can be seen in Fig. 13, where we show three profiles for B_{φ}: the dotted curve is the W&C profile for our model, as resulting from Eq. (1) but where we have used our profile for the angular velocity near to the equatorial plane (i.e. Eq. (23) at θ = 88°), the solid curve shows the shell average for our reference configuration, and the dashed curve is for the configuration with constant η (10^{11} cm^{2} s^{1}).
Fig. 13
Comparison of three different profiles of the toroidal field: (1) the solid curve is the shell average for the reference configuration, (2) the dashed curve is the equivalent one for the configuration with constant η (10^{11} cm^{2} s^{1}), (3) the dotted curve is given by the W&C formula with our Ω profile at θ = 88°. All of the curves have been normalised so as to have their peak at 1. 

Open with DEXTER 
The properties of having additional extrema of B_{φ} and of varying the location and magnitude of the two standard extrema, are not seen in the W&C models.
5.2. Role of B_{p} and D_{m}
In Fig. 14, we show the quantity ΔΩ B_{θ}, where ΔΩ = Ω_{disc} − Ω_{star} and B_{θ} is the θcomponent of the poloidal field as obtained from our numerical calculations for the reference configuration. This quantity has one global maximum and one global minimum. The maximum is located very close to where B_{φ} and Ω have their first maxima, and the minimum is at a radial location close to that of the first minimum of B_{φ} (only a few r_{g} smaller – see Fig. 5 and Table 1). This shows that the quantity ΔΩ B_{θ} is still relevant for grasping the fundamental properties of the toroidal component of the field, although care must be taken in choosing the profile of B_{θ}. Some differences between the predictions of the W&C models and our numerical results can, in fact, be explained in terms of the distortion of the θ component of the magnetic field. However, in Fig. 14 there is no evidence for the additional maximum and minimum, and so this is not the whole story. We need to consider the distortions of the field in more detail, and not focus only on the quantity ΔΩ B_{θ}.
Fig. 14
Contour plot of ΔΩ B_{θ} for the reference configuration. 

Open with DEXTER 
In Paper I we have shown that the distortions of the poloidal component due to the plasma motion can be well described by a generalisation of the magnetic Reynolds number (which we have called the magnetic distortion function), defined as^{4}(30)We studied D_{m} in detail in Paper I (see Sect. 4 of that paper); here we just recall that the magnitude of the peak in its radial derivative is a measure of the degree of accumulation of poloidal field lines in its vicinity. For the cases with constant η, the field amplification caused by this distortion is negligible (the distortion is more homogeneous and the field lines do not accumulate) and D_{m} is just proportional to v_{r} in all of the disc (so that there is no peak at all in the radial derivative). For the reference case, instead, ∂_{r}D_{m} does have a peak and its radial location (r ~ 28.5 r_{g}) is near to that of the additional maximum (r ~ 30 r_{g}). Increasing η_{0} reduces the magnitude of the peak and also the additional maximum of B_{φ} becomes less pronounced, eventually disappearing for ∂_{r}D_{m} ≤ 0.037 (see Fig. 15). We can therefore draw the conclusion that it is the radial derivative of the magnetic distortion function which is responsible for the additional maximum in the toroidal field profile.
Fig. 15
Radial derivative of the magnetic distortion function D_{m} near the equatorial plane (at θ = 89°) for three configurations: (1) the solid line is for the reference configuration with η_{0} = 10^{10} cm^{2} s^{1}, (2) the dotted line is for the configuration with η_{0} = 4 × 10^{10} cm^{2} s^{1}, (3) the dashed line is for η_{0} = 10^{11} cm^{2} s^{1}. 

Open with DEXTER 
5.3. Our picture for B_{φ}
In this subsection, we develop our alternative picture for the structure of the toroidal component of the magnetic field, following the same general approach as W&C, but with a more detailed representation of the poloidal magnetic field and the velocity field, and a twodimensional model. This is still a very simplified picture but we believe that it can represent a useful step forward, giving some additional insights. Our starting point is the φ component of the induction equation: (31)The advective term is the one responsible for generating the field, while the diffusive term causes field losses. We then define the following scalar quantities: In a steady state ∂_{t}B_{φ}_{ + } = ∂_{t}B_{φ}_{ − } so that G = L and ∂_{t}B_{φ} = 0.
We focus first on the generation term and split it into two parts: one involving the poloidal motion (G_{p}) and the other involving the toroidal motion (G_{φ}): (34)In spherical coordinates, these two terms are written as From these expressions we can see that, in general, the generation rate depends on several quantities and not only on the vertical gradient of the angular velocity (as in the W&C models). In fact, all of the components of the velocity field and magnetic field are present. However, usually in discs v_{φ} ≫ v_{r},v_{θ}^{5} and one can expect the second term to dominate. An important difference between G_{p} and G_{φ} is that the former cannot generate any toroidal field from zero, but can only modify B_{φ} once it has already been produced by some other means, e.g. coming from advection by the azimuthal motion (i.e. from G_{φ}).
We have calculated the ratio G_{φ}/G_{p} for our configurations and have found that the toroidal term always dominates, even if the ratio is not as large as would have been expected just by comparing the velocities (because one should consider also the magnetic field). However even neglecting the contribution from the poloidal advection, the generation rate for B_{φ} is still considerably different from the one considered in the W&C models. In fact, it contains also the radial derivative of the term r v_{r} B_{r} and the vertical derivative of B_{θ}.
So far, we have focused only on the generation term. In order to obtain the profile of the toroidal field we have to equate it to the loss term, for which we have If we consider only typical values ( and ), then we can write (39)In a steady state, when L = G ≃ G_{φ}, one has , and one gets (40)In our configuration, η has two main characteristic values: η_{0} in the main disc region and η_{c} = 10^{2} η_{0} in the corona and in the inner part of the disc. If in Eq. (40) we replace with the actual profile of the magnetic diffusivity, and choose , then we find that the variation of B_{φ} can be well approximated by that of the coefficient a_{0} in Eq. (8): (41)where B_{0} is the reference unit of the magnetic field and we recall that in Eq. (8) the fields are dimensionless (while here they are dimensional).
This formula for B_{φ} is rather well confirmed by our numerical calculations. In Fig. 16 we show the contour plot for a_{0}: this is very similar to that for the toroidal field (compare with Fig. 4), the differences being due to the approximations made in evaluating the loss term L, and to neglecting G_{p} in the generation term.
Fig. 16
Contour plot of a_{0} in Cartesian coordinates, in the region r = [10−115] r_{g}. This should be compared with Fig. 4. 

Open with DEXTER 
5.4. Consistency check with the W&C models
From Eq. (41) we can see that the profile of the toroidal field is basically determined by four factors: (1) the radial derivative of B_{r}, (2) the radial derivative of r v_{φ}, (3) the vertical derivative of B_{θ} and (4) the vertical derivative of v_{φ}.
In the W&C models, the first three of these are neglected because B_{r} is taken to be zero everywhere and B_{θ} is supposed not to vary with θ. If one adds also the other assumptions used in their models (about the velocity profile and the disc thickness), one then finds that Eq. (41) reduces to (42)which is the same as Eq. (1) with γ_{a} = (ϖ/h) and . We can therefore recover the expression appearing in the earlier analytic models from our result.
We now want to check a posteriori whether or not the simplifications made in the W&C models are still valid in our more general 2D model and, if so, in which parts of the disc. One can do this by calculating the ratio (Q) between the two terms on the righthandside of Eq. (41): (43)Clearly, a necessary condition required for matching with the W&C models is that Q ≫ 1. In Fig. 17 we show a contour plot of this quantity for the reference configuration. We show in red all of the regions where Q > 1, while regions having the condition clearly being violated are colourcoded to show by how much Q is smaller than 1.
Fig. 17
Contour plot of Q for the reference configuration. All locations with Q > 1 are marked in red so as to highlight the regions where Q < 1. 

Open with DEXTER 
There are large parts of the disc where this necessary condition is not met, with Q even reaching values as small as 10^{3}. Moreover, even in regions where Q ≫ 1, the further necessary condition ∂_{θ}B_{θ} ≪ ∂_{θ}Ω may not be met. In the main disc region the angular velocity is almost constant with θ (the transition from Keplerian rotation to corotation happens in the corona), and here we are exactly in an opposite regime, i.e. ∂_{θ}B_{θ} ≫ ∂_{θ}Ω ~ 0. Only in the corona and in the upper part of the disc are the vertical gradients of Ω not negligible.
For a pure dipolar field, Eq. (43) becomes (44)which, for Keplerian motion, gives exactly 2/5. Therefore the necessary condition is never met for a pure dipole and Keplerian rotation, which is a good description for the parts of our discs near to the mid plane. This should not surprise us, since two key assumptions made in the W&C models hold only in different parts of our discs and not together: the field was supposed to have both B_{r} and ∂_{θ}B_{θ} vanishing and decaying as for a dipole (which we have only very close to the equatorial plane or for large η) and the transition to corotation in the angular velocity was supposed to be very sharp (which we have just beneath the disc surface, where Ω has to become equal to Ω_{s}, far from the equatorial plane). We have calculated Q for a configuration with constant η (10^{12} cm^{2} s^{1}): in this case the diffusivity is so strong that deviations away from the dipole are quite small and so changes in Q come almost entirely from the angular velocity (according to Eq. (44)). As expected, in the lower part of the disc we find Q ~ 0.4.
Therefore having a larger value of η does not necessarily imply that the necessary conditions hold in a larger region of the disc, it just implies that the field is closer to a dipole (which is what we are imposing at the boundaries). In order to get Q to go to infinity (which is what was assumed in the W&C models) not only does one need B_{r} = 0 and ∂_{θ}B_{θ} = 0 ^{6} but also that the vertical gradient of the angular velocity, at the same location, has to be nonzero, or much larger than the vertical derivative of B_{θ}.
6. Conclusions
In this paper we have considered a system consisting of a rotating neutron star, having a dipole magnetic field aligned with the rotation axis and surrounded by an accretion disc. The disc is truncated at the Alfven radius and has a coronal layer above and below it. The region outside the corona is taken to be vacuum and we impose dipolar boundary conditions on all of the boundaries.
Our aim was to improve on the analytic models developed by Wang (1987) and Campbell (1987) (W&C). As in those models, we have made the kinematic approximation and have looked for an axisymmetric stationary solution of the induction equation, but we have gone beyond those models in solving for all of the components of the magnetic field and not assuming the poloidal component to be dipolar within the disc. We have also retained all of the components of the velocity field rather than putting v_{r} and v_{θ} to zero everywhere. We have performed a fully twodimensional calculation, without making any vertical average or Taylor expansion in h (the semiheight of the disc). Finally we have neglected dynamo action but have included a turbulent magnetic diffusivity.
The analysis of the poloidal component of the magnetic field has been presented in a previous paper (Naso & Miller 2010, Paper I); in the present paper we have focused on the toroidal component. We have solved the φ component of the induction equation numerically and have shown that the profile obtained for B_{φ} can be very different from that in the earlier analytic models.
In the W&C models, the toroidal field strength was taken to be proportional to the relative angular velocity between the disc and the central object multiplied by the vertical field, which was taken to be dipolar. However in Paper I we found that, when calculated consistently, the poloidal field component was often far from being dipolar. Therefore a first improvement with respect to the earlier models was to use the poloidal field as obtained in our calculations, i.e. a field dragged inwards by the plasma motion. This behaviour explains why we then find different intensities for the toroidal field, and also a different location for its global minimum.
When the turbulent magnetic diffusivity η increases or the radial velocity decreases one expects the field to be progressively less distorted by the plasma motion. This is indeed what we have found both here and in Paper I. Our results show that when the diffusivity η_{0} is larger than about 10^{12} cm^{2} s^{1}, with the characteristic velocity v_{0} being of the order of 10^{5} cm s^{1}, then the field is barely modified. Therefore whenever we expect the magnetic field to deviate from the stellar dipole, η should not be larger than 10^{7} v_{cgs} cm^{2} s^{1}, where v_{cgs} is the characteristic magnitude of the radial velocity expressed in cm s^{1}.
When the turbulent magnetic diffusivity is not constant throughout the disc and corona, two additional extremal points may well appear: if the radial derivative of the magnetic distortion function D_{m} is larger than a critical value (about 0.04 in the equatorial plane), there is an additional maximum and minimum, and in some cases B_{φ} can even become positive again outward of the corotation point, so that there are additional locations where B_{φ} = 0. It is clear that under these circumstances the picture of which regions of the disc tend to spin the star up or down has to be radically redrawn (this will be the subject of a future investigation). However we should emphasise here that there are still many uncertainties among experts about which profile of η should be used and we have therefore made very simple choices here in line with our stepbystep approach.
We have presented a new suggestion for the B_{φ} profile, which reduces to that of W&C if one imposes B_{r} = ∂_{θ}B_{θ} = 0 and Ω = Ω_{K}. In general there are large parts of the disc where the additional terms included in our new picture for B_{φ} dominate over the one retained by W&C (see Fig. 17). Our simplified expression (Eq. (41)) reproduces the numerical results quite well (compare Figs. 4 and 16), the differences being due to approximations made in calculating the generation and loss terms for B_{φ}.
Summarising, in our calculations we have found that B_{φ} can have two maxima and two minima (see Fig. 4). The first maximum (positive and inward of the corotation point) and the first minimum (negative and outward of the corotation point) can be explained referring to the quantity ΔΩ B_{θ}, which has two extrema at the same locations as for B_{φ} (see Fig. 14). These extrema appear also in the W&C models, where the toroidal field is, in fact, taken to be proportional to ΔΩ B_{z}. There is a fundamental difference however: in W&C B_{z} is a pure dipole, whereas the B_{θ} which we consider here is that of a field being dragged inward by the motion of the accreting material. When η is not constant, there is an additional maximum whose magnitude and sign depend on the diffusivity in the disc, and whose radial location is always outward of the first minimum, coinciding with that of the maximum in the radial derivative of the magnetic distortion function D_{m}. Outward of this maximum, the field tends to come back to the profile that it would have had if η were constant, and this produces the last minimum (compare Figs. 9 and 11).
The main conclusion of this analysis is that, when the poloidal component of the magnetic field is treated selfconsistently in the calculations, the profile for B_{φ} can be significantly different from that obtained by W&C, and the magnetic torque generated by it would then be different as well. Moreover, when the turbulent magnetic diffusivity η is not constant throughout the disc and corona, some additional unexpected features can appear (such as a region of positive B_{φ} outward of the corotation point). In the present work, we have retained the very simple Keplerian rotation law in the main part of the disc. Even within a purely hydrodynamical treatment, more complicated velocity fields than this are expected (see Kluzniak & Kita 2000; Jiao & Wu 2011) and further changes are expected when backreaction from the magnetic field on the velocity field is included. The effects of this will be another topic for investigation in subsequent stages of our stepbystep approach.
Note that the errorfunction widths in the radial and angular directions are larger than those used in Paper I (there d_{r} = 2 r_{g} and d_{θ} = 10^{3} radians = 0.057°). The reasons for this are explained in Appendix A.2.
Acknowledgments
This work has been partially supported by CompStar, a Research Networking Programme of the European Science Foundation, and by the National Natural Science Foundation of China (40890161, 2011CB811403, 11025315, 10873020 and 10921303). L.N. is currently supported by a Chinese Academy of Sciences fellowship for young international scientists (grant No. 2010Y2JB12) and would also like to thank the Department of Physics (Astrophysics) of the University of Oxford for support granted during the development of this work, CAMK (Centrum Astronomiczne im. M. Kopernika) in Warsaw, which provided partial support through Polish Ministry of Science grant N N203 381436, and SISSA (Trieste) whose high performance computing facilities have been used for running the numerical calculations.
References
 Campbell, C. G. 1987, MNRAS, 229, 405 [NASA ADS] [Google Scholar]
 Campbell, C. G. 1992, GApFD, 63, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, P., & Lamb, F. K. 1979, ApJ, 232, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Jiao, C.L., & Wu, X.B. 2011, ApJ, 733, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Kluzniak, W., & Kita, D. 2000, unpublished [arXiv:astroph/0006266] [Google Scholar]
 Naso, L., & Miller, J. C. 2010, A&A, 521, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rekowski, M. V., Ruediger, G., & Elstner, D. 2000, A&A, 353, 813 [NASA ADS] [Google Scholar]
 Romanova, M. M., Ustyugova, G. V., Koldoba, A. K., & Lovelace, R. V. E. 2011, MNRAS, submitted [arXiv:1102.1089v1] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Szuszkiewicz, E., & Miller, J. C. 2001, MNRAS, 328, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, H. R., & Masada, Y. 2011, ApJ, 727, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Tessema, S. B., & Torkelsson, U. 2010, A&A, 509, 45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tessema, S. B., & Torkelsson, U. 2011, MNRAS, 412, 1650 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M. 1987, A&A, 183, 257 [NASA ADS] [Google Scholar]
 Wang, Y.M. 1995, ApJ, 449, L153 [NASA ADS] [Google Scholar]
Appendix A: Testing of the code
In this Appendix we discuss some of the tests that we have performed on the numerical code used to solve Eq. (8). For a description of the GaussSeidel relaxation procedure and of the discretization scheme see Appendix A.1 of Paper I.
Before describing the tests, we should underline a difference in the boundary conditions with respect to the code used for the poloidal analysis. In Paper I we were not imposing the dipolar boundary conditions on the magnetic field directly but rather on the magnetic stream function . The poloidal magnetic field was then calculated by differentiating . Because of this, B_{r} and B_{θ} were not precisely dipolar on the boundaries. For the calculations in the present paper, we have introduced a row of ghost points, where we set the field to be exactly dipolar, i.e. the poloidal component is a pure dipole and the toroidal component is zero.
We divide the tests into two groups. In the first group we chose the coefficients of Eq. (8) in such a way that it was possible to find an analytic solution, while in the second group we used the values given by Eqs. (9)–(13). Here is a schematic description of these tests:

Tests with analytic solutions. In addition tofixing the coefficients, one also has to choose the boundaryconditions. We considered three subcases:
 1.1
all coefficients constant and set to 1.0. There is then the following analytic solution: (A.1)
 1.2
all coefficients set to 0.0 except for The general solution is then: (A.4)where k_{1} is a function of b and k_{2}. We chose h = 10, k_{2} = 36^{2} and b = 2 − k_{2}, so as to include a complete period of the angular part of the solution within our angular domain. The solution is then: (A.5)
 1.3
the same choice of coefficients as in test 1.2 but with different boundary conditions: h = 10, k_{2} = 0 and b = 70. The analytic solution is then: (A.6)
 1.1

Testing the model setup. In this group of tests we used a setup which was very similar to that used for our actual physical analysis but varied some numerical parameters so as to test the code. We performed four tests aimed at:
 2.1
checking convergence by changing the number ofgridpoints;
studying the dependence of the solution on the location of the radial outer boundary;
studying dependence on the initial estimate for the solution;
optimising the iteration step size by simple benchmarking.
 2.1
A.1. Tests with analytic solutions
For all of the tests in this category, we considered the code stability and convergence. We used grids with different numbers of points in both directions and compared the analytic errors and the solutions.
In all cases, the stability of the code was related to the size of the iteration step, the code being stable for values smaller than a certain threshold.
To check convergence, we considered how the maximum of the analytic errors changed with varying the total number of iterations. We considered both absolute and relative errors and analysed them in the region of interest (i.e. for r < r_{lc}) by calculating their maximum and looking at their 2D profile. As the iteration procedure continues, the errors decrease and at a certain point they saturate, so that making more iterations no longer leads to smaller errors. In addition to the errors, we have also considered the evolution of the root mean square (rms) of the solution.
The size of the errors at saturation depends on the grid resolution, being smaller for grids with more points (the size of the domain is fixed, so that increasing the number of points means increasing the resolution). Moreover the improvement obtained when increasing the number of gridpoints becomes progressively smaller, as it should do in a convergent regime. We calculated an effective order of convergence p_{eff} by considering the maximum value of the relative error at the final iteration, and how it changed with the grid size. Doing this we obtained p_{eff} ~ 1.5.
All of the tests gave satisfactory results, confirming stability and giving convergence within ~10^{4} iterations for tests 1.1 and 1.2, and within ~10^{7} iterations for test 1.3. The maximum saturation error was ~10^{9}% for test 1.3 and ~10^{3}% for test 1.1. As regards test 1.2, since the solution has some zeros in the considered domain, we could not estimate the error by considering the relative error (since it diverges at the zeros). We instead considered the maximum of the absolute error and compared it with the rms of the solution, finding that it was less than 1%.
A.2. Tests with the model setup
For these tests we had no analytic solutions, and so could not calculate analytic errors but only residuals. When considering stability, we looked at the residuals, while for testing convergence we considered the change in the θaverage of the solution (or in its rms) when changing the number of gridpoints.
We recall that we cannot change the number of gridpoints freely. In fact, one of the coefficients (a_{0}) is not calculated from an analytic function but comes instead from the numerical results of Paper I for the poloidal field. This coefficient is therefore directly defined only on the grid used in that calculation, which had 1001 × 21 points. We consider this grid as our reference one. When using a grid with fewer points, we have to choose them as a subset of our reference grid, while when using a larger number of points, we have to interpolate. Similarly when changing the size of the domain (by reducing r_{out}), we must also decrease the number of grid points N_{i} accordingly, as explained at the end of Sect. 3.4 (see also Eq. (29)).
A very important result of test 2.1 is that if the transition in η between the disc and the corona is not wellenough resolved, then we find convergence to a different solution. More precisely, using an angular width for the transition of 5 × 10^{2} radians we need to have at least 20 grid points overall in the θ direction in order to converge to the correct solution, i.e. to the same solution as for grids with larger values of N_{j} (we tested with N_{j} = 39 and N_{j} = 77). This gives a minimum number of angular zones required for convergence of about 5 within the transition region.
Test 2.2 tells us that the solution is not very dependent on the location of the radial outer boundary, contrary to the situation in Paper I for the poloidal field. The θaveraged solution obtained using r_{out} = 290 differs from that obtained with r_{out} = 750 by less than 0.02% in all of the region of interest. If we consider the solution rms, the difference is even smaller, being less than 5 × 10^{4}%.
For test 2.3, we ran calculations with quite different, and even unphysical, initial estimates for B_{φ}, in order to see whether the solution still converged correctly. We used (1) a constant value B_{φ} = 1; two decaying profiles: (2) B_{φ} = (10/x)^{3} and (3) B_{φ} = (10/x); (4) a growing profile and (5) a Gaussian profile in both directions (centred at x = 100 r_{g}, θ = 85° with widths 15 r_{g} and 2°). We obtained the same final solution for all of the cases; profiles (3) and (5) converged after ~2 × 10^{6} iterations, while all of the others converged after ~4 × 10^{6} iterations.
Finally in test 2.4 we changed the iteration step size Δt, looking for the largest possible value still giving stability. As in Paper I, we found that the maximum value depends more sensitively on N_{j} than on N_{i}. For N_{j} = 81 we found Δt_{max} = 4.7476 × 10^{4}.
All Tables
All Figures
Fig. 1
Schematic representation of our model (not to scale). We use r_{in} = 10 r_{g}, r_{tr} ~ 22 r_{g} and r_{lc} ~ 115 r_{g}. The opening angles are 8° for the disc alone and 10° for the disc plus corona. The outer disc extends much further out than the main region shown here: the grid continues until r_{out} = 380 r_{g}. 

Open with DEXTER  
In the text 
Fig. 2
Contour plot of Ω(r,θ). This is the same for all of the configurations. The straight white line indicates the boundary between the disc and the corona. The corotation point is at r = 18.8 r_{g}. 

Open with DEXTER  
In the text 
Fig. 3
Poloidal field lines from the numerical solution (solid) and those for a dipole (dotted). In this configuration, v_{0} = 10^{5} cm s^{1} while the diffusivity η_{0} = 10^{10} cm^{2} s^{1} with a transition width in θ direction of 2° (3.5 × 10^{2} radians). Straight lines indicate the top surface of the corona (solid) and the boundary between disc and corona (dashed). 

Open with DEXTER  
In the text 
Fig. 4
Contour plot of the toroidal field in Cartesian coordinates, in the region r = [10−115] r_{g} and θ = 80°−90°. The black solid lines are lines of positive B_{φ}, while black dashed lines are used for negative B_{φ}. B_{φ} = 0 is shown with triple dotteddashed white lines. The straight white line indicates the boundary between disc and corona. 

Open with DEXTER  
In the text 
Fig. 5
The same as in Fig. 4 but showing an expanded view of the region r = [10−35] r_{g}. The corotation point is at r = 18.8 r_{g}. 

Open with DEXTER  
In the text 
Fig. 6
Toroidal field strength plotted against r at different values of θ, in the region r = [10−50] r_{g} (it is equal to zero on the equatorial plane). Negative values mean that B_{φ} is pointing backwards with respect to the disc rotation. The thick solid line is the shell average of B_{φ}. 

Open with DEXTER  
In the text 
Fig. 7
Poloidal magnetic field lines. Comparison between two configurations having values of η_{0} larger than that of the reference configuration. 

Open with DEXTER  
In the text 
Fig. 8
Contour plots of B_{φ}. Comparison between two configurations having values of η_{0} larger than that of the reference configuration. 

Open with DEXTER  
In the text 
Fig. 9
Shell average of the toroidal component. Comparison between two configurations having turbulent magnetic diffusivity, η_{0}, larger than that of the reference configuration. The values of η_{0} are 4 × 10^{10} cm^{2} s^{1} for the solid curve and 10^{11} cm^{2} s^{1} for the dashed curve. 

Open with DEXTER  
In the text 
Fig. 10
Poloidal magnetic field lines. Comparison between two configurations with constant η. 

Open with DEXTER  
In the text 
Fig. 11
Shell average of the toroidal component for η = 10^{11} cm^{2} s^{1}. For η = 10^{12} cm^{2} s^{1}, the curve has exactly the same shape, but the values are about 100 times smaller (in absolute value). 

Open with DEXTER  
In the text 
Fig. 12
Contour plots of B_{φ}. Comparison between two configurations with constant η. 

Open with DEXTER  
In the text 
Fig. 13
Comparison of three different profiles of the toroidal field: (1) the solid curve is the shell average for the reference configuration, (2) the dashed curve is the equivalent one for the configuration with constant η (10^{11} cm^{2} s^{1}), (3) the dotted curve is given by the W&C formula with our Ω profile at θ = 88°. All of the curves have been normalised so as to have their peak at 1. 

Open with DEXTER  
In the text 
Fig. 14
Contour plot of ΔΩ B_{θ} for the reference configuration. 

Open with DEXTER  
In the text 
Fig. 15
Radial derivative of the magnetic distortion function D_{m} near the equatorial plane (at θ = 89°) for three configurations: (1) the solid line is for the reference configuration with η_{0} = 10^{10} cm^{2} s^{1}, (2) the dotted line is for the configuration with η_{0} = 4 × 10^{10} cm^{2} s^{1}, (3) the dashed line is for η_{0} = 10^{11} cm^{2} s^{1}. 

Open with DEXTER  
In the text 
Fig. 16
Contour plot of a_{0} in Cartesian coordinates, in the region r = [10−115] r_{g}. This should be compared with Fig. 4. 

Open with DEXTER  
In the text 
Fig. 17
Contour plot of Q for the reference configuration. All locations with Q > 1 are marked in red so as to highlight the regions where Q < 1. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.