An investigation of magnetic field distortions
in accretion discs around neutron stars
Key Words.:
accretion, accretion disks – magnetic fields – magnetohydrodynamics (MHD) – turbulence – methods: numerical – Xrays: binariesWe report results from calculations investigating stationary magnetic field configurations in accretion discs around magnetised neutron stars. Our strategy is to start with a very simple model and then progressively improve it, providing complementary insight into results obtained with large numerical simulations. In our first model, presented here, we work in the kinematic approximation and consider the stellar magnetic field as being a dipole aligned with the stellar rotation axis and perpendicular to the disc plane, while the flow in the disc is taken to be steady and axisymmetric. The behaviour in the radial direction is then independent of that in the azimuthal direction. We investigate the distortion of the field caused by interaction with the disc matter, solving the induction equation numerically in full 2D. The influence of turbulent diffusivity and fluid velocity on the poloidal field configuration is analysed, including discussion of outflows from the top and bottom of the disc. We find that the distortions increase with increasing magnetic Reynolds number (calculated using the radial velocity). However, a single global parameter does not give an adequate description in different parts of the disc and we use instead a ‘magnetic distortion function’ (a magnetic Reynolds number defined locally). Where (near to the inner edge of the disc) there is little distortion, but where (most of the rest of the disc), there is considerable distortion and the field becomes weaker than the dipole would have been. Between these two regions, there is a transition zone where the field is amplified and can have a local minimum and maximum. The location of this zone depends sensitively on the diffusivity. The results depend very little on the boundary conditions at the top of the disc.
1 Introduction
With this paper, we begin a study of the properties of accretion discs around magnetised neutron stars. We are interested in two kinds of system: Xray pulsars and old neutron stars in the process of being spunup (recycled) to become millisecond pulsars. Here we focus mostly on these recycled pulsars, for which the distortions of the magnetic field are larger and occur at smaller distances from the central object, making the effects easier to see.
The study of accretion onto magnetised neutron stars was pioneered by Ghosh and Lamb who, in a series of 3 papers (Ghosh et al. GL77 (); Ghosh & Lamb GL79a (), hereafter referred to as GL; Ghosh & Lamb GL79b ()), investigated the flow of accreting matter and the magnetic field configuration. They were able to estimate the total torque exerted on the neutron star by the accretion disc but, for doing this, they needed to make a number of quite drastic assumptions. In particular they used a dipolar profile for the poloidal component of the magnetic field and an ad hoc prescription for the azimuthal component, which was later shown to be inconsistent with having a stable disc beyond the corotation point (Wang W87 ()).
Wang (W87 ()) and Campbell (C87 (), C92 ()) obtained an analytic expression for the azimuthal component of the magnetic field, making the assumption that the poloidal component is dipolar. They solved the induction equation to find a steady solution for an axisymmetric magnetic field and found that the generation of toroidal field is due only to the vertical gradient of the azimuthal velocity . By further assuming that the disc is Keplerian^{1}^{1}1Campbell also considered nonKeplerian flow in the inner part of the disc. () and that the magnetosphere is corotating with the star (with angular velocity ), they then obtained .
Miller & Stone (MS97 ()) performed some numerical simulations with a 2D model, including the effects of nonzero resistivity. Their simulations confirmed that there is not a total screening of the magnetic field from the accretion disc (as had initially been assumed by Ghosh et al. GL77 ()) and also that the poloidal component can be very different from a dipolar field. Agapitou & Papaloizou (AP00 ()) also found that the poloidal component of the field can differ significantly from the dipole field of the central star, and that the magnetic torque can be much smaller than that estimated assuming .
Our aim here is to improve on these analyses and to find a stationary configuration for the magnetic field inside the disc and in the surrounding corona, without making any leading order expansion or vertical integration in the induction equation, without making any assumptions about the poloidal component of the magnetic field inside the disc and by using a more general profile of the velocity field (with all of the components being allowed to be nonzero).
Many investigators have simulated the magnetosphere–disc interaction, solving the full set of the MHD equations (see e.g. Romanova et al. Rom02 (); Kulkarni & Romanova KR08 ()). The work that we are presenting here should not be seen as being in competition with these analyses, but rather as being complementary to them, by providing both a useful test case and a means to better understand the underlying processes. Our strategy is to use a succession of simplified models (of which this paper presents the initial one) becoming progressively more sophisticated, and to proceed step by step so as to fully understand the effect of each successive additional feature as it is introduced. In largescale numerical works, one sees the results of an interacting set of inputs within the scope of the adopted model assumptions and numerical techniques. Deconstructing this, so as to have a clear conceptual understanding of the role of each of the different components, remains a valuable thing to do and an approach which needs to be carried on alongside the largescale simulations. The conceptual papers from the 70s and 80s mentioned above, continue to be widely quoted and used as the basis for new research (see e.g. Kluzniak & Rappaport KR07 ()) and our work here stands in the tradition of refining these approaches.
The structure of this paper is as follows. After the present Introduction, the details of the model are given in Sect. 2, and the equations are presented in Sect. 3. As we will show there, it is possible to separate the calculation for the poloidal component of the magnetic field from that for the toroidal component. In this paper we present results for the poloidal field; the toroidal one will be discussed in a subsequent paper. The numerical code used is described briefly in Sect. 3, and then in more detail in the Appendix (only in the online version of the paper). In Sect. 4 we present our results and Sect. 5 contains conclusions.
2 Model
In this study we consider disc accretion by a neutron star having a dipolar magnetic field. We assume that the star is rotating around its magnetic axis, and that this axis is perpendicular to the disc plane; 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. Turbulent diffusivity is included and the velocity field is not constrained to be purely rotational but it is allowed to be nonzero also along the other directions. We use spherical coordinates (, , ), with the origin being at the centre of the neutron star.
We assume that at large radii the magnetic pressure is negligible with respect to the gas pressure and that the disc can be described there as a standard disc. As one moves inwards, the magnetic field becomes progressively stronger and, for the cases considered here, eventually dominates over the gas pressure. It is convenient to divide the disc into two regions: an outer zone, where the magnetic field does not greatly influence the flow, and an inner zone, where the magnetic stresses are dominant over the viscous ones and matter begins to leave the disc vertically. For a sufficiently strong magnetic field, most of the matter eventually leaves the disc following magnetic field lines and the disc is disrupted. We note that some equatorial accretion continues to be possible even with a rather high stellar magnetic field, as shown by Miller & Stone (MS97 ()).
The precise location of the inner edge of the disc, , is still an open issue: several prescriptions have been suggested, but none of them differs very much from the Alfven radius calculated using a dipolar magnetic field. For our present model we take this as giving the inner edge of the disc; for subsequent models, a possible improvement will be to calculate the Alfven radius using the stationary magnetic field obtained in this work instead of the dipolar one.
We note here that, if we take the innermost part of the disc to have very high diffusivity, we naturally obtain the same structure as that of the GL model, i.e. a narrow inner region followed by a broad outer one. However the behaviour of the field in these regions is quite different from that in the GL model, as will be shown in the results section.
We suppose that all around the pulsar there is vacuum, except for where we have the disc and the corona, and that the field remains dipolar everywhere in this vacuum. In reality, this is not completely correct, both because the density in the magnetosphere is not zero and because between the star and the disc there is the matter which is accreting onto the neutron star. For the former we are introducing a low density corona in order to have a transition zone between the disc and the vacuum. We also allow the velocity and the diffusivity to have different values in the corona from those in the disc. For the latter, we suppose that the matter flow is perfectly aligned with field lines and that it has very low density, so that we can neglect its influence on the magnetic field structure. We are aware that in real astrophysical systems, imposing the dipole condition at the boundaries is rather drastic since the field is distorted well before reaching the disc. However at our level of analysis, the results are not very sensitively dependent on this, and we are here focusing on studying the influence on the magnetic configuration of the velocity field and diffusivity. That this is reasonable is confirmed by the fact that the magnetic configuration which we obtain is rather similar to that given by the numerical simulations of other authors (Miller & Stone 1997), which had a different treatment of the boundary conditions. We will comment more on this in Sect. 4.
Summarizing, we divide the surroundings of the central object into four parts (see Fig. 1): (1) the inner disc, extending from out to a transition radius (where the diffusivity changes); (2) the outer disc, extending from to an outer radius ; (3) a corona, above and below the disc; and (4) everything else, which we take here to be vacuum. Within the outer disc, we focus on the main region, extending from out to the light cylinder at ; the outer edge of our grid, , is very much further away than this, so that conditions imposed there do not affect the solution in the region of interest.
We take the ratio to be constant (where is the semithickness of the disc), and so the entire upper surface of the disc is located at a single value of (as also is the case for the corona). We take an opening angle of for the disc (measuring from the equatorial plane to the top of the disc), and for the disc plus corona. The disc then has .
Once the magnetic diffusivity and the poloidal velocity are specified, the poloidal component of the induction equation can be solved to obtain the configuration of the poloidal field without entering into details of the toroidal component. As regards the turbulent magnetic diffusivity, we take this to have a constant value in the bulk of the outer disc, and a higher value in the corona, the inner disc and near to the outer boundary (we join the different parts smoothly, using error functions). A similar approach was used by Rekowski et al. (R00 ()). We take higher values in regions (1) and (3) because the density is lower there and so we expect the effects of turbulence to be enhanced. Moreover in both regions we expect the field not to be influenced very much by the plasma flow and having a larger diffusivity is a way to reduce this influence. However as we move away from the corona into the vacuum, the density drops to zero and turbulence eventually disappears^{2}^{2}2It is only the turbulent diffusivity which disappears, the Ohmic one will always be present, therefore in vacuum ..
As regards the radial component of the velocity, , we use the expression given for the middle region of discs ^{3}^{3}3 For the parameter values which we are using, both and lie within the “middle region” of the disc model. by Shakura & Sunyaev (1973). In Sect. 3.2, we find that whenever a dipole field is a solution of the induction equation, a precise relation must hold between and . We use this relation to calculate in the corona, so as to be consistent with the dipole boundary conditions, and put to zero inside the disc.
3 Equations
If one considers mean fields, the induction equation can be written as:
(1) 
where 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 , where the term generates the socalled effect. In the present paper we neglect this effect, which could however be included in a subsequent dynamo model. This is in line with our approach, which is aiming at understanding, one at a time, the effects of the various elements which characterize the system of an accretion disc around a magnetised neutron star.
The induction equation then reduces to:
(2) 
where , because the turbulent diffusivity is much larger than the Ohmic one.
After imposing stationarity (i.e ) and axisymmetry (i.e ), the three components of Eq. (2) become:
(3)  
(4)  
The first two of these equations have the same expression inside the large square brackets and we call this . We can then reduce these two equations to a single equation and write the entire system as:
(6)  
where is a generic constant.
We notice that contains only the turbulent magnetic diffusivity , the poloidal components of the velocity field and the and components of the magnetic field. It is clear therefore that Eq. (6) on its own (plus Maxwell’s equation ) governs the poloidal part of the magnetic field and is independent of any azimuthal quantity, while to obtain the toroidal component of the magnetic field one has to solve the further equation (3).
In this paper we concentrate only on solving for the poloidal field, using Eq. (6), for which no knowledge of behaviour in the direction is required. In a forthcoming paper we will use the results obtained here to solve Eq. (3) and calculate the toroidal field component.
In order to guarantee that is satisfied and to be able to calculate the magnetic field lines easily, we write the magnetic field components in terms of the magnetic stream function , which is implicitly defined by:
(8) 
With this definition, the axisymmetric field is always solenoidal and isolines of represent magnetic field lines. Substituting these expressions into Eq. (6), we obtain an elliptic partial differential equation for in terms of and :
(9) 
where is the constant introduced in Eq. (6) and , and are nonconstant coefficients. This equation can be solved once boundary conditions and values for the coefficients have been specified. As described in Sect. 2 we are interested in a configuration where the magnetic field is a pure dipole at the boundaries. This implies and , where is the magnetic dipole moment.
3.1 Velocity field and turbulent diffusivity
In Sect. 2, we justified and qualitatively described the profiles of velocity and diffusivity that we are using. Here we give the precise expressions.
For the velocity, we use the expression given by Shakura & Sunyaev (SS73 ()):
(10) 
where the radius is expressed in units of the Schwarzschild radius, is given in solar mass units, is in units of the critical Eddington rate and is the standard ShakuraSunyaev viscosity coefficient. Using typical values (, and ) one obtains:
(11) 
For the other component of the poloidal velocity, , we set this to zero in the disc and use a nonzero profile in the corona. The formula for this is calculated in Sect. 3.2; we anticipate here the result:
(12) 
For the diffusivity, as mentioned in Sect. 2, we are using a value in the bulk of the outer disc, and a higher value, , in the corona, the inner disc and near to the outer boundary. We use combinations of error functions for giving a smooth join between the different regions:
(13) 
where
(14)  
(15) 
where on the upper surface of the disc, is the radius of the boundary between the inner disc and the main region, and is located shortly before , far away from the zone of interest. The widths of the error functions used for the angular and radial profiles are given by and respectively. We have used and .
3.2 Dipolar solution, an analytic constraint
In order to obtain a profile for we consider the situation when, from the top surface of the corona (at ) down to some , the stationary magnetic field is dipolar, i.e. where is the magnetic dipole moment. Then Eq. (2) becomes:
(16) 
Following a procedure similar to that used for obtaining Eqs. (6) and (3), we go to spherical coordinates, write out the three component equations and group the poloidal terms. This gives:
(17)  
(18) 
where is a generic constant. In order to calculate it we consider a path with constant , e.g. ; along this path equation (17) gives:
(19) 
where all of the terms in the brackets are constant. We do not know the exact profiles of and , but it is not plausible that the left hand side increases as and so we need to choose . Therefore from the last equation we get:
(20) 
Equation (20) implies not only that if there is a nonzero radial velocity then there must be a nonzero vertical velocity as well, but also that the vertical velocity is larger than the radial one calculated at the same location (for one has ).
The behaviour of the azimuthal velocity can be investigated using Eq. (18). As expected is a possible solution but is not, meaning that having a dipolar field is not consistent with having a Keplerian angular velocity profile, whereas it is consistent with having no rotation at all. It is interesting that there are also some non trivial profiles which are solutions. For example , which is equivalent to , gives a set of possible solutions. Therefore if one supposes that decreases with as a power law, then it must have also a dependence on in order for the magnetic field to be dipolar. We recall that in the works of Campbell (C87 (), C92 ()) and Wang (W87 ()) it is precisely the vertical gradient of the angular velocity that produces the toroidal field. Here we have shown that it is possible to have a nonzero vertical gradient of the angular velocity and still have a zero toroidal magnetic field. However, we stress that we are not giving physical explanations for having these kinds of velocity profiles. We will comment further on the profile of the angular velocity in the forthcoming paper where we will solve Eq. (3) to find the toroidal component of the magnetic field.
3.3 Solution method
Before solving Eq. (9) we put it into a dimensionless form, by scaling the quantities in the following way:
(21) 
where the hat quantities are dimensionless, is the Schwarzschild radius and is a reference value for the stream function, calculated as the value for a dipolar field on the equator of the neutron star. Substituting into Eq. (9) with we get:
(22) 
(23) 
where we go from (22) to (23) by dividing both sides by . We rename the variables ( and ) and obtain the following dimensionless form for the equation:
(24) 
where cm s, cm, cm s, so that and are dimensionless coefficients.
We have solved Eq. (24) with the GaussSeidel relaxation method, which uses a finite difference technique, approximating the operators by discretizing the functions over a grid. At any given iteration step, the values of the stream function 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 (only in the online version of the paper).
Before applying the numerical scheme to the real 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 for the velocity and diffusivity, initial estimates for the stream function, locations for the outer radial boundary of the grid and values for the iteration time step. In this way we have checked the stability and convergence, have optimized the iteration procedure and have determined where to place the outer radial boundary of the grid (which needs to be far enough away so that the outer boundary conditions do not significantly influence the solution in our region of interest).
All of the results presented in the next section have been obtained using a grid and with total iterations. With these settings, we have always obtained residuals of the order of or less, and the resolution is and .
4 Results
In this section we describe how the magnetic field configuration changes when we modify the velocity field and the turbulent diffusivity.
In Figs. 2 and 3 we show the magnetic field lines calculated with four values of , i.e. for different accretion rates^{4}^{4}4For the ShakuraSunyaev model, the radial velocity is proportional to if the mass of the central object and the viscosity are kept fixed.. For facilitating the comparison we show also a dipolar magnetic field (dotted). The field lines are labelled with the radial coordinate where the dipole field imposed at the top boundary would cross the equatorial plane if not distorted. We can see that if were zero, the field would not be distorted at all from the dipolar configuration and increasing the velocity then creates progressively more distortion. The degree of distortion depends on the location in the disc: in the inner part, where the field is strongest, it is most able to resist distortion; further out, the field is weaker and it becomes progressively more distorted.
According to the behaviour of the magnetic field lines, we can divide the disc into three regions: (1) an inner region, where the lines are not distorted very much away from the dipole; (2) an outer region, where the distortion can be very large and (3) the region inbetween the two, which we call a transition region, where there is an accumulation of field lines. Consequently in the transition region there is an amplification of the magnetic field (see Fig. 4).
In addition to varying the radial velocity, we also considered the role of the diffusivity. The results show that when we change , we get the opposite behaviour to that seen when varying the velocity, i.e. a larger gives a smaller distortion. Actually what really matters is not the velocity or the diffusivity alone but their ratio. This is not surprising since in the equation which we are solving (Eq. (24)) the quantities only appear in this ratio (bearing in mind that is taken to be either zero or proportional to ). In fact the magnetic Reynolds number , which describes the general solution of the induction equation (2), is built from them: it is defined as , where , and are respectively a characteristic length, velocity and diffusivity. This parameter gives the relative importance of the two terms on the righthand side of the induction equation. For large , we are in the regime of ideal MHD with the magnetic field and plasma being frozen together; for low , instead, the field and plasma are almost decoupled and the field simply diffuses. In accretion discs the radial velocity is usually many orders of magnitude smaller than the azimuthal velocity. In our case the Reynolds numbers calculated using the two velocities differ by about five orders of magnitude, if one takes the Keplerian velocity as the characteristic velocity. For our present calculations however only the radial motion is relevant, since in the disc and is proportional to in the corona, while does not appear in the equation that we are solving now. The value of reported in Figs. 2 and 3 is therefore the one calculated taking the characteristic velocity to be the radial velocity.
The panels in these figures are clearly showing that the distortion of the field is proportional to the magnetic Reynolds number calculated in this way. This happens because with increasing the freezing condition gets progressively stronger so that any fluid motion perpendicular to the magnetic field lines encounters more and more resistance. Therefore, since the velocity field is fixed, the magnetic field has to change. Figures 2 and 3 not only show that modifications in the magnetic field lines increase with , but also that their shape is consistent with what is expected from considering the flux freezing condition in the case of a conical flow (which is what we have in the disc).
However the actual value of the magnetic Reynolds number is somewhat arbitrary, because in general there is no unambiguous way of defining the characteristic length, velocity and diffusivity of a given system. In our case we choose to be the radius of the inner edge of the disc, to be the radial velocity at the inner edge of the disc and to be the value of the diffusivity in the main disc region. One could also make a different choice for the characteristic length , such as taking this to be the radius of the star or the average height of the disc; the trend of having larger distortions for larger values of would of course be seen in all cases, but the switching on of the distortions would occur at different threshold values of .
We have already noted that the distortion varies with position, and so it is clear that a single global parameter cannot give a sufficiently detailed description in all parts of the system. It is therefore convenient to introduce a new quantity which we call the “magnetic distortion function” . We define this in the same way as the magnetic Reynolds number but, instead of taking characteristic values for the velocity and diffusivity, we take the local values:
(25) 
This function gives the relative importance of the two terms on the righthand side of the induction equation at every point of the disc, rather than giving just a global measure as with the standard magnetic Reynolds number. We then expect the advection term () to dominate if and the diffusion term () to dominate if . This then explains why we find three regions inside the disc: the inner region corresponds to the zone where , the main region to , and the transition region to intermediate values of . This correspondence is made clear in Fig. 4, where we show the component of the magnetic field, the dipolar profile and the magnetic distortion function, all on the equatorial plane. We recall that the jump in follows from the profile chosen for , i.e. we use a larger value of the diffusivity in the inner part of the disc.
Another important aspect of the magnetic distortion function is that its definition is less arbitrary than that for the standard magnetic Reynolds number, since it is defined using only one characteristic value, . In addition there is a quite natural way for choosing . By looking at Eq. (24) one can see that, if we choose to be equal to our unit of length, , then the magnetic distortion function is already there in the equation (it is the coefficient of the partial derivative of with respect to ). We can then think of as a quantity needed to make the ratio dimensionless, and the most natural choice for this is the characteristic scale being used as the unit length.
Summarizing, we can describe the magnetic field configuration in the accretion disc by saying that magnetic field lines that enter the disc in the main region () are pushed towards the central object, whereas those which enter the disc in the inner region () are almost unmodified. The result is that in between these two regions there is an accumulation of field lines, and so there is an amplification of the magnetic field there, as can be seen in Fig. 4.
In order to test the dependence of the results on the boundary conditions, we have experimented with several different profiles for the magnetic stream function outside the disc. We have used three additional profiles: one which gives a magnetic field with spherical field lines, one which gives a magnetic field with vertical field lines and another one which gives field lines inclined at an arbitrary angle to the vertical axis. We have chosen these profiles by comparison with the results from the simulations by Miller & Stone (MS97 ()). For magnetic field lines entering the disc at the same locations, their shape within the disc varies hardly at all in the different cases (although the actual value of the field strength can be different). We conclude that the distortion of the field lines does not depend sensitively on the boundary conditions used and is instead mainly governed by the magnetic distortion function .
In order to better understand the influence of the magnetic distortion function on the magnetic field structure, we varied and saw how the field changed. We used three new profiles for and considered the previous one as a reference. In the first profile we increased the value of in the inner disc and left the rest unmodified, in the second one instead we lowered in the outer disc and did not change the inner part, and in the last one we just changed the width of the transition between the low and high values of . We then calculated the poloidal magnetic field and the results are presented in Fig. 5, where the top panel shows the different profiles of the magnetic distortion function and the bottom one shows the component of the magnetic field, referring to the equatorial plane in both cases.
Considering this figure, we can summarize the influence of the magnetic distortion function with four comments: (1) changing the value of in the inner part by a factor of 5 leaves the magnetic field almost unchanged; (2) on the other hand, the magnetic field is very sensitive to the width of the transition in and to its value in the outer disc; in particular (3) the position of the peak in is related to the width of the transition; and (4) the deviations away from the dipole field are mainly governed by the value of in the outer disc. We can go further and consider the radial derivative of , which is shown in Fig. 6 for all of the profiles used. From this we can see that the position of the peak of is strongly connected with the position of the maximum in , and that the maximum amount of magnetic distortion is related to the height of the peak in the derivative of .
5 Conclusions
In this paper we have begun a systematic study of the magnetic field configuration inside accretion discs around magnetised neutron stars, which is intended as being complementary to the large numerical simulations being carried out elsewhere. We have assumed that the star itself has a dipolar magnetic field, whose axis is aligned with the rotation axis, which is perpendicular to the disc plane. We have also assumed that the flow is steady and has axial symmetry everywhere. Our strategy was to start the analysis with a very simple model, where we made the kinematic approximation, included turbulent magnetic diffusivity and solved the induction equation numerically in full 2D, without making any leading order expansion. This initial model will subsequently be improved by including the magnetic backreaction on the fluid flow.
We have shown that it is possible to separate the calculation for the configuration of the poloidal magnetic field from any azimuthal quantities. This is a key point and has the consequence that the effective magnetic Reynolds number is that calculated with the radial velocity rather than the azimuthal one (which is very much larger). We have here considered only the poloidal component; the toroidal one will be addressed in a subsequent paper.
We have modelled the surroundings of the neutron star as being composed of four regions (see Fig. 1): the inner and the outer disc, the corona (taken to be a layer above and below the disc) and all of the rest, which is taken to be vacuum. We suppose that the stellar magnetic field remains dipolar until it reaches the corona. At that point it begins to feel the presence of the fluid flow and the magnetic field lines are pushed inwards, thus creating distortions away from the purely dipolar field.
We have studied the response of the magnetic field to changes in the velocity and the diffusivity, finding distortions away from dipolar increasing with the radial infall velocity and decreasing with increasing diffusivity. The underlying behaviour is that the distortions increase together with the magnetic Reynolds number where the ratio appears.
However a single value of cannot take into account any large changes in the magnitudes of the velocity and the diffusivity through the disc, since it is defined using single characteristic values. Therefore in order to have a sufficiently detailed description of the system, we have introduced a magnetic distortion function , based on local values of the quantities concerned, so that in regions where or one should expect to have large or small distortions respectively. We expect the turbulence to be enhanced in the regions of lower density (the corona and the inner part of the disc), therefore in our model we use a larger value of in these zones (similarly to Rekowski et al. R00 ()), giving a smaller value for . Moreover in these regions the magnetic field should be less sensitive to the plasma flow and, given that we are in the kinematic approximation, the only way of achieving this is to use a larger value of the diffusivity. As clearly shown in the panels of Figs. 2 and 3, the disc can be divided into three parts: (1) an inner region, where and the distortions are negligible; (2) a transition region, where is rapidly increasing and magnetic field lines accumulate; (3) an outer region, in the inner part of which and the the distortions can be very large. Considering can be very useful when analysing results of large numerical simulations.
Comparing our results with previous literature, we can confirm the idea of dividing the disc into two principle regions: an inner part, where the magnetic field is strongest, and an outer part, where the magnetic field is weaker and gently decaying. However the behaviour that we find for the field in these regions is very different from that of the Ghosh & Lamb model (GL79a ()) and we find it convenient to include a third zone, to be considered as a transition between the two principle ones (see Fig. 4). In the inner boundary layer of the GL model, the magnetic field is reduced by screening currents by a factor of , while in our case the field is barely modified in the first region. In the transition region instead we see a new effect: the field is amplified and has a local maximum. The properties of the maximum (i.e. its location, the peak magnitude of the field and its behaviour in the neighbourhood) are well described by the magnetic distortion function , in particular by its maximum value, by the width of the transition between the low and high values and by the behaviour of the radial derivative (i.e. by the location and the magnitude of its maximum). The behaviour in the outer zone is instead quite similar to that in the GL model, with the field decaying and being smaller than the dipole one at any given location.
As regards the magnetic field geometry, our results resemble rather closely those obtained by Miller & Stone (MS97 ()) (compare our Fig. 3 with the top panels of their figure 3), despite the fact that they solved the full set of MHD equations whereas we have solved just the induction equation and with different conditions at the top of the disc. Moreover we have found that the distortion of the field lines inside the disc depends very little on the profiles outside it (i.e. on the boundary conditions). We conclude that the distortion pattern seen for the poloidal component of the magnetic field should be rather a general result, related only to the fundamental aspects of the system, and that the distortion is not at all negligible for typical values of the accretion rate (see Figs. 2 and 3).
In our next paper, we will calculate the toroidal component of the magnetic field by solving Eq. (3) and will study its dependence on the angular velocity and other components of the magnetic field. Afterwards additional elements will be included in the model, including backreaction on the plasma flow and the dynamo effect.
Acknowledgments
It is a pleasure to thank Alfio Bonanno, Claudio Cremaschini and Kostas Glampedakis for stimulating discussions; this work was partly supported by CompStar, a Research Networking Programme of the European Science Foundation.
References
 (2000) Agapitou V. & Papaloizou J.C.B. 2000, MNRAS, 317, 273
 (1987) Campbell C.G. 1987, MNRAS, 229, 405
 (1992) Campbell C.G. 1992, GApFD, 63, 179
 (1977) Ghosh P., Lamb F.K. & Pethick C.J. 1977, ApJ, 217, 578
 (1979a) Ghosh P. & Lamb F.K. 1979a, ApJ, 232, 259
 (1979b) Ghosh P. & Lamb F.K. 1979b, ApJ, 234, 296
 (2007) Kluzniak W. & Rappaport S. 2007, ApJ, 671, 1990
 (2008) Kulkarni A. K. & Romanova M. M. 2008, MNRAS, 386, 673
 (1961) Mestel L. 1961, MNRAS, 122, 473
 (1997) Miller K. A. & Stone J. M. 1997, 489, 890
 (2002) Romanova M. M., Ustyugova G. V., Koldoba A. V. & Lovelace R. V. E. 2002, ApJ, 578, 420
 (2000) Rekowski M.v., Ruediger G. and Elstner D. 2000, A&A, 353, 813
 (1973) Shakura N.I. & Sunyaev R.A. 1973, A&A, 24, 337
 (1987) Wang Y.M. 1987, A&A, 183, 257
Appendix A The Code
In this Appendix we describe the numerical code that we have used to solve Eq. (24) and discuss some of the tests that we have performed on it.
a.1 Description of the code
In order to solve the 2D elliptic PDE (24) we use the GaussSeidel relaxation method. If we call the elliptic operator and the right hand side , then the original equation becomes: . We turn this elliptic equation into a hyperbolic one by adding a pseudo time derivative; we can then consider the iterative procedure as a time evolution and write: . In our case and .
We approximate the operators by discretizing the functions over a grid. The scheme that we use for discretizing the derivatives is as follows:
(26) 
(27) 
(28) 
where the indices and refer to grid points along the and coordinate directions respectively, while represents the pseudotime or iteration step. We use expressions (26)(28) to discretize Eq. (24) and then, isolating the term , we get the iterative algorithm that we use in our code:
(29)  
We solve this proceeding from , ; on the right hand side the terms that have already been calculated (i.e. the terms at positions and ) are taken at the current iteration . Once , and have been specified, a solution can then be obtained for any initial estimate of .
The magnitude of the central dipole field and the accretion rate do not enter Eq. (29) directly, but they are used to calculate the location of the inner edge of the disc .
For the mass and radius of the neutron star, we use the canonical values, and km respectively. We fix the accretion rate as (in units of ), giving a magnetospheric radius of about when G, as typical for a millisecond pulsar.
a.2 Testing of the code
In order to check the code for stability and convergence, to estimate errors and to optimize the iteration procedure by choosing an appropriate iteration step, we performed a number of tests, some of which are now described.
During this test phase we used the following values for the parameters:

magnetic field at the stellar surface: G;

size of the domain: , , ,

and ;

radial velocity at inner edge: cm s

which is the value obtained from Eq. (10) when

, , or , ;

diffusivity: , , cm s and

cm s;

initial estimate for the magnetic stream function:

(for a dipolar field );

iteration time step: .
The tests can be divided into two main groups: with and without a known analytic solution. For the latter, we can estimate errors by calculating the residuals and comparing the solutions obtained with different grid resolutions, while for the former we also have the difference between the computed values and the exact results.
a.2.1 Test with an analytic solution
We consider two cases. The first has dipolar boundary conditions and either no poloidal motion, or a velocity profile consistent with condition (20). In this case the poloidal component of the field must be dipolar everywhere (we refer to this test as D, for dipolar). The second case has the boundary conditions for set to zero. In this configuration, regardless of the profile used for the velocity, is a solution in all of the domain (we call this test Z, for zero). This last test allows us to test the code by including all of the terms that will be present when solving for the cases of interest (i.e. including , and ).
In both cases, we test two different configurations (which we call D1, D2, Z1 and Z2) by changing the velocity profile. In test D1 we set all velocities to zero; in test D2 is zero only in the disc while in the corona it is given by Eq. (10) and is given by Eq. (20). In test Z1 we consider the same velocity profile as the one that we will use for our cases of interest, given by Eqs. (10) and (20) and finally in test Z2 we use the same velocity profile as in test D2.
In all of these tests we follow the same procedure: we verify the stability of the code, we estimate the error and see if it scales correctly, checking the convergence of the solution. We do this by studying how the numerical solution changes when varying the number of grid points ( and ) and the number of iterations.
We use five grids in total. When testing the dependence on we use: , and ; while when testing the dependence on we use: , and . For each of these grids we calculate: (i) the absolute difference and (ii) the relative difference, between the numerical solution and the analytic one at each point of the grid; and (iii) the root mean square (rms) of the numerical solution at each iteration step. The results obtained are very similar for all of the five grids and for each of the four tests and can be summarized with the following four statements: (1) both the absolute error and the relative error have a maximum near to the inner edge and then decrease quite rapidly. For the grid, the maximum relative error is , while for the grid it is and for the grid it is ; (2) changing does not produce any visible effect: while increasing by a factor of 4 (from to ) decreases the maximum relative error only slightly (), changing from to has a much greater effect, giving a decrease in the error of two orders of magnitude; (3) the reduction in the rms and in the maximum error becomes progressively smaller with increasing , thus showing that we have convergence of the numerical solution; (4) using a sparser grid gives smaller errors at the beginning and during the relaxation process, however if one keeps iterating until the saturation level is reached, then the error with sparser grids is larger than with denser grids (suggesting that this problem could be suited for a multigrid approach). Regarding statements (2), (3) and (4), see figure 7.
a.2.2 Test with an unknown solution
Here we use the same configuration as the one that we will use for our cases of interest, i.e. dipolar boundary conditions and velocity field given by Eqs. (10) and (20). Even if we do not know the solution for this setup, we know from Eq. (9) that it has to approach a dipolar field when the coefficients and both go to zero. In order to test this we have considered five configurations, each with a different value of ranging between and cm s. Fig. 8 shows clearly that for increasing , the rms of the numerical solution is approaching that for a dipole calculated on the same grid.
For these five configurations we also performed the tests previously described, i.e. the ones regarding changing the grid and comparing the errors and the rms. The results are again similar and confirm the four statements made earlier.
a.2.3 Other tests
We used the configuration of test D2 for checking three further aspects: (1) determining the importance of the initial estimate for , investigating which values to choose; (2) for the location of the outer radial boundary; and (3) for the iteration step.
The kind of algorithm which we are using to solve Eq. (9) needs an initial estimate for the solution. According to how good or bad this estimate is with respect to the correct solution, one needs a smaller or larger number of iterations for completing the process. In order to show this and also to demonstrate that the final solution does not depend on the initial profile, we used four initial estimates for the magnetic stream function : (1) a constant value; (2) a Gaussian profile (centred on and with a width of ); (3) a profile increasing with (this gives and increasing linearly with ); and (4) the profile which gives a dipolar magnetic field. As expected, in all of the cases the final solution is the same, even for configuration (3), and the number of iterations required to reach saturation changes and goes from for case (4) to for cases (1) and (2) and to for case (3).
As mentioned in Sects. 2 and A.1, our region of physical interest goes from the inner edge of the disc out to the light cylinder . Since we do not want the solution in this region to be influenced by the outer boundary condition, we ran some tests using different values for the radius of the outer edge and then compared the numerical solutions in the region of interest. We used the same setup as for the real problem that we were wanting to solve and six values of (, , , , and ). We found that the difference within the region of interest between the numerical solutions obtained using two subsequent values of became progressively smaller, until one could barely distinguish the different solutions. We decided to put the outer boundary at for the physical analysis; this gives results differing from those with by less than about .
Finally we considered varying the iteration step size, i.e. the in Eq. (29), that is written as . There is no simple argument of principle that can be used to determine the best value for , therefore we determined it experimentally. We considered the same configuration as in test D2 and ran it several times varying only the value of , going from upwards. We found that the final asymptotic error was always the same, but that the number of iterations required to relax to the final solution was changing, decreasing as increased. However there is an upper limit: when the numerical solution diverges. Transferring this condition to , we obtain . We can then change the way in which the iteration step size is calculated in the code and write: , with always smaller than . We find that using the value is a good compromise in minimizing the number of iterations and preserving the code stability.