© ESO, 2014

1. Introduction

Dipole magnetic fields have two important properties that contribute to their success in the modelling of a pulsar magnetosphere: dipole fields dominate higher order multipole fields at large distances from the neutron star, and they are computationally simpler. Mostly based on the consideration of spin-up lines in the P − Ṗ diagram, Arons (1993) showed that low-altitude magnetic fields of pulsars are dominated by their dipole component, the non-dipole component not exceeding 40% of the dipole field. However, several observations tend to show that multipole magnetic fields cannot be neglected in every pulsar.

Gotthelf et al. (2013) measured period derivatives for the pulsar PSR J0821-4300. It is a central compact object (CCO) in a supernova remnant. They found exceptionally weak dipole magnetic field components for a young neutron star, about 10¹⁰ G. Antipodal surface hot spots with different temperatures and areas were deduced from the X-ray spectrum and pulse profiles. Such non-uniform surface temperature appears to require strong crustal magnetic fields, probably toroidal or quadrupolar components much stronger than the external dipole.

The pulsar J2144-3933, with a period of 8.51 s, is beyond the death-line in the P − B_s diagram and according to standard emission models, the pulsar should not emit radio-waves. Death-line models strongly depend on magnetic field lines curvatures. For a given surface magnetic dipole field strength, pulsars with strong multipolar field components have a highly curved field near the stellar surface that might permit the radio emission of the pulsar J2144-3933 (Young et al. 1999). Indeed, Harding & Muslimov (2011) has shown that a simply offset dipole field can increase the pair-cascade efficiency, and lower the death-line in the P − Ṗ diagram.

Multi-wavelength observations of intermittent radio emissions from rotation-powered pulsars beyond the pair-cascade death line, of the pulse profile of the magnetar SGR 1900+14 after its 1998 August 27 giant flare and of the X-ray spectral features of PSR J0821-4300 and SGR 0418+5729, suggest that the magnetic fields of non-accreting neutron stars are not purely dipolar and may contain higher order multipoles (Mastrano et al. 2013).

Güver et al. (2011) analysed upper bound on the spin-down rate and the high signal-to-noise ratio XMM-Newton spectra of the soft gamma-ray repeater SGR 0418+5729. They found a low surface magnetic field in comparison to other magnetars: 10¹⁴ G. In connection to the spin-down limits, this implies a significantly multipole structure of the magnetic field.

Most attempts to model the pulses profiles of pulsars are based on dipole magnetic fields. But some features of these profiles resist the models. Let us consider, for instance, the brightest pulsar A of the two pulsars binary system PSR J0737-3039. The radio pulse profiles of PSR J0737-3039A consist of two peaks shown in Fig. 1 (Kramer & Stairs 2008). Geometrical models have been produced with best-fit one and two pole models (Ferdman et al. 2013), two poles caustics (TPC), outer gap (OG; Guillemot et al. 2013), and a retarded vacuum dipole polar cap (Perera et al. 2014). With these models, one can reconstruct the main angles defining the orbital plane, rotation axis, and magnetic inclination of the pulsar, as well as the general shape of the pulses. For instance Ferdman et al. (2013) could reconstruct a Gaussian fit, and Perera et al. (2014) considered pulse width at four intensity levels. All these models involve a dipole magnetic field. They found that the two peaks are more likely to be associated with the two poles. But the peaks (especially the less intense one) show sub-structures that do not enter into their models. The sub-structures (a spiky plateau above the 75% intensity level before the main maximum, and a plateau at the 10% level after the main maximum) occupy a significant proportion of the total phase angle. It is quite possible that these sub-structures are associated with multipole components of the electromagnetic environment of the neutron star.

Fig. 1 Pulse profiles of PSR J0737-3039A at various radio frequencies. From Kramer & Stairs (2008).

Another example can be seen in the gamma rays profile of the Vela pulsar revealed by the Fermi-LAT telescope (Abdo et al. 2009) displayed in Fig. 2. Again, this profile contains large sub-structures. It also exhibits shorter sub-structures, which are visible in the enlarged insets. Here again, multipole components might be a cause of pulses sub-structures.

As recalled by Perera et al. (2014), in general, pulsar magnetosphere models are constructed at the following two limits: (a) a vacuum limit (Deutsch 1955); and (b) a force-free magnetohydrodynamics (MHD) limit with a plasma-filled magnetosphere (Spitkovsky 2006). However, a true magnetosphere operates between these two limits. We could expect that the MHD solutions are more realistic, but Harding & Muslimov (2011) found that the rotating dipole magnetosphere in vacuum, in many cases, provides better fits to observed gamma rays pulse profiles than the force-free magnetosphere. This, for instance, is what they found for Vela. This shows that the vacuum magnetosphere is still a useful approximation in pulsar physics.

Considering this general remark and the possible relevance of multipole electromagnetic fields to pulsar models, we present an analytically exact model of the vacuum magnetosphere where the neutron star magnetic field is expanded in multipole components.

Suitable boundary conditions are taken into account for an oblique rotator with a conducting surface. This algorithm allows us to describe electromagnetic fields with l,m quantum numbers as high as 100 (l ≥ m). This algorithm is a generalisation of the one described by Deutsch (1955) for a simple magnetic dipole (l = 1).

In Sect. 2, we present the method of resolution of the Maxwell equations and their boundary conditions. In the Sect. 3, we present the parallel solutions (m = 0) for any value of l. Section 4 contains the general solution of the Maxwell equation with the required boundary conditions for given quantum numbers l,mm ≥ 1. The numbers m> 0 correspond to the perpendicular case, i.e. where the axis of the mutipole is orthogonal to the axis of the neutron star. Thanks to the linearity of the Maxwell equation, the general solution is a linear combination of the perpendicular and parallel solutions. The matching conditions are applied in Sect. 5. Details of the analytical calculations are presented in Appendices A, B.

After this derivation, two applications of multipoles are suggested. The first concerns the problem of the pulsar current closure, and the second concerns the pulse profile of pulsars such as PSR J0737-3039A and Vela.

Fig. 2 Vela broadband (E = 0.1 − 10 GeV) pulse profile. Two pulse periods are shown. The dashed line shows the background level, as estimated from a surrounding annulus during the off-pulse phase. Insets show the pulse shape near the peaks and in the off-pulse region (from Abdo et al. 2009).

2. Methods

Vectors are expanded in spherical coordinates of axis z parallel to the angular velocity vector Ω of the neutron star.

Following Bonazzola et al. (2007) let us define the components of the magnetic field as the usual radial component B_r, and two scalar fields, η and μ, such that (1)The magnetic fields, and μ and η, are also related through the relations (2)The magnetic field associated with μ, noted B^TM is transverse/toroidal (i.e. ) as well as the electric derived from η, noted E^TE (i.e. ), defining the poloidal electromagnetic field. The magnetic and electric fields B, E around a magnetised spinning star B, can be considered the sum of the poloidal field (B^TE, E^TE) and a toroidal field (B^TM, E^TM), (3)The vectors B and E must be solutions of the Maxwell equations in the vacuum (4)The Maxwell Eqs. (4) expressed in term of B_r and of the coefficients μ and η become

where (8)is the angular Laplacian. When it is time dependent, the electric field is deduced from η and μ through the Faraday equation and the relations (9)After separation of the variables, it is found that the angular solution of Eqs. (5), (6) can be expanded in spherical harmonic functions , where are the associated Legendre functions. The scalars r, θ, φ, are the spherical coordinates the z axis being the star spin axis.

Two cases must be treated separately, depending on m. When m = 0, the solution is axially symmetric, and not time dependent. The parts depending on r in Eqs. (5)–(7) are simple differential equations with elementary solutions. When m ≠ 0, the solution is time dependent, and the parts of Eqs. (5)−(7) that depend on r can be converted into Bessel equations of the normalised variable x = mωr/c.

We solve the TM and TE solutions separately. The TE solution is derived from η and has a finite radial magnetic component B_r given by Eq. (5). This equation is solved directly (see the following sections). Then, using the divergence of the magnetic field (10)and the fact that with a angular dependence , we find η. Then, from Eq. (1), (11)The TM magnetic field is derived from μ. The field μ is found directly by resolution of Eq. (6).

The computation of the electric field is different in the cases m = 0 and m ≠ 0, which is detailed in Sects. 3 and 4. The outgoing solution of the Maxwell Eqs. (4) must also satisfy the boundary conditions (BC) (12)at the surface of the star, where n is the unit vector orthogonal to the surface of the star, and R is the radial vector connecting the centre of the star to the point of interest on its surface.

Let be R the radius of the spherical neutron star (NS). Inside of the NS r ≤ R, the magnetic field is generated by internal currents. Let be the l,m component of the spectral decomposition of B. The electric field E^< inside the NS is (13)The matching conditions are (see Eq. (12)) (14)\label{Brmatching}and (15)where E^> is the field in the vacuum.

3. Axially symmetric solutions and their matching conditions

In this section, we compute the multipole electromagnetic field around a rotating neutron star satisfying axially symmetric BC. In terms of spherical harmonics, they correspond to m = 0. When m = 0, there is a finite TM solution derived from Eq. (6), but the curl of this magnetic field is finite too. This means that there is either a time varying electric field, or an electric current density. Because m = 0, a time varying electric field is discarded. Since we are looking for a vacuum solution, a current density is discarded too. Therefore, only a TE electromagnetic field is retained in the axially symmetric case m = 0.

Following the method exposed in Sect. 2, it is found that the components of the vacuum TE magnetic field are: (16)where and is the Legendre polynomial of order l. If the interior of the rotating star is a perfect conductor, the internal electric field E vanishes in the co-rotating frame. Consequently, the electric field in the inertial frame is E^< = (Ω ∧ r) ∧ B and (17)The values of , and at the surface of the star r = R determine the boundary condition for the external field. Outside the star ( r ≥ R), the electric field must be the gradient of an harmonic potential Φ (no charge in the vacuum, steady magnetic field) (18)and its component E_θ must match the components E_θ inside the star: (19)This imposes a series of constraints on the coefficients C_l′. The relation (20)is deduced from the derivative of Eq. (8.914.2) in Gradshteyn et al. (2007) and the differential equation defining the Legendre functions (Eq. (8.820), same reference). Equation (20) is used to deduce the values of the C_l′ coefficients from Eq. (19). Finally, (21)We have added in Φ and E_r the effect of a possible global electric charge Q of the NS.

4. The non-axially symmetric solutions

The solutions corresponding to magnetic fields with an inclination i = 90° over the z axis correspond to m ≠ 0. They are developed in this section.

The TE solution includes a magnetic field with a finite radial component B_r. The solution of Eq. (5) is (22)where Ω = ∥ Ω ∥, c is the light velocity,* , and is the associated Legendre polynomial of order l,m. The function h_l(x) is the spherical Hankel function (23)where H⁽1)_{l + 1/2}(x) is the Bessel function of semi-integer order l + 1/2 and (24) and are constant numbers. Because the solutions involving are associated with an incoming wave, we do not keep them, and for simplicity, we use the notation h_l for .

For the TM solution, μ is derived from Eq. (6). Considering only a single l,m term, (25)where is a constant number. The electric field is derived from the Faraday equation and Eqs. (9). The solutions for the TM and TE components of the electric and magnetic field are:

(26) (27) (28) (29)

5. Matching conditions for the non-axially symmetric solutions

Let be the l,m component of the internal field at the surface of the star. Taking into account the elementary expression of the external magnetic field given by Eq. (26), the matching conditions described by Eq. (2) determine the coefficient in Eq. (26). We have (30)where x_s = mΩR/c. Note that the above B.C. is not sufficient to determine the magnetic field uniquely: in fact, an arbitrary toroidal component B^TM defined by can be added to the poloidal component in a such a way that the electric counterpart E^TM allows us to satisfy the boundary conditions

This equation reads

(31)

The details of its resolution are given in Appendix A. Only two coefficients remain in the right-hand side of Eq. (31); they are (32)and (33)where the coefficient D_l is defined in Eq. (A.6). Finally, the r and θ component of the total electric field ℰ_lm are: (34)The electric field computed above, satisfies the boundary condition ℰ_φ,lm(R,θ,φ,t) = 0. It is shown in Appendix B that it also fits the boundary condition given by Eq. (14).

The magnetic counterpart ℬ is given by (see Eqs. (26), (28)) (35)For l = 1, m = 1 we obtain the result given in Deutsch (1955).

6. A pulsar that extracts electrons from one pole and protons from the other

With dipole pulsar magnetosphere, the open field lines above the two opposite poles present vertical electric fields and Goldreich-Julian currents of the same sign. Therefore, the particles that are extracted from the two poles of the neutron star have the same electric charge. With pulsar dipole magnetosphere model ending with a wind, there is a continuous flux of emitted particles, and it is necessary to close the currents, otherwise the neutron star would accumulate electric charges. Charge accumulation cannot be indefinite, and it is generally assumed that the wind particles (of both positive and negative charges) come from pair creations. The pairs need a continuous flux of primary particles, however, and the question of charge neutrality, i.e. current closure, remains with the primary particles.

Static pulsar electrospheres (Pétri et al. 2002b) are models that do not involve charge circulation. Unfortunately, they do not create a wind either, and they are not expected to radiate. Aligned electrospheres have a dome of charged particles of one sign above each pole, and an equatorial belt of particles of the opposite sign. In that configuration, a dicotron instability can develop. The dicotron effect tends to modulate the shape of the equatorial belt, and it can expel some of its particles (Pétri et al. 2002a). Then, particles of the two signs can be ejected from the neutron star, and this solves the problem of charge neutrality and current closure.

In the present section, we present an alternative to electrospheres and dicotron instability that solves the charge neutrality problem. It consists of a neutron star with a multipole magnetic field. For simplicity, we consider only an aligned dipole and a quadrupole component.

Fig. 3 Normalised radial electric field Er as a function of the normalised electric charge Q at the north pole (thick continuous line) and south pole (thick dashed line). The grey area represent the domain where particles of opposite charges can be extracted from opposite poles.

Fig. 4 A dipole field , and a quadrupole field , and a total electric charge Q = 0.5. The colour code represents the radial electric field Er plotted on the NS surface (within the circle that delimits the surface) and in a meridian plane perpendicular to the line of sight (outside the circle that delimits the NS surface). Magnetic field lines with a foot on the surface in the same meridian plane are plotted as well.

Fig. 5 Goldreich-Julian density nGJ with the same mutipole components as in Fig. 4.

We have (36)(37)where α characterises the quadrupole component amplitude.

The radial electric field is (38)where and Q is an integration constant depending on the total charge of the star.

In what follows, we consider a dynamical process: the electrons are extracted and accelerated from the north pole () and at time t = 0, Q(t) = 0.

Electrons are accelerated above the star surface and they create electron-positron pairs. The vacuum electric field is then progressively screened by the pairs. If the protons remain attached in the vicinity of the star, the magnetosphere charges as long as the electrons are extracted and accelerated. Then, the total electric charge Q of the star increases. Figure 3 illustrates the evolution of the radial electric field at the two poles. It is shown that if Q increases (beware of signs, the normalised charge decreases) the radial electric fields on the north pole is less negative, and those on the south pole becomes more positive. Provided that α> 5/2, a finite range of values of Q (highlighted by a grey rectangle) allows for radial electric fields of opposite signs at the opposite poles.

Figure 4 shows a numerical example of superimposed and aligned dipole and quadrupole fields. The only finite multipole components are characterised by the coefficients (here purely real) , and and the total electric charge is Q = 0.5 C. The electromagnetic field is computed on a spherical grid extending from the star surface to a distance of 716.2 star radii (15 light cylinder radii). The figure only represents the area very close to the star, where the quadrupole component is noticeable. The magnetic field has the intensity 10⁵ T on the surface, and the dipole angle with the rotation axis is null. The period of rotation of the star is 10 ms, which corresponds to a rotation frequency 628 s^-1 and to a light cylinder radius 0.47 × 10⁶ m. We can see (colour code) the radial electric field E_r on the left-hand side as well as magnetic field lines. Because of the quadrupole component, the radial electric field does not have the same value on the two poles. Its high negative value on the north pole is appropriate for the acceleration of electrons out of the star. On the north pole, the positive electric field can accelerate positive ions.

In comparison to Q = 0 (not shown on a figure), the radial electric field amplitude with Q = 0.5 is reduced (but still negative) in the north pole and more positive in the south pole, where protons can be accelerated. Then, the ability of the proton to create pairs is increased, while those of the electron to create pairs remains high. When pairs are created above the two opposite poles, a stationary regime is attained where both electrons and protons are extracted from the star, the total charge reaches an asymptotic value Q₀, and the pulsar can be active. Of course, in this regime, and especially if the NS surface is hot, the Goldreich-Julian density is an important marker of primary charge extraction. We can see in Fig. 5 that the Goldreich-Julian density n_GJ = ∇·⁽B × V_Ω⁾/ 4πc also has opposite signs at opposite poles (V_Ω is the rotational velocity).

Of course, above the pair creation fronts, the electromagnetic field cannot correspond to the vacuum model derived in this paper. But this model is useful below the pair creation front, where the flux of primary particles is not expected to induce currents that could significantly change the magnetic field topology.

With this example, we do not argue that multipole fields are the most common solution to the pulsar current closure problem, but they represent at least one possibility.

At the opposite limit to vacuum approximation, the force-free equations of a magnetosphere were solved in a way that resolves the current system closure. This was done in 2D for an axially symmetric pulsar magnetosphere (Contopoulos et al. 1999; Gruzinov 2005, 2007) and for a 3D dissipative force-free magnetosphere where the magnetic axis is not necessarily aligned with the rotation axis (Spitkovsky 2006; Kalapotharakos & Contopoulos 2009). Those models are based on a dipole magnetic field at the NS surface. The current closes through an equatorial current sheet where the current is opposite to that carried in the open field lines regions. Since force-free equations do not include the plasma transport equations (no explicit equation of density and momentum, for instance), the force-free models do not say much about the nature of the particles that carry currents. It is generally argued that the equatorial return current is carried by electrons that were launched in open field line regions, as well as by positrons moving to the opposite direction, which result from pair creation cascades initiated by primary accelerated electrons.

At a large distance from the NS, a vacuum solution associated with a multipolar electromagnetic field is not different from that associated with a dipole field. This probably holds with a plasma filled magnetosphere. Force-free magnetosphere associated with a surface multipole field might be very analogous to those with dipole fields at distances larger than a fraction of the light-cylinder radius, but the current sheet could be different near the star. As we will see in the next section, this can affect the pulse shape.

7. Pulse shape

In the standard model of the magnetosphere, the strong electric field at the surface of the star r = R extract and accelerate electrons from the crust at relativistic energy. The current density is J ~ en_GJc where n_GJ is the Goldreich-Julian density. Those primary electrons follow the lines of the magnetic field B radiate high-energy γ rays via curvature radiation, and the gamma rays produce electron positron pair via the magnetic field B, if B is strong enough, or by γ rays and crust thermal background x-ray mechanism. Electron positron pairs are supposed to generate the observed radio high energy emission. The main consequence of this mechanism is that the observed pulse shapes depends strongly on the Goldreich-Julian density at the surface of the star.

As mentioned in the introduction, the most often invoked heuristic models to explain pulse shapes are the polar cap model, the slot gap and caustics models, and the outer gap model. In most of these models, a critical area where n_GJ determines the pulse shape is the curve drawn on the NS surface that corresponds to the feet of the last open field lines. Figure 6 shows n_GJ at the NS surface at the feet of the last open magnetic field lines for a dipole field with an inclination i = 40 deg. Its variations are very simple, symmetric, with a single maximum and a single minimum. Multipole components are now added to this dipole field. Their coefficients are displayed in Table 1. The Goldreich-Julian density n_GJ on the NS surface is displayed in Fig. 7. We can see the inclined dipole structure and the superimposition of smaller scale structures with a significant azimuthal modulation. The line corresponding to the feet of the last open field lines is displayed (for the northern hemisphere). The values of n_GJ are displayed in Fig. 8 as a function of the abscissa along the line. We can see that it is more complex than the dipole profile in Fig. 6. The curve has secondary extrema and it shows a higher range of values. Without entering into the detail of pulse shape theories (it is not the topics of the present paper), we can expect that the multipole field can be associated with irregular pulse shapes like those displayed in Figs. 1 and 2.

Fig. 6 Values of the Goldreich-Julian density on the NS surface at the foot of the last open magnetic field lines for a dipole magnetic field of inclination i = 40 deg.

Fig. 7 Values of the Goldreich-Julian density on the NS surface for the multipole magnetic field described in Table 1. The line drawn on the star surface corresponds to the foot (on northern hemisphere) of the last open field lines.

Fig. 8 Values of the Goldreich-Julian density on the NS surface at the foot (on northern hemisphere) of the last open magnetic field lines for the same multipole as in Fig. 7.

8. Conclusion

We have developed an analytical formalism allowing us the most general solution for an electromagnetic field in vacuum fulfilling the boundary conditions on the surface of a rotating magnetised star. This solution, based on an expansion on spherical harmonics is the linear combination of two types of contributions: axially symmetric fields (azimuthal number m = 0) given by Eqs. (16), (21), and non-axially symmetric fields (m ≠ 0) given by Eqs. (34), (35).

Of course, NS are well known to extract plasma in their immediate vicinity, therefore this solution cannot be used as is. Nevertheless, we showed in Sect. 6 that the presence of a quadrupole component of the magnetic field can solve the problem of the current closure in the pulsar magnetosphere. As suggested in Sect. 7, this formalism can also be useful in modelling the observed pulse shapes in pulsars emission.

This solution can be used as a benchmark for codes solving the electromagnetic field equations in the surrounding of a rotating magnetised star.

Pétri (2013) has built numerical solutions of the electromagnetic field surrounding a star with a dipole field in the context of general relativity. Our model does not include gravitational effects, but it is possible that a numerical solution can be developed as well. The present solution can be used as a test when strong gravitational effects are neglected.

Moreover, the vacuum electromagnetic solution can be the first step in an iterative process to find more suitable pulsar models, where a plasma is (numerically) progressively introduced into the NS environment.

References

Appendix A: Derivation of the coefficients

The coefficients are computed by taking the proprieties of the associated Legendre functions into account. We obtain: (Eq. (8.733,1) in Gradshteyn et al. 2007) (A.1)or equivalently (A.2)By using the expression (Eq. (8.731,2) in Gradshteyn et al. 2007) (A.3)and Eq. (A.2) reads (A.4)By replacing the above value of in Eq. (31) we obtain

(A.5)where (A.6)By multiplying both sides of Eq. (A.5) by , after the integration on θ between 0 and π, and on φ between 0 and 2π, with Eq. (A.4) and the orthogonality properties of the associated Legendre functions , (A.7)only two coefficients survive,we obtain Eq. (32). Equation (33) is derived in an analogous way.

Appendix B: Proof that the last boundary condition is fulfilled

We obtained a solution that fulfils the condition E_φ(R) = 0. Does it fit the last condition imposed by the boundary condition E_θ(R) = −(rΩ /c)B_rsinθ? From Eqs. (26)−(30), this requirement is equivalent to

(B.1)When the coefficients and are expressed using Eqs. (32) and (33), the requirement becomes Considering the derivative of Eq. (A.4) relatively to θ, the condition becomes (B.2)By definition, the Lagrange polynomials are the solutions of the differential equation (B.3)With x = cosθ, this differential equation results in the nullity of the expression in Eq. (B.2). This proves that the boundary condition E_θ(R) = −(rΩ /c)B_rsinθ is fulfilled, and that the electromagnetic field derived in Sect. 5 is a consistent solution of the problem.

All Tables

All Figures

	Fig. 1 Pulse profiles of PSR J0737-3039A at various radio frequencies. From Kramer & Stairs (2008).
In the text

	Fig. 2 Vela broadband (E = 0.1 − 10 GeV) pulse profile. Two pulse periods are shown. The dashed line shows the background level, as estimated from a surrounding annulus during the off-pulse phase. Insets show the pulse shape near the peaks and in the off-pulse region (from Abdo et al. 2009).
In the text

	Fig. 3 Normalised radial electric field Er as a function of the normalised electric charge Q at the north pole (thick continuous line) and south pole (thick dashed line). The grey area represent the domain where particles of opposite charges can be extracted from opposite poles.
In the text

Fig. 4 A dipole field , and a quadrupole field , and a total electric charge Q = 0.5. The colour code represents the radial electric field Er plotted on the NS surface (within the circle that delimits the surface) and in a meridian plane perpendicular to the line of sight (outside the circle that delimits the NS surface). Magnetic field lines with a foot on the surface in the same meridian plane are plotted as well.

In the text

	Fig. 5 Goldreich-Julian density nGJ with the same mutipole components as in Fig. 4.
In the text

	Fig. 6 Values of the Goldreich-Julian density on the NS surface at the foot of the last open magnetic field lines for a dipole magnetic field of inclination i = 40 deg.
In the text

	Fig. 7 Values of the Goldreich-Julian density on the NS surface for the multipole magnetic field described in Table 1. The line drawn on the star surface corresponds to the foot (on northern hemisphere) of the last open field lines.
In the text

	Fig. 8 Values of the Goldreich-Julian density on the NS surface at the foot (on northern hemisphere) of the last open magnetic field lines for the same multipole as in Fig. 7.
In the text