General solution for the vacuum electromagnetic field in the surroundings of a rotating star.

Many recent observations of pulsars and magnetars can be interpreted in terms of neutron stars (NS) with multipole electromagnetic fields. As a first approximation, we investigate the multipole magnetic and electric fields in the environment of a rotating star when this environment is deprived of plasma. We compute a multipole expansion of the electromagnetic field in vacuum for a given magnetic field on the conducting surface of the rotating star. Then, we consider a few consequences of multipole fields of pulsars. We provide an explicit form of the solution. For each spherical harmonic of the magnetic field, the expansion contains a finite number of terms. A multipole magnetic field can provide an explanation for the stable sub-structures of pulses, and they offer a solution to the problem of current closure in pulsar magnetospheres. This computation generalises the widely used model of a rotating star in vacuum with a dipole field. It can be especially useful as a first approximation to the electromagnetic environment of a compact star, for instance a neutron star, with an arbitrarily magnetic field.

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-\dot{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^{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-\dot{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^{14}$ 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.

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 is for instance 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\geq m$). This algorithm is a generalisation of the one described by Deutsch (1955) for a simple magnetic dipole ($l=1$).

In section 3, we present the method of resolution of the Maxwell equations and their boundary conditions. In the section 4, we present the parallel solutions ($m=0$) for any value of $l$. Section 5 contains the general solution of the Maxwell equation with the required boundary conditions for given quantum numbers $l,m\;m\geq 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 section 6. 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.

Vectors are expanded in spherical coordinates of axis $z$ parallel to the angular velocity vector ${\bf\Omega}$ 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, $\eta$ and $\mu$, such that

	$\displaystyle B_{\theta}$	$\displaystyle=$	$\displaystyle\partial_{\theta}\eta-\frac{1}{\sin\theta}\partial_{\phi}\mu,$
	$\displaystyle B_{\phi}$	$\displaystyle=$	$\displaystyle+\frac{1}{\sin\theta}\partial_{\phi}\eta+\partial_{\theta}\mu.$		(1)

The magnetic fields, and $\mu$ and $\eta$, are also related through the relations

	$\displaystyle\nabla_{\theta\phi}^{2}\eta$	$\displaystyle=$	$\displaystyle\partial_{\theta}B_{\theta}+\frac{\cos\theta}{\sin\theta}B_{\theta}+\frac{1}{\sin\theta}\partial_{\phi}B_{\phi},$
	$\displaystyle\nabla_{\theta\phi}^{2}\mu$	$\displaystyle=$	$\displaystyle\partial_{\theta}B_{\phi}+\frac{\cos\theta}{\sin\theta}B_{\phi}-\frac{1}{\sin\theta}\partial_{\phi}B_{\theta}$		(2)
		$\displaystyle=$	$\displaystyle\frac{1}{\sin\theta}[\partial_{\theta}(B_{\phi}\sin\theta)-\partial_{\phi}B_{\theta}].$

The magnetic field associated with $\mu$, noted $B^{TM}$ is transverse/toroidal (i.e. $B_{r}^{TM}=0$) as well as the electric derived from $\eta$, noted $E^{TE}$ (i.e. $E_{r}^{TE}=0$), defining the poloidal electromagnetic field. The magnetic and electric field ${\bf B}$, ${\bf E}$ around a magnetised spinning star ${\mathbf{B}}$, can be considered the sum of the poloidal field (${\mathbf{B}}^{TE}$, ${\mathbf{E}}^{TE}$) and a toroidal field (${\mathbf{B}}^{TM}$, ${\mathbf{E}}^{TM}$ ),

$${\mathbf{B}}={\mathbf{B}}^{TE}+{\mathbf{B}}^{TM},\qquad\qquad\qquad{\mathbf{E}}={\mathbf{E}}^{TE}+{\mathbf{E}}^{TM}.$$

(3)

The vector ${\vec{B}}$ and ${\vec{E}}$ must be solutions of the Maxwell equations in the vacuum

	$\displaystyle\frac{1}{c}\frac{\partial{\bf B}}{\partial\,t}$	$\displaystyle=$	$\displaystyle-\nabla\wedge\,{\bf E},\;\;\;\;\frac{1}{c}\,\frac{\partial{\bf E}}{\partial\,t}=\nabla\wedge\,{\bf B},$
	$\displaystyle\nabla\cdot{\bf B}$	$\displaystyle=$	$\displaystyle 0,\;\;\;\;\nabla\cdot{\bf E}=0.$		(4)

The Maxwell equations (4) expressed in term of $B_{r}$ and of the coefficients $\mu$ and $\eta$ become

	$\displaystyle(\frac{\partial^{2}}{\partial r^{2}}$	$\displaystyle+$	$\displaystyle\frac{2}{r^{2}}+\frac{4}{r}\frac{\partial}{\partial r})B_{r}+\nabla^{2}_{\theta,\phi}B_{r}-\frac{1}{c^{2}}\frac{\partial}{\partial t}B_{r}=0,$		(5)
	$\displaystyle\partial_{r^{2}}^{2}\mu$	$\displaystyle+$	$\displaystyle\frac{2}{r}\partial_{r}\mu+\frac{1}{r^{2}}\nabla^{2}_{\theta,\phi}\mu+k^{2}\mu=\mu_{S}$		(6)
	$\displaystyle\partial_{r^{2}}^{2}\eta$	$\displaystyle+$	$\displaystyle\frac{2}{r}\partial_{r}\eta+\frac{1}{r^{2}}\nabla^{2}_{\theta,\phi}\eta+k^{2}\eta+\frac{2}{r^{2}}B_{r}=\eta_{S}$		(7)

where

	$\displaystyle\nabla^{2}_{\theta,\phi}$	$\displaystyle=$	$\displaystyle\partial_{\theta^{2}}^{2}+\cot\theta\partial_{\theta}+\frac{1}{\sin^{2}\theta}\partial_{\phi^{2}}^{2}$		(8)
		$\displaystyle=$	$\displaystyle\frac{1}{\sin\theta}[\partial_{\theta}(\sin\theta\,\partial_{\theta})+\frac{1}{\sin\theta}\partial_{\phi^{2}}^{2}]$

is the angular Laplacian. When it is time dependent, the electric field is deduced from $\eta$ and $\mu$ through the Faraday equation and the relations

	$\displaystyle(\nabla\times\vec{B})_{r}$	$\displaystyle=$	$\displaystyle\frac{1}{r}\nabla^{2}_{\theta,\phi}\mu$
	$\displaystyle(\nabla\times\vec{B})_{\theta}$	$\displaystyle=$	$\displaystyle\frac{1}{r\sin\theta}[\partial_{\phi}B_{r}-(I+r\partial_{r})\partial_{\phi}\eta-\sin\theta(I+r\partial_{r})\partial_{\theta}\mu]$
	$\displaystyle(\nabla\times\vec{B})_{\phi}$	$\displaystyle=$	$\displaystyle\frac{1}{r\sin\theta}[-\sin\theta\partial_{\theta}B_{r}+\sin\theta(I+r\partial_{r})\partial_{\theta}\eta$		(9)
			$\displaystyle-(I+r\partial_{r})\partial_{\phi}\mu].$

After separation of the variables, it is found that the angular solution of Eqs. 5-6 can be expanded in spherical harmonic functions $Y^{m}_{l}(\theta,\phi)=P_{l}^{m}(\cos\theta)e^{im\phi}$, where $P^{m}_{l}(\cos\theta)$ are the associated Legendre functions. The scalars $r$, $\theta$, $\phi$, 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\neq 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\omega r/c$.

We solve the $TM$ and $TE$ solutions separately. The $TE$ solution is derived from $\eta$ 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

$$\nabla\cdot\vec{B}=\frac{\partial}{\partial r}B_{r}+\frac{2}{r}B_{r}+\frac{1}{r}\nabla^{2}_{\theta,\phi}\eta,$$

(10)

and the fact that with a $Y_{l}^{m}$ angular dependence $\nabla^{2}_{\theta,\phi}\eta=-l(l+1)\eta$, we find $\eta$. Then, from Eq. (3),

	$\displaystyle B_{\theta}$	$\displaystyle=$	$\displaystyle\frac{1}{l(l+1)}\left(x\frac{\partial^{2}B_{r}}{\partial x\partial\theta}+2\frac{\partial B_{r}}{\partial\theta}\right),$
	$\displaystyle B_{\phi}$	$\displaystyle=$	$\displaystyle\frac{im}{l(l+1)\sin\theta}\left(x\frac{\partial^{2}B_{r}}{\partial x\partial\phi}+2\frac{\partial B_{r}}{\partial\phi}\right).$		(11)

The $TM$ magnetic field is derived from $\mu$. The field $\mu$ is found directly by resolution of Eq. (6).

The computation of the electric field is different in the cases $m=0$ and $m\neq 0$, which is detailed in sections 4 and 5. The outgoing solution of the Maxwell Eqs.(4) must also satisfy the boundary conditions (BC)

$${\bf E}\times{\bf n}=\left[\frac{1}{c}({\bf\Omega}\times{\bf R})\wedge{\bf B}\right]\times{\bf n}$$

(12)

at the surface of the star, where ${\bf n}$ is the unit vector orthogonal to the surface of the star, and $\bf{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\leq R$, the magnetic field is generated by internal currents. Let be ${\bf B}^{<}_{rlm}(r,\theta,\phi)$ the $l,m$ component of the spectral decomposition of ${\bf B}$. The electric field ${\bf E^{<}}$ inside the NS is

$$E^{<}_{r}=\Omega\frac{r}{c}B^{<}_{\theta}\sin\theta,\;\;E^{<}_{\theta}=-\Omega\frac{r}{c}B^{<}_{r}\sin\theta,\;\;E^{<}_{\phi}=0.$$

(13)

The matching conditions are (See Eq.(12))

$$B^{>}_{r}(R,\theta,\varphi)=B^{<}_{r}(R,\theta,\varphi),\;\;E^{>}_{\theta}(R,\theta,\varphi)=E^{<}_{\theta}(R,\theta,\varphi),$$

(14)

and

$$E^{>}_{\phi}(R,\theta,\varphi)=0,$$

(15)

where ${\bf E}^{>}$ is the field in the vacuum.

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 section 3, it is found that the components of the vacuum $TE$ magnetic field are:

	$\displaystyle B_{r}$	$\displaystyle=$	$\displaystyle B_{l}^{0}\left(\frac{R}{r}\right)^{l+2}\,P^{0}_{l}(\theta),$
	$\displaystyle B_{\theta}$	$\displaystyle=$	$\displaystyle-B_{l}^{0}\frac{1}{l+1}\left(\frac{R}{r}\right)^{l+2}\frac{d\,P^{0}_{l}(\theta)}{d\,\theta},$
	$\displaystyle B_{\phi}$	$\displaystyle=$	$\displaystyle 0,$		(16)

where $P^{0}_{l}(\theta)=\mathcal{P}^{0}_{l}(\cos\theta)$ and $\mathcal{P}^{0}_{l}$ is the Legendre polynomial of order $l$. If the interior of the rotating star is a perfect conductor, the internal electric field $\vec{E}$ vanishes in the co-rotating frame. Consequently, the electric field in the inertial frame is $\vec{E}^{<}=(\vec{\Omega}\land\vec{r})\land\vec{B}$ and

$$E_{r}^{<}=+\Omega\frac{r}{c}\sin\theta\,B_{\theta}^{<},\;\;\;E_{\theta}^{<}=-\Omega\frac{r}{c}\sin\theta B_{r}^{<}.$$

(17)

The value of $B^{<}_{r}(R)$ , $E^{<}_{\theta}(R)$ and $E^{<}_{\phi}(r)$ at the surface of the star $r=R$ determine the boundary condition for the external field. Outside the star ( $r\geq R$), the electric field must be the gradient of an harmonic potential $\Phi$ (No charge in the vacuum, steady magnetic field)

$$\Phi=\sum_{l^{\prime}}\frac{C_{l^{\prime}}}{r^{(l^{\prime}+1)}}P_{l^{\prime}}^{0}(\theta)$$

(18)

and its component $E_{\theta}$ must match the components $E_{\theta}$ inside the star:

$$E_{\theta}(R)=\frac{1}{R}\frac{\partial\,\Phi}{\partial\,\theta}=-\Omega\frac{R}{c}\sin\theta\,B_{l0}^{0}\,P^{0}_{l}.$$

(19)

This imposes a series of constraints on the coefficients $C_{l^{\prime}}$. The relation

$$\frac{d}{d\theta}\left[P^{0}_{l+1}(\theta)-P^{0}_{l-1}(\theta)\right]=-(2l+1)P^{0}_{l}(\theta)\sin\theta$$

(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^{\prime}}$ coefficients from Eq. (19). Finally,

	$\displaystyle\Phi$	$\displaystyle=$	$\displaystyle\frac{B_{l0}^{0}}{2\,l+1}\left[\left(\frac{R}{r}\right)^{l+2}P_{l+1}^{0}(\theta)-\left(\frac{R}{r}\right)^{l}P_{l-1}^{0}(\theta)\right]\Omega\frac{R^{2}}{c}-\frac{Q}{r}$
	$\displaystyle E_{r}$	$\displaystyle=$	$\displaystyle=\frac{B_{l0}^{0}}{2l+1}\left[(l+2)P^{0}_{l+1}\left(\frac{R}{r}\right)^{l+3}-l\,P^{0}_{l-1}\left(\frac{R}{r}\right)^{l+1}\right]\Omega\frac{R}{c}+\frac{Q}{r^{2}}$
	$\displaystyle E_{\theta}$	$\displaystyle=$	$\displaystyle\frac{-B_{l0}^{0}}{2\,l+1}\left[\frac{dP^{0}_{l+1}}{d\theta}\left(\frac{R}{r}\right)^{l+3}-\frac{dP^{0}_{l-1}}{d\theta}\,\left(\frac{R}{r}\right)^{l+1}\right]\Omega\frac{R}{c}$
	$\displaystyle E_{\phi}$	$\displaystyle=$	$\displaystyle 0.$		(21)

We have added in $\Phi$ and $E_{r}$ the effect of a possible global electric charge $Q$ of the NS.

The solutions corresponding to magnetic fields with an inclination $i=90^{o}$ over the $z$ axis correspond to $m\neq 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

$$B_{r}=\sum_{l=1}^{\infty}\sum_{-l\leq m\leq l}\left[C_{lm}^{(1)}\frac{h_{l}^{(1)}(x)}{x}+C_{lm}^{(2)}\frac{h_{l}^{(2)}(x)}{x}\right]P_{lm}(\cos\theta)e^{im(\phi-\Omega t)},$$

(22)

where $\Omega=\|{\bf\Omega}\|$, $c$ is the light velocity,* $P_{l}^{m}(\theta)=\mathcal{P}_{l}^{m}(\cos\theta)$, and $\mathcal{P}_{l}^{m}$ is the associated Legendre polynomial of order $l,m$. The function $h_{l}(x)$ is the spherical Hankel function

$$h_{l}(x)=\sqrt{\frac{\pi}{2x}}H^{1}_{l+1/2}(x)$$

(23)

where $H^{(}1)_{l+1/2}(x)$ is the Bessel function of semi-integer order $l+1/2$ and

$$x=m\,\frac{r\,\Omega}{c},\;\;\;\varphi=\phi-\Omega\,t.$$

(24)

$C_{lm}^{(1)}$ and $C_{lm}^{(2)}$ are constant numbers. Because the solutions involving $h_{l}^{2}$ are associated with an incoming wave, we do not keep them, and for simplicity, we use the notation $h_{l}$ for $h_{l}^{(1)}$.

For the $TM$ solution, $\mu$ is derived from Eq. (6). Considering only a single $l,m$ term,

$$\mu=\sigma_{l,m}C_{lm}^{\prime}h_{l}(x)Y_{lm}(\theta,\phi),$$

(25)

where $C_{lm}^{\prime}$ 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:

	$\displaystyle B_{r;\,lm}^{TE}(r,\theta,\phi,t)=A^{TE}_{lm}\ \frac{h_{l}(x)}{x}\ P_{l}^{\,m}(\theta)\ e^{im\varphi}$		(26)
	$\displaystyle B_{\theta;\,lm}^{TE}(r,\theta,\phi,t)=\frac{A^{TE}_{lm}}{l(l+1)}\left(\frac{1}{x}\frac{d}{dx}\left(x\,h_{l}(x)\right)\right)\frac{dP_{l}^{\,m}(\theta)}{d\theta}\ e^{im\varphi}$
	$\displaystyle B_{\phi;\,lm}^{TE}(r,\theta,\phi,t)=i\,\frac{mA^{TE}_{lm}}{l(l+1)}\left(\frac{1}{x}\frac{d}{dx}\left(x\,h_{l}(x)\right)\right)\frac{P_{l}^{\,m}(\theta)}{\sin\theta}\ e^{im\varphi}$

	$\displaystyle E_{r;\,lm}^{TE}(r,\theta,\phi,t)=0$		(27)
	$\displaystyle E_{\theta;\,lm}^{TE}(r,\theta,\phi,t)=-m\,A^{TE}_{lm}\,\frac{h_{l}(x)}{l(l+1)}\,\frac{P_{l}^{\,m}(\theta)}{\sin\theta}\ e^{im\varphi}$
	$\displaystyle E_{\phi;\,lm}^{TE}(r,\theta,\phi,t)=-i\,A^{TE}_{lm}\frac{h_{l}(x)}{l(l+1)}\,\frac{dP_{l}^{\,m}(\theta)}{d\theta}\ e^{im\varphi}$

	$\displaystyle B_{r;\,lm}^{TM}(r,\theta,\phi,t)=0$		(28)
	$\displaystyle B_{\theta;\,lm}^{TM}(r,\theta,\phi,t)=+mA^{TM}_{lm}\,\frac{h_{l}(x)}{l(l+1)}\,\frac{P_{l}^{\,m}(\theta)}{\sin\theta}\ e^{im\varphi}$
	$\displaystyle B_{\phi;\,lm}^{TM}(r,\theta,\phi,t)=+i\,A^{TM}_{lm}\ \frac{h_{l}(x)}{l(l+1)}\ \frac{dP_{l}^{\,m}(\theta)}{d\theta}\ e^{im\varphi}$

	$\displaystyle E_{r;\,lm}^{TM}(r,\theta,\phi,t)=\,A_{lm}^{TM}\frac{h_{l}(x)}{x}\,P_{l}^{\,m}(\cos\theta)\,e^{im\varphi}$		(29)
	$\displaystyle E^{TM}_{\theta;\,lm}(r,\theta,\phi,t)=\frac{A^{TM}_{lm}}{l(l+1)}\left(\frac{1}{x}\frac{d}{dx}\left(x\,h_{l}(x)\right)\right)\frac{dP_{l}^{\,m}(\theta)}{d\theta}\,e^{im\varphi}$
	$\displaystyle E^{TM}_{\phi;\,lm}(r,\theta,\phi,t)=i\,\frac{m\,A^{TM}_{lm}}{l(l+1)}\!\left(\frac{1}{x}\frac{d}{dx}\left(x\,h_{l}(x)\right)\right)\!\!\frac{P_{l}^{\,m}(\theta)}{\sin\theta}\,e^{im\varphi}.$

Let be $B^{<}_{r;\,lm}(R)P^{m}_{l}(\theta)e(im\varphi)$ the $l,m$ component of the internal field $B_{r}^{<}$ 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.(14) determine the coefficient $A^{TE}_{lm}$ in Eq.(26). We have

$$A^{TE}_{lm}=B^{<}_{r;\,lm}(R)\frac{x_{s}}{h_{l}(x_{s})}$$

(30)

where $x_{s}=m\Omega R/c$. Note that the above B.C. is not sufficient to determine the magnetic field uniquely: in fact, an arbitrary toroidal component ${\bf B}^{TM}$ defined by $B_{r}^{TM}=0$ can be added to the poloidal component in a such a way that the electric counterpart ${\bf E}^{TM}$ allows us to satisfy the boundary conditions

$$E^{TE}_{\phi}(R,\theta,\phi,t)+E^{TM}_{\phi}(R,\theta,\phi.t)=0.$$

This equation reads

	$\displaystyle{\cal E}_{\phi}(R,\theta,\phi,t)=\sum_{lm}-A^{TE}_{lm}\,\frac{h_{l}(x_{s})}{l(l+1)}\,\frac{dP_{l}^{m}(\theta)}{d\theta}\ e^{im\varphi}$		(31)
	$\displaystyle+\sum_{l^{\prime}m^{\prime}}\ \frac{m^{\prime}\,A^{TM}_{l^{\prime}m^{\prime}}}{l^{\prime}(l^{\prime}+1)}\left(\frac{1}{x}\frac{d}{dx}\left(x\,h_{l^{\prime}}(x)\right)\right)_{r=R}\!\frac{P_{l^{\prime}}^{m^{\prime}}(\theta)}{\sin\theta}\,e^{im^{\prime}\varphi}\ =\ 0.$

The details of its resolution are given in Appendix A. Only two coefficients $A^{TM}_{l^{\prime}m^{\prime}}$ remain in the right-hand side of Eq. (31); they are

$$A^{TM}_{l+1\,m}=A^{TE}_{l\,m}\frac{h_{l}(x_{s})}{D_{l+1}}\frac{(l-m+1)(l+2)}{m(2l+1)}$$

(32)

and

$$A^{TM}_{l-1\,m}=-A^{TE}_{lm}\,\frac{h_{l}(x_{s})}{D_{l-1}}\frac{(l+m)(l-1)}{m(2l+1)},$$

(33)

where the coefficient $D_{l}$ is defined in Eq. (44). Finally, the $r$ and $\theta$ component of the total electric field ${\bf{\cal E}}_{lm}$ are:

	$\displaystyle{\cal E}_{r,lm}(r,\theta,\phi)$	$\displaystyle=$	$\displaystyle A^{TM}_{l+1\,m}\frac{h_{l+1}(x)}{x}P^{m}_{l+1}(\theta)e^{i\,m\varphi}$
			$\displaystyle+A^{TM}_{l-1\,m}\frac{h_{l-1}(x)}{x}\,P^{m}_{l-1}(\theta)\,e^{i\,m\varphi},$

	$\displaystyle{\cal E}_{\theta,lm}$	$\displaystyle=$	$\displaystyle-mA^{TE}_{lm}\frac{h_{l}(x)}{l(l+1)}\frac{P^{m}_{l}(\theta)}{\sin\theta}\,e^{im\varphi}$
			$\displaystyle+\frac{A^{TM}_{l+1\,m}}{(l+1)(l+2)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l+1}(x))\right)\frac{dP^{m}_{l+1}(\theta)}{d\theta}\,e^{i\,m\varphi}$
			$\displaystyle+\frac{A^{TM}_{l-1\,m}}{l(l-1)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l-1}(x))\right)\frac{dP^{m}_{l-1}(\theta)}{d\theta}\,e^{im\,\varphi},$

	$\displaystyle{\cal E}_{\phi,lm}$	$\displaystyle=$	$\displaystyle-i\,A^{TE}_{lm}\frac{h_{l}(x)}{l(l+1)}\frac{dP^{m}_{l}(\theta)}{d\theta}\,e^{i\,m\varphi}$
			$\displaystyle+i\,\frac{mA^{TM}_{l+1\,m}}{(l+1)(l+2)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l+1}(x))\right)\frac{P^{m}_{l+1}(\theta)}{\sin\theta}\,e^{im\varphi}$
			$\displaystyle+i\frac{m\,A^{TM}_{l-1m}}{l(l-1)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l-1}(x))\right)\frac{P^{m}_{l-1}(\theta)}{\sin\theta}\,e^{i\,m\varphi}.$

The electric field computed above, satisfies the boundary condition ${\cal E}_{\phi,lm}(R,\theta,\phi,t)$ =0. It is shown in section B that they also fit the boundary condition given by Eq.(14)

The magnetic counterpart ${\cal B}$ is given by (See Eq. 28, Eq. 5, Eq. 5)

$${\cal B}_{r,lm}(r,\theta,\phi,t)=A^{TE}_{lm}\frac{h_{l}(x)}{x}\,P^{m}_{l}(\theta)\,e^{i\,m\varphi},$$

	$\displaystyle{\cal B}_{\theta,lm}(r,\theta,\phi,t)$	$\displaystyle=$	$\displaystyle\frac{A^{TE}_{lm}}{l\,(l+1)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l}))\right)\frac{dP^{m}_{l}(\theta)}{d\theta}\,e^{i\,m\varphi}$
			$\displaystyle+A^{TM}_{l+1\,m}\frac{m}{(l+1)(l+2)}h_{l+1}(x)\frac{P^{m}_{l+1}(\theta)}{\sin\theta}\,e^{i\,m\varphi}$
			$\displaystyle+A^{TM}_{l-1\,m}\frac{m}{l\,(l-1)}h_{l-1}(x)\frac{P^{m}_{l-1}(\theta)}{\sin\theta}\,e^{i\,m\varphi},$

	$\displaystyle{\cal B}_{\phi,lm}(r,\theta,\phi,t)$	$\displaystyle=$	$\displaystyle i\frac{mA^{TE}_{lm}}{l\,(l+1)}\left(\frac{1}{x}\frac{d}{dx}(xh_{l}(x))\right)\frac{P^{m}_{l}(\cos\theta)}{\sin\theta}\,e^{im\varphi}$
			$\displaystyle+i\,A^{TM}_{l+1,m}\frac{h_{l+1}}{(l+1)(l+2)}\,\frac{d\,P^{m}_{l+1}(\cos\theta)}{d\,\theta}\,e^{i\,m\varphi}$
			$\displaystyle+iA^{TM}_{l-1,m}\frac{h_{l-1}(x)}{(l-1)\,l}\,\frac{d\,P^{m}_{l-1}(\cos\theta)}{d\theta}\,e^{i\,m\varphi}.$

For $l=1$, $m=1$ we obtain the result given in Deutsch (1955).

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.

We have

$$B_{r}=B_{0}\left[\left(\frac{R}{r}\right)^{3}\cos\theta+\frac{\alpha}{2}(3\cos\theta^{2}-1)\left(\frac{R}{r}\right)^{4}\right]$$

(36)

$$B_{\theta}=B_{0}\left[\frac{1}{2}\left(\frac{R}{r}^{3}\right)^{3}\sin\theta+\alpha\left(\frac{R}{r}\right)^{4}\cos\theta\sin\theta\right]$$

(37)

where $\alpha$ characterises the quadrupole component amplitude.

The radial electric field is

	$\displaystyle E_{r}$	$\displaystyle=$	$\displaystyle+\tilde{\Omega}B_{0}\left[\frac{1}{2}(3\cos\theta^{2}-1)\left(\frac{R}{r}\right)^{4}\right]+\frac{Q}{R^{2}}\left(\frac{R}{r}\right)^{2}$
		$\displaystyle+$	$\displaystyle\tilde{\Omega}B_{0}\left[+\frac{2}{5}\alpha\left((5\cos\theta^{3}-3\cos\theta\,)\left(\frac{R}{r}\right)^{5}-\cos\theta\left(\frac{R}{r}\right)^{3}\right)\right]$

where $\tilde{\Omega}=\Omega\,R/c$ 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 ($\tilde{\Omega}B_{0}<0$) 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 $\alpha>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.

Fig. 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) $A^{TE}_{1,0}=$-1, and $A^{TE}_{2,0}=$-2.5 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^{5}$ 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^{6}$ 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 remain 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_{0}$, 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 on Fig. 5 that the Goldreich-Julian density $n_{GJ}=\nabla\cdot\left(\vec{B}\times\vec{V}_{\Omega}\right)/4\pi\,c$ also has opposite signs at opposite poles ($V_{\Omega}$ 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.

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\sim en_{GJ}c$ where $n_{GJ}$ is the Goldreich-Julian density. Those primary electrons follow the lines of the magnetic field $\vec{B}$ radiate high-energy $\gamma$ rays via curvature radiation, and the gamma rays produce electron positron pair via the magnetic field $\vec{B}$, if $\vec{B}$ is strong enough, or by $\gamma$ 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.

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. Eq. (16, 21), and non-axially symmetric fields ($m\neq 0$) given by Eq. (6, 6).

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 section 7 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 section 8, 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.