Issue |
A&A
Volume 521, October 2010
|
|
---|---|---|
Article Number | A31 | |
Number of page(s) | 10 | |
Section | Astrophysical processes | |
DOI | https://doi.org/10.1051/0004-6361/200913867 | |
Published online | 15 October 2010 |
An investigation of magnetic field distortions in accretion discs around neutron stars
I. Analysis of the poloidal field component![[*]](/icons/foot_motif.png)
L. Naso1,2 - J. C. Miller1,3
1 - SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy
2 - Department of Physics, University of Padova, via Marzolo 8, 35131 Padova, Italy
3 - Department of Physics (Astrophysics), University of Oxford, Keble Road, Oxford OX1 3RH, UK
Received 14 December 2009 / Accepted 7 June 2010
Abstract
We 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.
Key words: accretion, accretion disks - magnetic fields - magnetohydrodynamics (MHD) - turbulence - methods: numerical - X-rays: binaries
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: X-ray pulsars and old neutron stars in the process of being spun-up (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. 1977; Ghosh & Lamb 1979a, hereafter referred to as GL; Ghosh & Lamb 1979b), 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 1987).
Wang (1987) and Campbell (1987, 1992) 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
(
)
and that the magnetosphere is
corotating with the star (with angular velocity
), they then
obtained
.
Miller & Stone (1997) performed some numerical simulations with a 2D
model, including the effects of non-zero 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.
1977) and also that the poloidal component can be very different from a
dipolar field. Agapitou & Papaloizou (2000) 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 non-zero).
Many investigators have simulated the magnetosphere-disc interaction, solving the full set of the MHD equations (see e.g. Romanova et al. 2002; Kulkarni & Romanova 2008). 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 large-scale 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 large-scale simulations. The conceptual papers from the 70 s and 80 s mentioned above, continue to be widely quoted and used as the basis for new research (see e.g. Kluzniak & Rappaport 2007) 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 non-zero also along the other
directions. We use spherical coordinates (r, ,
), 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 (1997).
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 Alfvén 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 Alfvén 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 h/r to be constant (where h is the semi-thickness
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. (2000).
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
.
As regards the radial component of the velocity, vr, we use the expression
given for the middle region of -discs
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 vr 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.
![]() |
Figure 1:
Schematic representation of our model (not to scale). We use
|
Open with DEXTER |
3 Equations
If one considers mean fields, the induction equation can be written as:![]() |
(1) |
where






The induction equation then reduces to:
where

After imposing stationarity (i.e.
)
and
axisymmetry (i.e.
), the three components of
Eq. (2) become:
The first two of these equations have the same expression inside the large square brackets and we call this

where k is a generic constant.
We notice that
contains only the turbulent magnetic
diffusivity
,
the poloidal components of the velocity field and the
r 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,
Eq. (7).
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. (7) 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:
With this definition, the axisymmetric field




where k is the constant introduced in Eq. (6) and vr,




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
(1973):
where the radius r is expressed in units of the Schwarzschild radius, m is given in solar mass units,




![]() |
(11) |
For the other component of the poloidal velocity,

![]() |
(12) |
For the diffusivity, as mentioned in Sect. 2, we are using a value


where
where







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:
Following a procedure similar to that used for obtaining Eqs. (6) and (7), we go to spherical coordinates, write out the three component equations and group the poloidal terms. This gives:
where
is a generic constant. In order to calculate it we
consider a path with constant
,
e.g.
;
along this path Eq. (17) gives:
where all of the terms in the brackets are constant. We do not know the exact profiles of vr and


Equation (20) implies not only that if there is a non-zero radial velocity then there must be a non-zero vertical velocity as well, but also that the vertical velocity is larger than the radial one calculated at the same location (for


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 r 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 (1987, 1992) and Wang (1987) 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 non-zero 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. (7) to find the toroidal component of the
magnetic field.
We note that Eq. (16) has been solved in the context
of stellar winds by Mestel (1961); his result was that the poloidal
magnetic field and velocity field need to be parallel. Our result that
(from Eq. (17) with
)
is in agreement with this.
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,


where we go from (22) to (23) by dividing both sides by



where
![$[v_r]=[v_\theta]=$](/articles/aa/full_html/2010/13/aa13867-09/img119.png)
![$[r_{{\rm g}}]=$](/articles/aa/full_html/2010/13/aa13867-09/img120.png)
![$[\eta]=$](/articles/aa/full_html/2010/13/aa13867-09/img121.png)


We have solved Eq. (24) with the Gauss-Seidel 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).
![]() |
Figure 2:
Magnetic field lines from the numerical solution (solid) and those for a dipole (dotted). The top panel refers to a case with no accretion, where the field remains exactly dipolar, while in the bottom panel v0=104 cm s-1, and the field is distorted. The diffusivity |
Open with DEXTER |
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).
![]() |
Figure 3:
The same as in Fig. 2, but with different values of v0. The top panel is for a case with v0=105 cm s-1 while the bottom one is for v0=106 cm s-1. If we use
|
Open with DEXTER |
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
10-14 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 v0, i.e. for different accretion
rates. 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 vr 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 in-between 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 vr).
In fact the magnetic Reynolds number
,
which describes the
general solution of the induction, Eq. (2), is built from
them: it is defined as
,
where
l0, v0 and
are respectively a characteristic length,
velocity and diffusivity. This parameter gives the relative importance of
the two terms on the right-hand 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 vr 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 l0 to be the radius of the inner edge of the disc, v0
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 l0, 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
.
![]() |
Figure 4:
Comparison between the numerical |
Open with DEXTER |
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 right-hand 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 (










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,
l0. In addition there is a quite natural way for choosing l0. By
looking at Eq. (24) one can see that, if we choose
l0 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
x). We can then think of l0 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 (1997). 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
.
![]() |
Figure 5:
Top panel: the magnetic distortion function in the equatorial
plane. Bottom panel: |
Open with DEXTER |
![]() |
Figure 6:
The radial derivative of the magnetic distortion function in the
equatorial plane. The profiles refer to the same case as in Fig. 5. Comparing with the bottom panel of Fig. 5 one can see that the peaks in the magnetic field
|
Open with DEXTER |
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 back-reaction 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. 2000), 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 (1979a)
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 80%, 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 (1997) (compareour 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. (7) 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 back-reaction on the plasma flow and the dynamo effect.
AcknowledgementsIt 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
- Agapitou, V., & Papaloizou, J. C. B. 2000, MNRAS, 317, 273 [NASA ADS] [CrossRef] [Google Scholar]
- Campbell, C. G. 1987, MNRAS, 229, 405 [NASA ADS] [Google Scholar]
- Campbell, C. G. 1992, GApFD, 63, 179 [Google Scholar]
- Ghosh, P., & Lamb, F. K. 1979a, ApJ, 232, 259 [NASA ADS] [CrossRef] [Google Scholar]
- Ghosh, P., & Lamb, F. K. 1979b, ApJ, 234, 296 [NASA ADS] [CrossRef] [Google Scholar]
- Ghosh, P., Lamb, F. K., & Pethick, C. J. 1977, ApJ, 217, 578 [NASA ADS] [CrossRef] [Google Scholar]
- Kluzniak, W., & Rappaport, S. 2007, ApJ, 671, 1990 [NASA ADS] [CrossRef] [Google Scholar]
- Kulkarni, A. K., & Romanova, M. M. 2008, MNRAS, 386, 673 [NASA ADS] [CrossRef] [Google Scholar]
- Mestel, L. 1961, MNRAS, 122, 473 [NASA ADS] [Google Scholar]
- Miller, K. A., & Stone, J. M. 1997, ApJ, 489, 890 [NASA ADS] [CrossRef] [Google Scholar]
- Romanova, M. M., Ustyugova, G. V., Koldoba, A. V. & Lovelace, R. V. E. 2002, ApJ, 578, 420 [NASA ADS] [CrossRef] [Google Scholar]
- Rekowski, M. V., Ruediger, G., & Elstner, D. 2000, A&A, 353, 813 [NASA ADS] [Google Scholar]
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
- Wang, Y.-M. 1987, A&A, 183, 257 [NASA ADS] [Google Scholar]
Online Material
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
A.1 Description of the code
In order to solve the 2D elliptic PDE (24) we use the
Gauss-Seidel relaxation method. If we call the elliptic operator
and the right hand side b, 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 b=0.
We approximate the operators by discretizing the functions over a grid. The scheme that we use for discretizing the derivatives is as follows:
where the indices i and j refer to grid points along the x and ycoordinate directions respectively, while t represents the pseudo-time or iteration step. We use expressions (A.1)-(A.3) to discretize Eq. (24) and then, isolating the term

We solve this proceeding from i=1, j=1; on the right hand side the terms that have already been calculated (i.e. the terms at positions i-1 and j-1) are taken at the current iteration t+1. Once vr,



The magnitude of the central dipole field and the accretion rate do not
enter Eq. (A.4) 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 10 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-1 which is the value obtained from Eq. (10) when
,
, 0.03) or (0.1, 0.07);
- diffusivity:
,
,
cm2 s-1 and
cm2 s-1;
- initial estimate for the magnetic stream function:
(for a dipolar field
);
- iteration time step:
.
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 vr,
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 vr 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 (Ni and Nj) and the number of iterations.
We use five grids in total. When testing the dependence on Nj we use:
,
and
;
while when testing
the dependence on Ni 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
15%, while for the
grid it is
and for the
grid it is
;
(2) changing Njdoes not produce any visible effect: while increasing Nj by a factor of 4
(from 20 to 80) decreases the maximum relative error only slightly (
7%), changing Ni
from 100 to 400 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 Ni, 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 Fig. A.1.
![]() |
Figure A.1:
Maximum analytic error for the test D2 with three grids, which differ only in Ni. Increasing Nj from 20 to 80 produces a very small decrease in the error, only |
Open with DEXTER |
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



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.
![]() |
Figure A.2:
The rms of the numerical solution for different values of |
Open with DEXTER |
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
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 r3 (this gives Br
and
increasing linearly with r); 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 0 for case
(4) to 220 for cases (1) and (2) and to 226 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 5%.
Finally we considered varying the iteration step size, i.e. the
in Eq. (A.4), that is written as
.
There is no simple argument of principle that can be
used to determine the best value for c, therefore we determined it
experimentally. We considered the same configuration as in test D2 and ran
it several times varying only the value of c, going from 0.025
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 c 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 n always smaller than 1. We find that using the value n=0.95 is a good compromise in minimizing the
number of iterations and preserving the code stability.
Footnotes
- ... component
- Appendix A is only available in electronic form at http://www.aanda.org
- ... Keplerian
- Campbell also considered non-Keplerian flow in the inner part of the disc.
- ... disappears
- It is only
the turbulent diffusivity
which disappears, the Ohmic one
will always be present, therefore in vacuum
.
- ...
-discs
- For the parameter
values which we are using, both
and
lie within the ``middle region'' of the
-disc model.
- ...
rates
- For the Shakura-Sunyaev model, the radial velocity is
proportional to
if the mass of the central object and the viscosity
are kept fixed.
All Figures
![]() |
Figure 1:
Schematic representation of our model (not to scale). We use
|
Open with DEXTER | |
In the text |
![]() |
Figure 2:
Magnetic field lines from the numerical solution (solid) and those for a dipole (dotted). The top panel refers to a case with no accretion, where the field remains exactly dipolar, while in the bottom panel v0=104 cm s-1, and the field is distorted. The diffusivity |
Open with DEXTER | |
In the text |
![]() |
Figure 3:
The same as in Fig. 2, but with different values of v0. The top panel is for a case with v0=105 cm s-1 while the bottom one is for v0=106 cm s-1. If we use
|
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Comparison between the numerical |
Open with DEXTER | |
In the text |
![]() |
Figure 5:
Top panel: the magnetic distortion function in the equatorial
plane. Bottom panel: |
Open with DEXTER | |
In the text |
![]() |
Figure 6:
The radial derivative of the magnetic distortion function in the
equatorial plane. The profiles refer to the same case as in Fig. 5. Comparing with the bottom panel of Fig. 5 one can see that the peaks in the magnetic field
|
Open with DEXTER | |
In the text |
![]() |
Figure A.1:
Maximum analytic error for the test D2 with three grids, which differ only in Ni. Increasing Nj from 20 to 80 produces a very small decrease in the error, only |
Open with DEXTER | |
In the text |
![]() |
Figure A.2:
The rms of the numerical solution for different values of |
Open with DEXTER | |
In the text |
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.