Impact of an off-centred dipole on neutron star binaries

TL;DR

This paper investigates how an off-centred magnetic dipole affects neutron star electromagnetic properties, using numerical simulations and analytical comparisons, with implications for binary neutron star eccentricities.

Contribution

It provides a detailed analysis of off-centred dipole magnetic fields and their effects, combining simulations with analytical models, which was not thoroughly explored before.

Findings

Good agreement between simulations and analytical models for large dipole displacements.

Off-centred dipoles influence neutron star spin-down and electromagnetic kicks.

Proposes an alternative explanation for binary neutron star eccentricity distributions.

Abstract

Neutron stars are strongly magnetized rotating compact objects. Therefore they also produce huge electric fields accelerating particles to ultra-relativistic energies. The simplest magnetic topology is a dipole traditionally located at the stellar centre. In this paper, we reinvestigate the consequences of an off-centred rotating magnetic dipole, showing accurate magnetic field line geometries, the associated spin-down luminosity as well as the corresponding electromagnetic kick and torque imprinted to the neutron star. Results are obtained by time dependent numerical simulations of Maxwell equations in vacuum using pseudo-spectral methods. We compare our results to known analytical expressions available to lowest order in the parameter $\epsilon = d/R$, where $d$ is the displacement of the dipole from the stellar centre and $R$ the neutron star radius. We found good agreement between our numerical computations and our analytical approximations even for well off-centred dipoles having large displacements with a sizeable fraction of the radius, i.e. $\epsilon \lesssim 1$. An explanation for binary neutron star eccentricity distribution functions is given with an emphasize on highly eccentric systems as an alternative scenario to traditional binary formation.

Equations

$$\mathbf{d}=d(\sin \delta \cos (\Omega t){\mathbf{e}}_{\mathrm{x}}+\sin \delta \sin (\Omega t){\mathbf{e}}_{\mathrm{y}}+\cos \delta {\mathbf{e}}_{\mathrm{z}})$$

$$\mathbf{m}=\sin \alpha \cos (\beta +\Omega t){\mathbf{e}}_{\mathrm{x}}+\sin \alpha \sin (\beta +\Omega t){\mathbf{e}}_{\mathrm{y}}+\cos \alpha {\mathbf{e}}_{\mathrm{z}}.$$

$$\mathbf{B}=\frac{B{R}^3}{\Vert \mathbf{r}-\mathbf{d}{\Vert }^3}\left[\frac{3\text{mn@boldsymbol}\mu \cdot (\mathbf{r}-\mathbf{d})}{\Vert \mathbf{r}-\mathbf{d}{\Vert }^2}(\mathbf{r}-\mathbf{d})-\text{mn@boldsymbol}\mu \right]$$

$$L_{m=1}$$

$$L_{m=2}$$

$$L_{\mathrm{dip}}=\frac{8\pi }{3{\mu }_0{c}^3}{\Omega }^4{B}^2{R}^6.$$

$$F_{m=1}$$

$$F_{m=2}$$

$$\mathbf{K}=R^3\iint [{\sigma }_{\mathrm{s}}\mathbf{n}\wedge \mathbf{E}+(\mathbf{B}\cdot \mathbf{n}){\mathbf{i}}_{\mathrm{s}}]d\Omega$$

$$K_{\mathrm{x}}^{\mathrm{E}}$$

$$3(a^4(3744{ϵ}^2+875)+a^2(1750-21870{ϵ}^2)+100(81{ϵ}^2+140))\sin \beta )$$

$$K_{\mathrm{y}}^{\mathrm{E}}$$

$$+200(27{ϵ}^2+70))\cos \beta +350a^3(a^2(243{ϵ}^2+5)+243{ϵ}^2)\sin \beta )$$

$$K_{\mathrm{x}}^{\mathrm{B}}$$

$$+(3a^4(75228{ϵ}^2-11375)$$

$$-1470a^2(66{ϵ}^2-25)-100(324{ϵ}^2+245))\sin \beta )$$

$$K_{\mathrm{y}}^{\mathrm{B}}$$

$$-100(216{ϵ}^2+245))\cos \beta$$

$$-700a(a^2(207{ϵ}^2-50)+50)\sin \beta )$$

$$K_{\mathrm{z}}^{\mathrm{B}}$$

$${\ddot{\mathbf{r}}}_1$$

$${\ddot{\mathbf{r}}}_2$$

$$\mathbf{R}=\frac{m_1{\mathbf{r}}_1+m_2{\mathbf{r}}_2}{m_1+m_2}$$

$$\ddot{\mathbf{R}}=\frac{{\mathbf{F}}_1+{\mathbf{F}}_2}{m_1+m_2}.$$

$${\text{mn@boldsymbol}\rho }_1$$

$${\text{mn@boldsymbol}\rho }_2$$

$${\ddot{\text{mn@boldsymbol}\rho }}_1$$

$${\ddot{\text{mn@boldsymbol}\rho }}_2$$

$$\ddot{\mathbf{r}}=-G(m_1+m_2)\frac{\mathbf{r}}{r^3}+\frac{{\mathbf{F}}_2}{m_2}-\frac{{\mathbf{F}}_1}{m_1}.$$

$$\mathbf{A}=\frac{{\mathbf{F}}_2}{m_2}-\frac{{\mathbf{F}}_1}{m_1}$$

$$e(t)=\left|\sin i_0\sin \left(\frac{3A}{2\Omega a}t\right)\right|$$

$$T_e=\frac{\pi \Omega a}{3A}=\frac{\pi }{3A}\sqrt{\frac{G(m_1+m_2)}{a}}.$$

$$L(t)=\frac{L_0}{(1+t/{\tau }_{\mathrm{c}}{)}^2}$$

$$e(t)=\left|\sin i_0\sin \left(\frac{3}{2\Omega a}{\int }_0^{t}A(t)dt\right)\right|.$$

$${\int }_0^{{\tau }_{\mathrm{c}}}A(t)dt=\frac{A(t=0){\tau }_{\mathrm{c}}}{2}.$$

$$e({\tau }_{\mathrm{c}})=\left|\sin i_0\sin \left(\frac{3A{\tau }_{\mathrm{c}}}{4\Omega a}\right)\right|.$$

$$e({\tau }_{\mathrm{c}})=\left|\sin i_0\sin \left(\frac{9{\pi }^{5/3}}{5\times 2^{1/3}}ϵ\frac{I{P}^{-2}{P}_{\mathrm{orb}}^{1/3}}{m_{1}c\sqrt{G(m_1+m_2)}}\right)\right|$$

$$e({\tau }_{\mathrm{c}})$$

Full text

††pubyear: 2019††pagerange: Impact of an off-centred dipole on neutron star binaries–Impact of an off-centred dipole on neutron star binaries

1 Introduction

Magnetic field and rotation play an important role in the life and evolution of stars. This is especially true for magnetized compact objects like white dwarfs and neutron stars. There the magnetic field is often assumed to be dipolar, which corresponds to the lowest order expansion of any magnetic field geometry not possessing magnetic monopoles. Moreover, for spherical stars, the dipolar magnetic moment is located exactly at the centre of the star. The main reason for this assumption is simplicity. However, such approximation is hardly supported by any physical explanation. Actually the probability to find such a coincidence seems rather weak. It is conceivable that currents in the stellar interior lead to field structures not centred and also not strictly dipolar. For strongly magnetized stars such as neutron stars, this asymmetry in the topology could have observable consequences. Harrison & Tademaru (1975) already mentioned the possibility to explain the high kick velocity above 100 km/s induced by the electromagnetic force imprinted to the newly born star. Some refinements were soon brought by Tademaru (1976) who added polarization effects and thus contributions from an additional electric field component.

In the case of a centred dipole, the torque applied on a neutron star in vacuum always predicts an alignment as explained by Davis & Goldstein (1970) and Goldreich (1970). This fact is also discussed by Philippov et al. (2014). Moreover, this electromagnetic torque has been computed following several expressions, either using the electromagnetic stress-energy tensor or the Lorentz force acting inside and on the surface of the neutron star. Both approaches agree on the regular electromagnetic torque but not on the anomalous torque. A single analytical formula has been found for this anomalous torque but the scaling factor depends on the expression used to compute the torque, from the stress-energy tensor or directly from the Lorentz force. Also the electric field contribution, sometimes omitted, should be included, see for instance Beskin & Zheltoukhov (2014) for a critical review of these discrepancies. These complications did not prevent Good & Ng (1985) from considering contributions from a quadrupole component to the total torque.

Decentred dipoles furnish simple explanations for the observations of planetary magnetic fields. Indeed, Komesaroff (1976) already recognized that a decentred dipole better fits the Pioneer data about Jupiter than a centred dipole. It was used to explain the asymmetries in the radio emission by minimizing the quadrupolar term. Landstreet (1980) presented a review about magnetic fields in non degenerated stars and concluded that a decentred oblique dipole can reasonably fit the magnetic topology. Hints for off-centred dipoles or dipole plus quadrupole fields are given by Putney & Jordan (1995) from polarization observations of white dwarfs. In our solar system some interesting results have been found for the outer planets. Indeed, inclination angles and displacements have been constrained for instance for Uranus with an inclination of $60\degr$ and an off-centring of $d=0.3\,R_{\rm U}$ (Ness et al., 1986) and Neptune with an inclination of $47\degr$ and an off-centring of $d=0.55\,R_{\rm N}$ (Ness et al., 1989) where $R_{\rm U}$ and $R_{\rm N}$ are Uranus and Neptune radius respectively. However, distinguishing an off-centred dipole from a combination of multipoles seems difficult on the base only of observations (Martin & Wickramasinghe, 1984).

Deviation from a pure and centred rotating dipole is also an attractive assumption to get more insight into neutron star and pulsar physics. Indeed, a simple prescription for a distorted dipolar field was proposed by Harding & Muslimov (2011) to enhance the pair production rate above polar caps. Non dipolar fields have also been suggested to explain the anomalous braking index of pulsars as discussed by Barsukov & Tsygan (2010). With the wealth of new and very accurate multi-wavelength observations of pulsed emission in pulsars including polarization, time is ripe to go further than the almost exclusively used centred dipole to sharpen our view on pulsar magnetospheres. An off-centred rotating dipole anchored in a perfectly conducting sphere with finite radius leads very naturally to multipolar components to any order with weighting coefficients solely related to the geometry of the dipole. We use this assumption as a starting point for the justification of multipolar components.

In this paper we compute accurate numerical solutions for the electromagnetic field in vacuum outside an off-centred rotating dipole, valid for any displacement $d$ such that $\epsilon=d/R<1$, where $d$ is the distance from the centre of the neutron star and $R$ its radius. In Sec. 2, we recall the basic model of an off-centred rotating dipole. Then some examples of field lines are presented in Sec. 3. Next, in Sec. 4 we compute the spin-down luminosity expected from such a system and compare it with previous works. The same comparison is done for the associated electromagnetic kick which is estimated in Sec. 5 and the electromagnetic torque in Sec. 6. Implications for binaries containing neutron stars are highlighted in Sec. 7. Conclusions are drawn in Sec. 8.

2 Rotating off-centred dipole

We remind the geometrical set-up for an off-centred dipole, following the notations introduced by Pétri (2016). The neutron star is a perfectly conducting sphere of radius $R$ in solid body rotation at an angular rate $\Omega$. Its magnetic moment $\mn@boldsymbol{\mu}$ is located inside this sphere at a point $M$ such that at any time $t$ its position vector is decomposed into

[TABLE]

where $d$ is the distance from the stellar centre and $\delta$ the colatitude and $(\mathbf{e}_{\rm x},\mathbf{e}_{\rm y},\mathbf{e}_{\rm z})$ is a cartesian orthonormal basis, Fig. 1. The relevant displacement is given by the normalised parameter $\epsilon=d/R<1$. Entrainment by the star is included in the phase term $\Omega\,t$. Meanwhile, the magnetic moment $\mn@boldsymbol{\mu}$ points toward a direction depicted by the angles $(\alpha,\beta)$ and given by the unit vector $\mn@boldsymbol{\mu}=\mu\,\mathbf{m}$ such that

[TABLE]

All important parameters are summarized in Fig. 1.

Inside the star, the magnetic field is given by a static dipole such that

[TABLE]

where $B$ is the surface magnetic field at the equator and $\mathbf{r}$ the position vector.

Starting from this static solution, we perform time dependent numerical simulations solving Maxwell equations in vacuum. We use the pseudo-spectral code developed and discussed in detail in Pétri (2014), applying it to flat spacetime. Boundary conditions at the stellar surface are given by the continuity of the radial component of the magnetic field $B_{r}$ and the tangential component of the electric field where $\mathbf{E}^{\prime}=\mathbf{0}$ in the corotating frame inside the star. At large distances, we enforce outgoing wave boundary conditions.

We run several sets of simulations to scan a full range of geometries, varying the angles $\alpha$, $\beta$ and $\delta$ and the displacement $d$. We next summarized our results by first showing some magnetic field line structures, then compute the spin down luminosities and eventually the associated electromagnetic kick and torque.

3 Field lines

As a first result of a rotating off-centred dipole, we plot magnetic field lines in the equatorial plane for a perpendicular rotator with $\alpha=\delta=90\degr$ for several values of the displacement $d$ and angle $\beta$. For small displacements $d\ll R$ the topology is very similar to the centred dipole given by Deutsch solution. In all cases, a two armed spiral develops, rotating at a constant speed equal to the rotation period of the star $\Omega$.

An example of field lines is shown for $\beta=0\degr$ and $\epsilon=0.3$ in Fig. 2 in red solid line and compared to the centred dipole given by Deutsch solution in blue dashed line. At large distances, in the wave zone we observed a shift in phase for the position of the two-armed spiral with respect to the centred case. This shift is proportional to $R/r_{\rm L}$, therefore too weak to be detectable for realistic pulsar parameters. Again, an off-centred dipole can only reasonably be detected when examining electromagnetic activity and radiation occurring close to the surface.

Another example is shown for $\beta=90\degr$ and $\epsilon=0.3$ in Fig. 3 in red solid line and compared to the centred dipole given by Deutsch solution in blue dashed line. Here again a two-armed spiral forms, reminiscent of the centred dipole. However, in this second case, at large distance, well outside the light-cylinder, both structure are very similar, only the dipolar part propagates significantly into the wave zone. It is hard to detect an off-centred dipole by inspecting the large scale field structure because the typical two-armed spiral structures overlap. The asymmetry is only clearly visible close to the neutron star surface.

Next we diagnose quantitatively the effect of an off-centering by computing relevant physical parameters that are the spindown luminosity, the electromagnetic kick and its associated torque.

4 Spindown luminosities

The spindown luminosity is an important characteristic of any rotating magnetic multipole. It determines the secular evolution of the rotational period. In Pétri (2016), we showed that a finite size off-centred dipole looses angular momentum mainly through its $m=1$ (dipole and quadrupole) and $m=2$ (quadrupole) modes such that each contribution splits into

[TABLE]

where $a=R/r_{\rm L}$ and the centred perpendicular dipole spindown is

[TABLE]

For brevity, we only showed expressions (4) valid for $\delta=90\degr$ although it is possible to give (lengthy) expressions for any angle $\delta$. Note that the quadrupolar contribution to the spindown is of the same order of magnitude as the perturbation in the dipole spindown. Both terms scale as $a^{2}\,\epsilon^{2}$ but show a different dependence with respect to the inclination angle $\alpha$, the $m=1$ mode adds a $\cos^{2}\alpha$ contribution whereas the $m=2$ mode adds a $\sin^{2}\alpha$ contribution. However, both corrections are independent of the angle $\beta$.

We compare these analytical approximations to the results of our numerical simulations where we also found that the spindown does not depend on the angle $\beta$ as expected from eq. (4). We performed a set of runs with relevant geometric parameters by varying the set $(\alpha,\beta,\delta,\epsilon)$ and choosing different rotation periods symbolized by the adimensionalized parameter $a$. For $a=0.2$ we summarize the simulation outputs for $\alpha=\{0\degr,30\degr,60\degr,90\degr\}$ in Fig. 4. The luminosity computed from the simulation $L(\epsilon)$ is compared to the exact vacuum centred dipole $L_{\rm exact}=L_{\rm dip}\,\sin^{2}\alpha$. The corresponding analytical expectations are shown by solid color curves. The agreement is good, the $\epsilon^{2}$ dependence is clearly retrieved and follows the formula given in eq. (4) even for a large off-centred distance with $\epsilon=0.5$. For such a high displacement the quadrupolar expansion is not sufficient to recover the field line topology close to the star but it suffices to compute accurately the spindown luminosity because higher order contributions to the spindown luminosities scale at least as $a^{3}\,\epsilon^{3}$. Therefore our approximate expression detailed in eq. (4) remains valid even for large displacements $\epsilon\lesssim 1$ as long as the spindown is concerned.

5 Electromagnetic kick

In Pétri (2016) we gave approximate expressions for the electromagnetic kick felt by the neutron star when radiating as an off-centred dipole. The dipolar $m=1$ and quadrupolar $m=2$ contributions are given respectively by

[TABLE]

These formula are again only valid for $\delta=90\degr$. The quadrupolar term remains negligible and a factor $a^{2}\,\epsilon^{2}$ weaker than the dipolar term. It can be dropped from the estimate without introducing much error in the estimate of the electromagnetic force. Indeed, we computed this force from our set of runs presented in the previous section. A compilation of the results is shown for $\alpha=30\degr$ in Fig. 5 and for $\alpha=60\degr$ in Fig. 6. The linear scaling with respect to $\epsilon$ is retrieved to good accuracy as seen from the color solid lines representing the analytical expectations.

Actually, both plots are very similar because they only differ by the expression $\cos\alpha\,\sin\alpha$ which are numerically identical if $\alpha=30\degr$ or $\alpha=60\degr$.

Next we show the dependence on the angle $\beta$ for $\alpha=30\degr$ in Fig.7 and for $\alpha=60\degr$ in Fig.8. The accuracy with only the dipolar term is already remarkable, pointing out the $\sin\beta$ dependence. Only the $\beta=0\degr$ geometry is not well reproduced. Again, both kicks are hardly distinguishable because of the $\cos\alpha\,\sin\alpha$ dependence and the particular value of $\alpha$ chosen.

We finish this theoretical study by the electromagnetic torque calculation.

6 Electromagnetic torque

Following Beskin & Zheltoukhov (2014), we compute the electromagnetic torque exerted on the neutron star surface according to the Laplace force given by

[TABLE]

where $\mathbf{n}$ is the unit normal to the surface, $\sigma_{\rm s}=\varepsilon_{0}\,[\mathbf{E}]\cdot\mathbf{n}$ represents the surface charge density and $\mu_{0}\,\mathbf{i}_{\rm s}=\mathbf{n}\wedge[\mathbf{B}]$ the surface current density. The notation $[\mathbf{F}]$ means the jump of the vector field $\mathbf{F}$ across the layer. Because of the perfect conductor assumption inside the star, there is no volume contribution to the torque as by definition $\rho\,\mathbf{E}+\mathbf{j}\wedge\mathbf{B}=\mathbf{0}$ where $(\mathbf{E},\mathbf{B})$ is the electromagnetic field, $\rho$ the charge density and $\mathbf{j}$ the current density inside the star. $\varepsilon_{0}$ and $\mu_{0}$ are the vacuum permittivity and permeability.

Expressions for the torque are lengthy for a decentred dipole. We only give the approximation for $\delta=90\degr$. In this special case, the components of the electric and magnetic torques, respectively denoted by $K^{\rm E}$ and $K^{\rm B}$, are given to lowest order in $a=R/r_{\rm L}$ by a decomposition onto the orthonormal basis $(\mathbf{e}_{\rm x},\mathbf{e}_{\rm y},\mathbf{e}_{\rm z})$ with components

[TABLE]

The electric torque along the $z$ axis always vanishes $K_{\rm z}^{\rm E}=0$.

We compare these analytical expressions (8) to the torque found by integration of eq. (7) directly from the simulations. The results are shown individually for the electric torque $K_{\rm x}^{\rm E}$ along the $x$ axis in Fig. 9, the electric torque $K_{\rm y}^{\rm E}$ along the $y$ axis in Fig. 10, the magnetic torque $K_{\rm x}^{\rm B}$ along the $x$ axis in Fig. 11, the magnetic torque $K_{\rm y}^{\rm B}$ along the $y$ axis in Fig. 12 and magnetic torque $K_{\rm z}^{\rm B}$ along the $z$ axis in Fig. 13. To a good accuracy, the simulation results agree with the analytical approximations even for high displacements $\epsilon\approx 0.5$.

After this rather theoretical investigation of the off-centred rotating dipole, showing the induced spindown luminosity perturbation, the kick and additional torque imparted to the star, we finish this paper by a last section about possible consequences for neutron stars in binary systems.

7 Impact on binary neutron star systems

Several hundreds of neutron star binary systems are known to date. The usual channels to form such binaries are well described in several papers like Phinney et al. (1992) or Postnov & Yungelson (2014). Moreover, Tauris et al. (2017) give an excellent review about the more specific double neutron star systems, pointing out their characteristics, formation and evolution. Neutron star binaries are expected to relax to almost circular orbits with very low eccentricities $e\approx 0$ due to mass transfer and tidal circularization. Large eccentricities $e\gtrsim 0.3$ can however be produced by a supernova explosion in double neutron stars when a significant fraction of the binary mass is lost. Here surprisingly, some double neutron stars show low eccentricities requiring alternative binary evolution scenarios. In this last section, we show how the electromagnetic kick produced by an off-centred magnetic dipole can modify the orbital eccentricity, sometimes generating moderate to large eccentricities in neutron star binaries.

Let us assume that both stars, labelled with subscript $1$ and $2$ are subject to an electromagnetic kick respectively denoted by $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$. The equations of motion in the inertial observer frame for this problem are given by

[TABLE]

where $\mathbf{r}=\mathbf{r}_{2}-\mathbf{r}_{1}$ and $\mathbf{r}_{i}$ is the position of star $i$ and dots meaning time derivation. The barycenter, defined by

[TABLE]

feels a net force according to Newton second law

[TABLE]

It therefore accelerates at a constant rate. In order to remove this uniform acceleration not relevant for the orbital parameters evolution, we express the motion in the barycenter frame by introducing two new position vectors

[TABLE]

These vectors satisfy the equations of motion

[TABLE]

The equivalent one body problem then reads

[TABLE]

This tow-body problem in gravitation, subject to an additional acceleration expressed by

[TABLE]

arising from an external body, an external pressure or produced by the stars themselves, is known as the Stark or accelerated Kepler problem (Namouni & Guzzo, 2007). Moreover, if this force is constant in direction and time, interestingly enough, the problem is fully integrable as shown in depth by Lantoine & Russell (2011). They give a complete set of analytical solutions involving parabolic coordinates and elliptic functions for elliptic, parabolic and hyperbolic trajectories. In the present work, we focus only on bound orbits for a constant in direction acceleration $\mathbf{A}$. Important but simple results for eccentricity excitation have been reported by Namouni (2005) for extrasolar planets and for jets by Namouni & Guzzo (2007). In all these cases, the accelerating field is derived from the gradient of a scalar function. Therefore, the time evolution of the orbital parameters are known to satisfy Lagrange planetary equations. A detailed derivation of the formalism and practical results can be found in Beutler et al. (2005). Next we discuss a straightforward application to binaries in circular orbit for which simple analytical expressions have been derived.

Indeed, let us describe the time dependence of the eccentricity $e(t)$ when starting from a circular orbits with $e(t=0)=0$. It varies periodically following a sinusoidal law expressed as

[TABLE]

where $a$ is the semi-major axis, $A=\mathbf{\|A\|}$ is the acceleration produced by the neutron stars themselves due to the electromagnetic kick $F_{i}$, $i_{0}$ is the inclination angle between the acceleration vector $\mathbf{A}$ and the orbital angular momentum vector and $\Omega=\sqrt{G\,(m_{1}+m_{2})/a^{3}}$ is the keplerian frequency. The typical timescale of eccentricity excitation is then

[TABLE]

Note the factor $2$ difference with respect to the $\sin$ argument because of the absolute sign, dividing the period by the same factor 2. This order of magnitude is easily reproduced by noting that generating a significant eccentricity requires the accelerating field $A$ to produce orbital velocities of the order of the keplerian velocity $v=\Omega\,a$, thus equating $A\,T_{e}=v$ which agrees with Eq. (17) within a factor unity. This timescale holds if the excitation $A$ is constant in time. However, the acceleration $\mathbf{A}$ depends on the electromagnetic kick $\mathbf{F}$ which itself depends linearly on the spindown luminosity $L$ because $F\propto L/c$, see previous sections. For a magnetic dipole rotating in vacuum, a braking index of $n=3$ (which still hold to good accuracy for an off-centred dipole) leads to a decreasing luminosity according to

[TABLE]

where $\tau_{\rm c}=P/2\dot{P}$ is the characteristic electromagnetic spindown time scale, the characteristic age of the pulsar. The eccentricity then follows from

[TABLE]

From the spindown decrease in Eq. (18), we easily integrate the action of the acceleration acting on a characteristic time scale $\tau_{\rm c}$ to find

[TABLE]

After a time equal to the characteristic age, the acceleration becomes negligible. Thus the estimate given for a constant in time acceleration is a good guess if time $t$ is replaced by $\tau_{\rm c}/2$. The actual binary eccentricity is therefore

[TABLE]

Eq. (16) shows that the maximal eccentricity achieved after a full excitation period $T_{e}$ is $\sin i_{0}$. This time scale must be compared to other typical time scales like the true age of the binary and the electromagnetic spindown time scale. It can be shown that the eccentricity depends only on $P$ and $P_{\rm orb}$ but not on $\dot{P}$. Straightforward calculations give

[TABLE]

$m_{1}$ being the mass of the pulsar and $m_{2}$ the mass of its companion. For low eccentricities, this expression reduces to

[TABLE]

showing the dependence

[TABLE]

The eccentricity in Eq. (21) can be compared with observations through the ATNF Pulsar Catalogue (http://www.atnf.csiro.au/research/pulsar/psrcat/) of Manchester et al. (2005). For concreteness, we set the inclination angle to $i_{0}=90\degr$, representing the most favourable geometry to efficiently excite the eccentric orbit. For the accelerating electromagnetic kick, we set an average of $\epsilon\,\cos\alpha\,\sin\alpha\,\sin\beta=0.01$ and the total mass of the binary to $m_{1}+m_{2}=2.5\,M_{\odot}$. The neutron star radius is 12 km and its moment of inertia is $I=10^{38}$ kg.m2.

Fig. 14 shows the period-eccentricity relation $(P,e)$ for the neutron star binaries contained in the ATNF catalog in red, compared to the model predictions in blue. Pulsars with period less than $P<0.1$ s are reasonably well reproduced with largest eccentricities $e\lesssim 0.1$, but longer period pulsar eccentricities predictions are orders of magnitude lower than those observed. There are obviously other means to achieve very large eccentricities $e\lesssim 1$ like for instance a supernova explosion leading to a substantial mass loss in the system. All possible perturbations of orbital parameters must be combined to extract realistic binary geometries. However this is out of the scope of this work that only emphasizes the non negligible role of the electromagnetic kick. Fig. 15 compares the observed and predicted orbital period-eccentricity relation $(P_{\rm orb},e)$. Observed eccentricities for orbital periods less than $P_{\rm orb}<100$ days are well reproduced, but longer orbital period eccentricities here also are orders of magnitude lower than observed. Obviously again, the formation and evolution of the binary system must account for those large eccentricities $e$. Finally, the eccentricity distribution function is plotted in Fig. 16, showing the range of eccentricities achieved by the electromagnetic kick in blue, starting from a perfect circular orbit and compared to observations in red. The peak arises at $e\approx 10^{-5}-10^{-4}$ with minimum value of $e_{\rm min}\approx 10^{-12}$ and maximum value of $e_{\rm max}\approx 10^{-1}$. The distribution function is well reproduced in the range $[10^{-7},10^{-2}]$. The deficit in high eccentricities for large spin periods and large orbital periods is clearly identified. Amplifying the excitation acceleration $A$ will produce eccentricities closer to 1. Indeed, when $A$ is increased by a factor 10, i.e. $\epsilon\,\cos\alpha\,\sin\alpha\,\sin\beta=0.1$, the eccentricity distribution function shown in Fig. 17 is shifted one order of magnitude to larger $e$ with $e_{\rm max}\gtrsim 0.1$. Precise values of the acceleration $A$ are geometry dependent and only barely constrained. Nevertheless, already small off-centered dipoles with $\epsilon\,\cos\alpha\,\sin\alpha\,\sin\beta\ll 1$ can accounted for the majority of low eccentricity neutron star binaries.

A complete scenario of binary neutron star formation and evolution must account for not only the dramatic birth event but also for secular change in the orbits due to continuous perturbations induced by the stars themselves. We hope that the present study will help to better understand the formation of binary neutron star systems.

8 Conclusions

We performed accurate time-dependent numerical simulations of an off-centred rotating dipole in vacuum, including contributions from high order multipoles, going beyond the quadrupole expansion as done in previous analytical works. As a diagnostic, we showed magnetic field line structures, spin-down luminosities, induced electromagnetic forces and torques. We demonstrated that the analytical approximations taking into account only perturbations up to the magnetic quadrupole are reliable for quickly computing with satisfactory accuracy the geometric dependence of the spin-down luminosity, the electromagnetic kick and torque according to the orientation and location of the magnetic dipole with respect to the centre of the star.

Other models able to produce naturally multipolar components are those incriminating an inhomogeneous magnetization inside the star. It is well known that a uniformly magnetized sphere produces in vacuum a perfectly centred dipole. If the homogeneity is broken, multipolar components will easily arise in vacuum outside this sphere. Therefore, such models offer an interesting alternative to the off-centred dipole, avoiding a singular magnetic field at the location of the magnetic moment, replacing it by a smooth magnetization distribution over the whole stellar volume.

Electromagnetic kicks and torques are important for understanding the secular evolution of isolated neutron stars, their proper motion and precession, but also for binary neutron stars. Indeed, when in a binary neutron star system, the electromagnetic force permanently perturbs the keplerian orbit and could explain the origin of high eccentricities in several of those binaries as shown in the previous section. The distribution of low eccentricities is well reproduced by an electromagnetic kick acting on the typical spin-down time scale. The electromagnetic torque leads to precession of an isolated pulsar that could be observed in the secular evolution of radio pulse profiles and maybe at higher energies like X-rays and gamma-rays. Such investigations however require knowledge about the radiation mechanisms and is left for future work.

An important extension to the present work is the inclusion of pair plasma within the magnetosphere, screening the electric field component parallel to the magnetic field in the force-free regime. Such models will be computed with our pseudo-spectral Maxwell solver in an upcoming paper.

Acknowledgements

This work has been published under the framework of the IdEx Unistra and benefits from a funding from the state managed by the French National Research Agency as part of the investments for the future program. This work is also supported by the CEFIPRA grant IFC/F5904-B/2018.