GR multipolar electromagnetic fields

# Multipolar electromagnetic fields around neutron stars: general-relativistic vacuum solutions

## 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:
magnetic fields - methods: analytical - methods: numerical - stars: neutron - stars: rotation - pulsars: general
12

## 1 Introduction

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.

## 2 Maxwell equations and their solutions

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  is expressed as

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

where , is the time coordinate or universal time and some associated space coordinates. The metric signature is given by . is the lapse function, the shift vector and the spatial metric of absolute space. By convention, latin letters from to are used for the components of vectors in absolute space, in the range , whereas latin letters starting from are used for four dimensional vectors and tensors, in the range . A fiducial observer (FIDO) is defined by its 4-velocity  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 . Its proper time  is measured according to . The relation between the determinants of the space-time metric  and the pure spatial metric  is given by . 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  and an orthonormal spatial basis . 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  through . This justifies the slow-rotation approximation. The spin parameter  is related to the angular momentum  by . It follows that has units of a length and should satisfy . Introducing the moment of inertia , we also have . 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  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 . 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)

is the vacuum permittivity and 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  and the fields. Space curvature and frame dragging have therefore an imprint on the constitutive relations eq. (9a), (9b). From the auxiliary fields we get the Poynting flux through a sphere of radius by computing the two dimensional integral on this sphere by

 L=∫ω{E}∧{H}r2dω (10)

where is the infinitesimal solid angle and the full sky angle of  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 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 and . 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  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 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 as shown in the numerical results in Sec. 5. For other multipoles with , 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 such that

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

We are looking for solutions describing outgoing waves that reduce to in flat space time, thus we introduce another unknown field  such that

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

where . The curved spacetime Helmholtz operator in terms of these new dependent variables  becomes

 Wℓ,m[fℓ,m]={1r3(1−Rsr)[r(r−Rs)u′′ℓ,m+(2ikmr2+Rs)u′ℓ,m−ℓ(ℓ+1)uℓ,m]+k2m[(1−ωΩ)2−1]uℓ,mr}eikmr∗. (17)

In terms of , the inhomogeneous Helmholtz equations become for the electric field

 1r[r(r−Rs)uDℓ,m′′+(2ikmr2+Rs)uDℓ,m′−ℓ(ℓ+1)uDℓ,m]+k2m[(1−ωΩ)2−1]ruDℓ,mα2=3ε0ω[uBℓ−1,m√(ℓ−1)(ℓ+1)Jℓ,m−uBℓ+1,m√ℓ(ℓ+2)Jℓ+1,m] (18a) and similarly for the magnetic field 1r[r(r−Rs)uBℓ,m′′+(2ikmr2+Rs)uBℓ,m′−ℓ(ℓ+1)uBℓ,m]+k2m[(1−ωΩ)2−1]ruBℓ,mα2=−3μ0ω[uDℓ−1,m√(ℓ−1)(ℓ+1)Jℓ,m−uDℓ+1,m√ℓ(ℓ+2)Jℓ+1,m]. (18b)

For one single multipole with potential specified at the surface by , 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 .

### 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 where . Neglecting the axisymmetric mode decreasing much faster, like , the electromagnetic field becomes in the limit of large distances

 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 at frequency . The time averaged Poynting flux is therefore

 S=Dw∧B∗w2μ0ε0 (21)

where is the complex conjugate of . 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  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  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.

## 3 Exact static multipole solutions

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 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 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 , the functions 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 to

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

Setting the parameters 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). is a constant of integration imposed by the boundary condition on the neutron star surface. The solution does not depend on the quantum number  in the static limit. The constant  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 .

### 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 ). This solution already found by Ginzburg & Ozernoy (1964) asymptotes to the flat space-time field at very large distances . 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 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 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 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 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)

 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 are given in an orthonormal basis by

 B^r =−35BR54R5scosϑ(5cos2ϑ−3)L3(x)≈BR5r5cosϑ(5cos2ϑ−3)[1+158Rsr] (41a) B^ϑ =105BR516R5s(sinϑ+5sin3ϑ)T3(x)≈316BR5r5cosϑ(5cos2ϑ−3)[1+2Rsr] (41b) B^φ =0 (41c)

where we introduced the longitudinal and transversal part by

 L3(x) =(12(6x2−20x+15)log(1−x)+x(x(x(x+12)−150)+180))x2 (42a) =−435x5[1+158x+o(x2)] (42b) T3(x) =(2−x)x((x−30)x+30)+12(1−x)((x−5)x+5)log(1−x)x2√1−x (42c) =135x5[1+2x+o(x2)]. (42d)

The static (3,1)-quadrupole frozen is normalized according to

 fB3,1=−32√3fB3,0 (43)

such that the mode becomes

 B^r =−35BR5R5seiφ(sinϑ+5sin3ϑ)L3(x)≈4BR5r5eiφ(sinϑ+5sin3ϑ)[1+158Rsr] (44a) B^ϑ =−35BR5R5seiφ(cosϑ+15cos3ϑ)T3(x)≈−BR5r5eiφ(cosϑ+15cos3ϑ)[1+2Rsr] (44b) B^φ =−70BR5R5sieiφ(5cos2ϑ+3)T3(x)≈−2BR5r5ieiφ(5cos2ϑ+3)[1+2Rsr]. (44c)

The static (3,2)-quadrupole frozen is normalized according to

 fB3,2=√215fB3,0 (45)

such that the mode gives

 B^r =−35BR54R5se2iφsin2ϑcosϑL3(x)≈BR5r5e2iφsin2ϑcosϑ[1+158Rsr] (46a) B^ϑ =35BR516R5se2iφ(sinϑ−3sin3ϑ)T3(x)≈BR516r5e2iφ(sinϑ−3sin3ϑ)[1+2Rsr] (46b) B^φ =−35BR54R5sie2iφsin2ϑT3(x)≈−BR54r5ie2iφsin2ϑ[1+2Rsr]. (46c)

The static (3,3)-quadrupole is normalized to

 fB3,3=8√5fB3,0 (47)

thus we get for the mode

 B^r =35BR5R5se3iφsin3ϑL3(x)≈−4BR5r5e3iφsin3ϑ[1+158Rsr] (48a) B^ϑ =105BR5R5se3iφsin2ϑcosϑT3(x)≈3BR5r5e3iφsin2ϑcosϑ[1+2Rsr] (48b) B^φ =105BR5R5sie3iφsin2ϑT3(x)≈3BR5r5ie3iφ