Abstract

Magnetic fields inside and around neutron stars are at the heart of pulsar magnetospheric activity. Strong magnetic fields are responsible for quantum effects, an essential ingredient to produce leptonic pairs and the subsequent broadband radiation. The variety of electromagnetic field topologies could lead to the observed diversity of neutron star classes. Thus it is important to include multipolar components to a presumably dominant dipolar magnetic field. Exact analytical solutions for these multipoles in Newtonian gravity have been computed in recent literature. However, flat spacetime is not adequate to describe physics in the immediate surrounding of neutron stars. We generalize the multipole expressions to the strong gravity regime by using a slowly rotating metric approximation such as the one expected around neutron stars. Approximate formulas for the electromagnetic field including frame dragging are computed from which we estimate the Poynting flux and the braking index. Corrections to leading order in compactness and spin parameter are presented. As far as spindown luminosity is concerned, it is shown that frame dragging remains irrelevant. For high order multipoles starting from the quadrupole, the electric part can radiate more efficiently than the magnetic part. Both analytical and numerical tools are employed.

keywords:

Our understanding of pulsar physics has benefited from recent advances in multi-wavelength observational campaigns as well as from developments of new numerical tools able to investigate carefully its magnetosphere and the subsequent radiation mechanisms from a theoretical point of view. As the quality and quantity of data increases regularly, theoreticians are forced to improve the physics of their models in conjunction with the precision of their predictions. Detailed analysis of pulsar phase-resolved spectra and polarisation properties requires accurate models dealing with all possible perturbations departing from a too simple description of pulsar magnetospheres. General relativity belongs to one of these additional physical ingredient compulsory to investigate properly neutron star electrodynamics. This fact became clear in the past years. Indeed efficient pair creation in force-free magnetospheres seems to require frame dragging effects close to the magnetic poles (Philippov et al., 2015). Pétri (2016a) investigated in depth general-relativistic force-free magnetospheres. In the same vain Ruiz et al. (2014) look at the spindown luminosity and attempted to match neutron star exterior described in the force-free regime to its interior described by relativistic MHD, see also Paschalidis & Shapiro (2013). Konno & Kojima (2000) investigated the impact of general relativistic corrections to curvature radiation. Oscillations of neutron stars in general relativity were also of interest (Kojima & Hosonuma, 2000). Morozova et al. (2010) studied the influence of neutron star oscillations in general relativity on the plasma density in the magnetosphere for a aligned rotator. Curvature of space time and frame dragging effects on the surrounding electromagnetic field was already emphasized by Beskin (1990) and Muslimov & Tsygan (1992).

Any rotating field can be expanded into multipolar components. Thus Bonazzola et al. (2015) and Pétri (2015) showed how to compute exact analytical solutions to multipoles in closed form for flat space-time. Extension to general relativity is highly desirable. The simplest case is a rotating dipole for which Kojima et al. (2004) gave an approximate solution in general relativity comparing their results with previous analytical work of Rezzolla et al. (2001). An analytical estimate for the dipole spindown has been given by Rezzolla & Ahmedov (2004). Recently these authors extended their analysis of oscillations by adding damping due to heating and Joule dissipation (Rezzolla & Ahmedov, 2016). In this paper, we show how to extend those results to any multipole and to any order of accuracy.

Maxwell equations in presence of strong gravity using the 3+1 formalism are used to solve for arbitrary electromagnetic field configurations in vacuum as presented in Sec. 2. Exact solutions for static multipolar magnetic fields in Schwarzschild background metric in terms of hypergeometric functions are reminded in Sec. 3. In the same gravitational field, exact analytical solutions are found in terms of local confluent Heun functions as explained in Sec. 4. Frame dragging is included in an approximate fashion by numerically solving the system of elliptic equations (Helmholtz equations) for the unknown fields as exposed in Sec. 5. Our analytical treatment is supported by time-dependent numerical simulations of Maxwell equations in general relativity as presented in Sec. 6. Conclusions are drawn in Sec. 7.

In this section the general formalism to solve Maxwell equations in general-relativistic vacuum is reviewed. For the 3+1 split of spacetime we follow the conventions and definitions given by Pétri (2013).

2.1 Split of space-time metric

We split space-time into a 3+1 foliation such that the background metric gik is expressed as

ds2=gikdxidxk=α2c2dt2−γab(dxa+βacdt)(dxb+βbcdt)

(1)

where xi=(ct,xa), t is the time coordinate or universal time and xa some associated space coordinates. The metric signature is given by (+,−,−,−). α is the lapse function, βa the shift vector and γab the spatial metric of absolute space. By convention, latin letters from a to h are used for the components of vectors in absolute space, in the range {1,2,3}, whereas latin letters starting from i are used for four dimensional vectors and tensors, in the range {0,1,2,3}. A fiducial observer (FIDO) is defined by its 4-velocity ni such that

ni

=dxidτ=cα(1,−\mn@boldsymbolβ)

(2a)

ni

=(αc,{0}).

(2b)

This vector is orthogonal to the hyper-surface of constant time coordinate Σt. Its proper time τ is measured according to dτ=αdt. The relation between the determinants of the space-time metric g and the pure spatial metric γ is given by √−g=α√γ. For a slowly rotating neutron star, the lapse function is

α=√1−Rsr

(3)

and the shift vector

c\mn@boldsymbolβ=

−ωrsinϑeφ

(4a)

ω=

Rsacr3.

(4b)

We use a spherical coordinate system (r,ϑ,φ) and an orthonormal spatial basis (er,eϑ,eφ). The metric of a slowly rotating neutron star remains close to the usual flat space, except for the radial direction. Indeed the components of the spatial metric are given in Boyer-Lindquist coordinates by

γab=⎛⎜⎝α−2000r2000r2sin2ϑ⎞⎟⎠.

(5)

For this slow rotation approximation, the spatial metric does not depend on the spin frequency of the massive body but only on its mass M through α. This justifies the slow-rotation approximation. The spin parameter a is related to the angular momentum J by J=Mac. It follows that a has units of a length and should satisfy a≤Rs/2. Introducing the moment of inertia I, we also have J=IΩ. In the special case of a homogeneous and uniform neutron star interior with spherical symmetry, the moment of inertia reads

I=25MR2.

(6)

Thus the spin parameter can be expressed as

aRs=25RRsRrL.

(7)

For the remainder of this paper, we will use this expression to relate the spin parameter intervening in the metric to the spin frequency of the neutron star. From the above expression, note that the parameter a/Rs remains usually small.

2.2 Maxwell equations in general relativity

Maxwell equations in 3+1 notation take a traditional form close to the one known in flat space-time except that, as the reader should keep in mind, the three dimensional space is curved and differential operators defined according to the spatial metric γab. In vacuum, the system reads

∇⋅{B}

=0

(8a)

∇×{E}

=−1√γ∂t(√γ{B%
})

(8b)

∇⋅{D}

=0

(8c)

∇×{H}

=1√γ∂t(√γ{D}% ).

(8d)

The three dimensional vector fields are not independent, they are related by two important constitutive relations, namely

ε0{E}

=α{D}+ε0c\mn@boldsymbolβ×{B}

(9a)

μ0{H}

=α{B}−β×{D}ε0c

(9b)

ε0 is the vacuum permittivity and μ0 the vacuum permeability. The curvature of absolute space is taken into account by the lapse function factor α in the first term on the right-hand side and the frame dragging effect is included in the second term, the cross-product between the shift vector \mn@boldsymbolβ and the fields. Space curvature and frame dragging have therefore an imprint on the constitutive relations eq. (9a), (9b). From the auxiliary fields ({E},{H}) we get the Poynting flux through a sphere of radius r by computing the two dimensional integral on this sphere by

L=∫ω{E}∧{H}r2dω

(10)

where dω is the infinitesimal solid angle and ω the full sky angle of 4π sr.
This integral can be computed analytically in the asymptotic regime of large distances given in the wave zone. The Poynting flux should then be interpreted as the power radiated as seen by a distant observer and not the intrinsic spindown as measured on the neutron star surface. Indeed, due to gravitational redshift, energy is degraded when photons propagate from the surface to the observer and this affects also the measured power in the wave zone. Our spindown luminosity is computed according to the normalization done by this distant observer.

2.3 General solution to Maxwell equations

It is formally possible to give arbitrary solutions to Maxwell equations for divergencelessness electric and magnetic fields in vacuum in a stationary regime. Indeed their expansion into vector spherical harmonics reads, neglecting a possible monopolar ℓ=0 contribution,

D(r,ϑ,φ,t)=

∞∑ℓ=1ℓ∑m=−ℓ(∇×[fDℓ,m(r)Φℓ,m]+gDℓ,m(r)Φℓ,m)e−imΩt

(11a)

B(r,ϑ,φ,t)=

∞∑ℓ=1ℓ∑m=−ℓ(∇×[fBℓ,m(r)Φℓ,m]+gBℓ,m(r)Φℓ,m)e−imΩt.

(11b)

which correspond to stationary solutions expressed in the frame of a distant observer. The Φℓ,m are vector spherical harmonics defined and introduced in recent works by Pétri (2013). The functions gDℓ,m and gBℓ,m are related to the function fBℓ,m and fDℓ,m according to Maxwell equations by a linear scaling. Indeed, there exists a simple algebraic relation between gDℓ,m and fBℓ,m on one side, and between gBℓ,m and fDℓ,m on the other side such that

αgDℓ,m=

+iε0m~ωfBℓ,m

(11c)

αgBℓ,m=

−iμ0m~ωfDℓ,m

(11d)

where ~ω=Ω−ω is the rotation rate as measured by a local observer. After substitution in Maxwell equations, we get inhomogeneous Helmholtz equations for the potentials fDℓ,m and fBℓ,m. Indeed as shown by Pétri (2013), introducing the Helmholtz operator in curved spacetime by

Wℓ,m[f]≡α2[1r∂∂r(α2∂∂r(rf))−ℓ(ℓ+1)r2f]+m2~ω2c2f

(12)

the potentials fDℓ,m must satisfy

Wℓ[fDℓ,m]=3ε0α2ωr[fBℓ−1,m√(ℓ−1)(ℓ+1)Jℓ,m−fBℓ+1,m√ℓ(ℓ+2)Jℓ+1,m]

(13a)

where Jℓ,m=√ℓ2−m24ℓ2−1 and similarly for the magnetic field fBℓ,m

Wℓ[fBℓ,m]=−3μ0α2ωr[fDℓ−1,m√(ℓ−1)(ℓ+1)Jℓ,m−fDℓ+1,m√ℓ(ℓ+2)Jℓ+1,m].

(13b)

The boundary conditions on the neutron star surface are imposed on the electric field in the following way

α2∂r(rfDℓ,m)=ε0r~ω[√(ℓ+1)(ℓ−1)Jℓ,mfBℓ−1,m−√ℓ(ℓ+2)Jℓ+1,mfBℓ+1,m].

(14)

This last expression shows that the electric field strength is proportional to ~ω which contains the frame dragging effect. This leads to a lowering of the actual rotation rate of the star as seen by a local observer on the surface. Thus frame dragging decreases the electric field intensity and therefore also the spindown luminosity corrections due to rotation of spacetime. However, as shown later in this work, for realistic neutron star parameters, these corrections remain negligible. A second correction is induced by the space curvature and implies a additional factor α−2 compared to flat spacetime. The constants of integration are magnified but this effect is sometimes completely cancelled by the general relativistic spherical Hankel functions when ℓ=m as shown in the numerical results in Sec. 5. For other multipoles with ℓ>m, compensation is only partial. These conclusions are discussed in detail in the numerical approximate solution Sec. 5.

The Helmholtz operator is conveniently rewritten by introducing the tortoise coordinate r∗ such that

r∗=r+Rsln(rRs−1).

(15)

We are looking for solutions describing outgoing waves that reduce to eikmr in flat space time, thus we introduce another unknown field uB/Dℓ,m such that

fB/Dℓ,m(r)=uB/Dℓ,m(r)eikmr∗r

(16)

where km=m/rL. The curved spacetime Helmholtz operator in terms of these new dependent variables uℓ,m becomes

For one single multipole with potential specified at the surface by fBℓ,m(R), the boundary conditions for the electric field reduces to

α2∂r(uDℓ−1,meikmr∗)

=−ε0r~ω√(ℓ−1)(ℓ+1)Jℓ,mfBℓ,m

(19a)

α2∂r(uDℓ+1,meikmr∗)

=+ε0r~ω√ℓ(ℓ+2)Jℓ+1,mfBℓ,m

(19b)

to be evaluated at the surface for r=R.

2.4 Wave zone and Poynting flux

The Poynting flux of a rotating multipole can be most easily computed in the asymptotic flat spacetime at very large distance. Because of energy conservation law, the flux leaving the star must reach infinity, there is no absorption layer in between. In the wave zone, the expressions (11) can be drastically reduced by the fact that the potential functions behave asymptotically as fℓ,m(r)≈u∞ℓ,meikmr/r where limr→+∞uℓ,m(r)=u∞ℓ,m. Neglecting the axisymmetric mode decreasing much faster, like r−(ℓ+1), the electromagnetic field becomes in the limit of large distances r≫rL

Bw

=∑ℓ≥1,m≠0−ieim(kr−Ωt)rkm(uB,∞ℓ,mΨℓ,m+μ0cuD,∞ℓ,mΦℓ,m)

(20a)

Dw

=∑ℓ≥1,m≠0−ieim(kr−Ωt)rkm(uD,∞ℓ,mΨℓ,m−ε0cuB,∞ℓ,mΦℓ,m)

(20b)

=

ε0cBw∧n.

(20c)

Equation (20c) shows that the solution behaves as a monochromatic plane wave propagating in the radial direction n=er at frequency Ω. The time averaged Poynting flux is therefore

S=Dw∧B∗w2μ0ε0

(21)

where B∗w is the complex conjugate of Bw. Integrating the radial component of the Poynting vector along the solid angle ω we get the power radiated, using the orthonormality of the vector spherical harmonics, such that

L=∫ωS^rr2dω=c2μ0∑ℓ≥1,m≠0k2m(|uB,∞ℓ,m|2+μ20c2|uD,∞ℓ,m|2).

(22)

The spin down luminosity L is independent of the radius as it should from the energy conservation law. Equation (22) represents the most general expression for the magneto-multipole losses from an arbitrary multipole magnetic field in general relativity. As soon as the constants of integration uD/B,∞ℓ,m are known, the full stationary electromagnetic field is determined and its subsequent properties such as spindown and magnetic topology. Our main goal is to fix these constants either with some analytical argument or more accurately via numerical integration of the elliptic problems related to the Helmholtz equations in curved spacetime. We next recall the exact analytical solutions of the static multipoles in general relativity and then look for the stationary rotating multipoles solutions found numerically.

Finding explicit exact solutions to the rotating multipole problem is difficult because there is no known analytical solution to Helmholtz equations in general relativity in the Schwarzschild metric. Nevertheless, in the static limit of non rotating neutron stars, it is possible to find exact close formula for the multipoles to any order.

From the expansion into vector spherical harmonics eqs. (11), each function fD/Bℓ,m has to satisfy a homogeneous second order linear differential equation in Schwarzschild space-time such that

α2r∂r(α2∂r(rfℓ,m))−α2ℓ(ℓ+1)r2fℓ,m=0.

(23)

Introducing the new unknown function ϕℓ,m=rfℓ,m we get the simple differential equation

∂r(α2∂r(ϕℓ,m))−ℓ(ℓ+1)r2ϕℓ,m=0

(24)

to be solved with appropriate boundary conditions, namely vanishing potentials at spatial infinity. Moreover, introducing the normalized inverse radial coordinate by x=Rs/r, the functions ϕℓ,m will be solution of

x2(1−x)ϕ′′ℓ,m+x(2−3x)ϕ′ℓ,m−ℓ(ℓ+1)ϕℓ,m=0

(25)

which reduces to the hypergeometric differential equation by the change of unknown function ϕℓ,m=xℓvℓ,m to

x(1−x)v′′ℓ,m+(2(ℓ+1)−(2ℓ+3)x)v′ℓ,m−ℓ(ℓ+2)vℓ,m=0.

(26)

Setting the parameters a=ℓ,b=ℓ+2,c=2(ℓ+1) we indeed find the standard form of the hypergeometric differential equation (Olver, 2010). The only solution vanishing at infinity is

ϕℓ,m(x)=Cxℓ2F1(ℓ,ℓ+2,2(ℓ+1),x)

(27)

which is the expression found by Muslimov & Tsygan (1986). C is a constant of integration imposed by the boundary condition on the neutron star surface. The solution does not depend on the quantum number m in the static limit. The constant C is chosen such that the asymptotic value of the potentials converge to their flat spacetime counterpart at large distance. This represents our normalization of the magnetic field strength throughout the paper. We now give the exact analytical expression for the first four multipoles ℓ∈{1,2,3,4}.

3.1 The magnetic dipole ℓ=1

We start our discussion with the general-relativistic magnetic dipole. Introducing the vector spherical harmonics expansion, the static aligned dipole frozen into the neutron star is conveniently written as

fB1,0=2√6πBR3R2s[ln(1−x)x+1+x2]≈−2√2π3BR3r2[1+34Rsr]

(28)

(recall that x=Rs/r). This solution already found by Ginzburg & Ozernoy (1964) asymptotes to the flat space-time field at very large distances r≫Rs. The aligned dipolar magnetic field components are given in an orthonormal basis by

B^r

=−3BR3R3scosϑL1(x)≈2BR3r3cosϑ[1+34Rsr]

(29a)

B^ϑ

=3BR3R3ssinϑT1(x)≈BR3r3sinϑ[1+Rsr]

(29b)

B^φ

=0

(29c)

where we introduced the longitudinal and transversal part by

L1(x)

=x(x+2)+2ln(1−x)≈−23x3[1+34x+o(x2)]

(30a)

T1(x)

=(2−x)x+2(1−x)ln(1−x)√1−x≈13x3[1+x+o(x2)].

(30b)

The static perpendicular dipole frozen into the neutron star is conveniently written with the normalization

fB1,1=−√2fB1,0

(31)

meaning inclining the previous aligned dipole to 90\degr with respect to the rotation axis leading to the magnetic field components

B^r

=−3BR3R3seiφsinϑL1(x)≈2BR3r3eiφsinϑ[1+34Rsr]

(32a)

B^ϑ

=−3BR3R3seiφcosϑT1(x)≈−BR3r3eiφcosϑ[1+Rsr]

(32b)

B^φ

=−3BR3R3sieiφT1(x)≈−BR3r3ieiφ[1+Rsr].

(32c)

Gravitational corrections to first order expressed by the coefficient Rs/r are shown to better grab the increase in magnetic field components. The amplification is different depending on the component under consideration thus the field line topology is also modified with respect to a flat spacetime dipole.

3.2 The magnetic quadrupole ℓ=2

Let us perform the same expansion to the general-relativistic magnetic quadrupole such that it can be expressed inside the star by

fB2,0

=−83√10π3BR4R3sx(x(x+6)−24)+6(3x−4)log(1−x)x2

(33a)

≈

−4√2π15BR4r3[1+43Rsr].

(33b)

The quadrupolar magnetic field components for the axisymmetric mode m=0 are given in an orthonormal basis by

B^r

=−103BR4R4s(3cos2ϑ+1)L2(x)≈BR4r4(3cos2ϑ+1)[1+43Rsr]

(34a)

B^ϑ

=20BR4R4ssin2ϑT2(x)≈2BR4r4sin2ϑ[1+32Rsr]

(34b)

B^φ

=0

(34c)

where we introduced the longitudinal and transversal part by

L2(x)

=−x(x(x+6)−24)+6(3x−4)log(1−x)x

(35a)

=−310x4[1+43x+o(x2)]

(35b)

T2(x)

=x((x−12)x+12)+6(x−2)(x−1)log(1−x)x√1−x

(35c)

=110x4[1+32x+o(x2)].

(35d)

The static (2,1) quadrupole frozen into the neutron star is conveniently written with the normalization

fB2,1=−12√32fB2,0

(36)

giving the m=1 mode by

B^r

=−5BR4R4seiφsin2ϑL2(x)≈32BR4r4eiφsin2ϑ[1+43Rsr]

(37a)

B^ϑ

=−10BR4R4seiφcos2ϑT2(x)≈−BR4r4eiφcos2ϑ[1+32Rsr]

(37b)

B^φ

=−10BR4R4sieiφcosϑT2(x)≈−BR4r4ieiφcosϑ[1+32Rsr].

(37c)

The static (2,2)-quadrupole frozen into the neutron star is conveniently written with the normalization

fB2,2=√32fB2,0

(38)

leading to the m=2 mode

B^r

=−10BR4R4se2iφsin2ϑL2(x)≈3BR4r4e2iφsin2ϑ[1+43Rsr]

(39a)

B^ϑ

=−10BR4R4se2iφsin2ϑT2(x)≈−BR4r4e2iφsin2ϑ[1+32Rsr]

(39b)

B^φ

=−20BR4R4sie2iφsinϑT2(x)≈−2BR4R4sie2iφsinϑ[1+32Rsr].

(39c)

3.3 The magnetic hexapole ℓ=3

Next the magnetic hexapole existing inside the neutron star is written as

fB3,0

=52√7π3BR5R4s12(6x2−20x+15)log(1−x)+x(x(x(x+12)−150)+180)x3

(40a)

≈−2√π21BR5r4[1+158Rsr].

(40b)

The hexapolar magnetic field components for the axisymmetric mode m=0 are given in an orthonormal basis by