EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 700 (August 2025) A&A, 700 (2025) A103 Full HTML
Open Access
Issue
A&A
Volume 700, August 2025
Article Number A103
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202554426
Published online 08 August 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The magnetic field of neutron stars provides the main route for their observation, via the generation of coherent radio emission (as is the case in rotation powered pulsars; Philippov et al. 2020), through its decay (Pons et al. 2007), or through particle acceleration and crust heating (Beloborodov 2009), which is relevant to thermally emitting X-ray neutron stars. Moreover, explosive and transient events in the form of flares (Mazets et al. 1979; Hurley et al. 1999, 2005), bursts (Coti Zelati et al. 2018), giant pulses (Karuppusamy et al. 2010), changes in timing properties (Woods et al. 2002; Archibald et al. 2015; Pintore et al. 2016; Scholz et al. 2017; Hu & Ng 2019; Levin et al. 2019), nulling, moding (Rankin et al. 2013; Lyne et al. 2010), and even some forms of timing noise (Tsang & Gourgouliatos 2013) are related to a greater or lesser extent to their magnetic field. This is mediated either through instabilities and major reconfigurations of its structure (Parfrey et al. 2012, 2013; Gourgouliatos & Hollerbach 2016; Gourgouliatos & Esposito 2018; Gourgouliatos & Pons 2019) or, possibly, through electric discharges (Soglasnov et al. 2004; Chen & Beloborodov 2014) and secular variations in the magnetic field (Pons et al. 2012).

The central role played by the observable properties of neutron stars has motivated thorough studies of the magnetic field structure both at the interior (Braithwaite & Nordlund 2006; Lander & Jones 2009; Ciolfi & Rezzolla 2013; Uryū et al. 2014; Armaza et al. 2015; Suvorov & Glampedakis 2023) and at the exterior. The first studies were based on axisymmetry (Contopoulos et al. 1999; Uzdensky 2003; Goodwin et al. 2004; Gruzinov 2005; Timokhin 2006) and then expanded to three dimensions (Spitkovsky 2006; Kalapotharakos & Contopoulos 2009), the use of techniques ranging from relaxation to force-free electrodynamics and particle in cell (Philippov & Spitkovsky 2014; Philippov et al. 2015; Kalapotharakos et al. 2018), and finally to approaches based on machine learning (Stefanou et al. 2022; Dimitropoulos et al. 2024). Despite the high level of complexity of these studies, there have been limited attempts to simultaneously address both the internal and the external field structure, and this has been done by simplifying the configuration, i.e. by considering a non-relativistic magnetosphere (Glampedakis et al. 2014; Fujisawa & Kisaka 2014; Akgün et al. 2016, 2018a,b) or a light cylinder sufficiently far from the surface of the star while focusing on the field structure at the interior and the near zone of the star (Uryū et al. 2014; Pili et al. 2015). Thus, the question of the impact of the internal field on the external and vice versa remains unresolved.

In this work we address the issue of linking the magnetic field of the interior of the star to the magnetosphere. We simultaneously solved for the internal field so that it obeys a barotropic magnetohydrodynamic (MHD) equilibrium and the magnetospheric field is in a relativistic force-free equilibrium. The system contains a toroidal field that is either in the form of a twisted torus confined within the star or extends at the magnetosphere in the region of the field lines that close within the light cylinder. We provide a global solution for the magnetic field structure of an axisymmetric relativistically rotating neutron star.

In Sect. 2 we present the mathematical setup and formulate the equations that describe the system. In Sect. 3 we present our solution strategy. In Sect. 4 we present the simulations and the results of the calculations. We discuss the properties of the solutions in Sect. 5, present some possible astrophysical applications in Sect. 6, and conclude in Sect. 7.

2. Mathematical setup

Here we present the derivation of the equations that will be solved for the equilibrium of the magnetic field at the interior and exterior of the neutron star. We used cylindrical coordinates (R, ϕ, z) and assumed axial symmetry throughout the system. The length units were chosen so that the light cylinder RLC = 1 = c/Ω, where c is the speed of light and Ω the angular frequency of the star.

We expressed the magnetic field through two scalar functions, Ψ(R, z) and I(R, z), as follows:

B = ∇ Ψ × ∇ ϕ + I ∇ ϕ , $$ \begin{aligned} \boldsymbol{B} = \nabla \Psi \times \nabla \phi + I\nabla \phi , \end{aligned} $$(1)

where ∇ ϕ = ϕ ̂ / R $ \nabla\phi = \hat{\phi}/R $. The function Ψ(R, z) represents the poloidal flux, while I(R, z) is proportional to the poloidal electric current passing through a circular disk parallel to the horizontal plane, whose radius is R and is centred on the axis of symmetry at a distance z from the horizontal plane. From Gauss’s law, the magnetic field has zero divergence, which is satisfied by construction given the definition of Eq. (1). We used the above expression in the appropriate framework to describe the equilibrium at the stellar interior and exterior.

2.1. Internal field

In an MHD equilibrium, the force equilibrium is given by the following expression (Reisenegger 2009):

1 c ρ ( j × B ) = 1 ρ ∇ P + ∇ Φ , $$ \begin{aligned} \frac{1}{c\rho }\left(\boldsymbol{j}\times \boldsymbol{B}\right) = \frac{1}{\rho }\nabla P + \nabla \Phi , \end{aligned} $$(2)

where j is the electric current density, P is the pressure, Φ is the gravitational potential, and ρ is the stellar density. Under the assumption of barotropicity, pressure P is a function of density ρ, P = P(ρ). Taking the curl of Eq. (2) and using Ampère’s law, we obtain

∇ × ( B × ∇ × B ρ ) = 0 . $$ \begin{aligned} \boldsymbol{\nabla }\times \left(\boldsymbol{B}\times \frac{\boldsymbol{\nabla }\times \boldsymbol{B}}{\rho }\right) = \boldsymbol{0}. \end{aligned} $$(3)

Since the above expression is a curl of a vector, equal to 0, we can set the quantity inside curl equal to a gradient of a scalar function S(R, z):

∇ S = 1 ρ B × ( ∇ × B ) . $$ \begin{aligned} \nabla S = \frac{1}{\rho }\boldsymbol{B}\times (\nabla \times \boldsymbol{B}). \end{aligned} $$(4)

From the requirement of axisymmetry, it is evident that

[ B × ( ∇ × B ) ] ϕ = 0 , $$ \begin{aligned} \left[\boldsymbol{B}\times (\nabla \times \boldsymbol{B})\right]_{\phi } = 0, \end{aligned} $$(5)

and thus any poloidal current must be parallel to the poloidal magnetic field, which leads to the following relation:

∇ Ψ × ∇ I = 0 ⇔ I = I ( Ψ ) . $$ \begin{aligned} \nabla \Psi \times \nabla I = 0 \Leftrightarrow I = I(\Psi ). \end{aligned} $$(6)

The poloidal component of Eq. (4) satisfies the following equation:

( ∂ 2 Ψ ∂ R 2 − 1 R ∂ Ψ ∂ R + ∂ 2 Ψ ∂ z 2 + I d I d Ψ ) ∇ Ψ + R 2 ρ ( r ) ∇ S = 0 , $$ \begin{aligned} \left(\frac{\partial ^2\Psi }{\partial R^2} - \frac{1}{R}\frac{\partial \Psi }{\partial R} + \frac{\partial ^2\Psi }{\partial z^2} + I\frac{\mathrm{d}I}{\mathrm{d}\Psi }\right)\nabla \Psi + R^2 \rho (r) \nabla S =\mathbf{0}, \end{aligned} $$(7)

where we have assumed that the density is a function only of the spherical radius r = R 2 + z 2 $ r = \sqrt{R^2+z^2} $. From Eq. (7), it is evident that the gradient of S is parallel to that of Ψ, ∇Ψ||∇S, leading to the conclusion that S = S(Ψ). Eventually, the barotropic equilibrium equation is expressed as follows:

Δ ∗ Ψ + I I ′ + R 2 ρ ( r ) S ′ = 0 , $$ \begin{aligned} \Delta ^*\Psi + II^{\prime } + R^2 \rho (r)S^{\prime } = 0, \end{aligned} $$(8)

where

Δ ∗ = ∂ 2 ∂ R 2 − 1 R ∂ ∂ R + ∂ 2 ∂ z 2 $$ \begin{aligned} \Delta ^* =\frac{\partial ^2}{\partial R^2} - \frac{1}{R}\frac{\partial }{\partial R} + \frac{\partial ^2}{\partial z^2} \end{aligned} $$(9)

is the Grad-Shafranov operator and a prime denotes differentiation with respect to Ψ.

We approximated the mass density, ρ, of the star using the relation

ρ ( r ) = ρ 0 r ns 2 − r 2 r ns 2 , $$ \begin{aligned} \rho (r) = \rho _{0}\frac{r_{\rm ns}^2-r^2}{r_{\rm ns}^2}, \end{aligned} $$(10)

where rns = 0.1 is the radius of the star and ρ0 is the density at the centre of the star. Equation (10) corresponds to the solution of the mass density of a neutron star, for an n = 1 polytrope, where P = ρ2 (Lattimer & Prakash 2001).

Despite the functional freedom of S(Ψ), as shown by Glampedakis & Lasky (2016), there are some restrictions on the applicable forms. Here, we assumed that S is a linear function of Ψ and thus that its derivative appearing in the final partial differential equation is a constant, and we set its value S′ = 1. This form is compatible with the requirements of the aforementioned work. The choice of a linear form for S(Ψ) restricts the solutions. However, it allows magnetic fields in the form of a dipole for I = 0 and assuming vacuum external boundary conditions, something that is not possible for non-linear expressions of S in Ψ (Gourgouliatos et al. 2013). Moreover, given the large number of free parameters appearing in the problem, we chose to adopt the linear form as a first approach and will explore the wider parameter space in future work.

2.2. External field

Considering the relativistically corotating, plasma-filled magnetosphere, we assumed that electromagnetic forces prevail over gravitational, inertial, and pressure gradient forces. Therefore, a system in equilibrium must have zero Lorentz force, specifically

ρ q E + 1 c j × B = 0 , $$ \begin{aligned} \rho _q \boldsymbol{E} + \frac{1}{c}\boldsymbol{j}\times \boldsymbol{B} = 0, \end{aligned} $$(11)

where ρq is the electric charge density. By assuming ideal MHD, the expression of Eq. (1) for the magnetic field and Ohm’s law for the electric field corresponding to a corotating magnetosphere, where the rotation axis coincides with the dipole axis, is

E = − Ω R c ϕ ̂ × B . $$ \begin{aligned} {\boldsymbol{E}} = -\frac{\Omega R}{c} {\hat{\phi }} \times \boldsymbol{B}. \end{aligned} $$(12)

From this, we obtain the axisymmetric pulsar equation (Scharlemann & Wagoner 1973):

( 1 − R 2 ) ( ∂ 2 Ψ ∂ R 2 − 1 R ∂ Ψ ∂ R + ∂ 2 Ψ ∂ z 2 ) − 2 R ∂ Ψ ∂ R = − I ( Ψ ) I ′ ( Ψ ) , $$ \begin{aligned} (1 - R^2)\left(\frac{\partial ^2\Psi }{\partial R^2} - \frac{1}{R}\frac{\partial \Psi }{\partial R} + \frac{\partial ^2\Psi }{\partial z^2}\right) -2R\frac{\partial \Psi }{\partial R} = -I(\Psi )I^{\prime }(\Psi ), \end{aligned} $$(13)

where the lengths have been suitably normalized to the light-cylinder radius.

Considering the relativistic Grad-Shafranov operator,

Δ ∗ = ( 1 − R 2 ) ( ∂ 2 ∂ R 2 − 1 R ∂ ∂ R + ∂ 2 ∂ z 2 ) − 2 R ∂ ∂ R , $$ \begin{aligned} \Delta _*= (1 - R^2)\left(\frac{\partial ^2}{\partial R^2} - \frac{1}{R}\frac{\partial }{\partial R} + \frac{\partial ^2}{\partial z^2}\right) -2R\frac{\partial }{\partial R}, \end{aligned} $$(14)

Eq. (13) can be written as

Δ ∗ Ψ = − I ( Ψ ) d I ( Ψ ) d Ψ · $$ \begin{aligned} \Delta _*\Psi = -I(\Psi )\frac{\mathrm{d}I(\Psi )}{\mathrm{d}\Psi }\cdot \end{aligned} $$(15)

In the standard pulsar solution (Contopoulos et al. 1999), I(Ψ) = 0 in the area of closed field lines, while for the open field it is I(Ψ)≠0, ensuring a smooth transition through the light cylinder region. Magnetic fields of minimal complexity have zero poloidal current inside regions of closed magnetic lines, based on the hypothesis that pulsar fields are not inherently twisted. The implementation of a non-zero current, I(Ψ)≡Itw ≠ 0, in this area, along with the concept of a twisted magnetosphere, is prevalent in magnetar models and, subsequently, in the strongly magnetized systems studied here. We refer to this current as the twist current, Itw, and it will appear in the right-hand-side of Eq. (13).

3. Solution strategy

We studied two main families of models. In the first, the toroidal twist field is allowed to exist only within the star. We refer to this type of twist as an internal twist. In the second family of models, closed magnetospheric fields lines and field lines at the interior of the star that are connected to the closed magnetospheric fields lines can be twisted. We refer to this twist as the global twist, since it extends both at the stellar interior and the magnetosphere. We remark that the open field lines have, in either model, a toroidal field, which is related to the smooth crossing of the light cylinder.

The solution domain is partitioned in five regions, as shown in Fig. 1. Region (I) is the part of the magnetosphere containing open magnetic field lines that cross the light cylinder. Region (II) is the area of the closed magnetic field lines in the magnetosphere containing closed magnetic field lines. Regions (III), (IV), and (V) are all the stellar interior, where region (III) contains field lines that cross the stellar surface and connect to the open magnetospheric field lines, region (IV) contains magnetic field lines that connect to the closed field lines of the magnetosphere, and region (V) contains field lines that close within the star; thus, they are not connected to the magnetosphere and form closed toroids.

thumbnail Fig. 1.

Various regions of the domain, labelled (I) to (V). The field lines are depicted in black, the last open field line of the magnetosphere in red, the stellar surface in green, and the light cylinder in dashed blue. The top-right inset provides a zoomed-in view of the star.

The magnetospheric field in region (I) requires a poloidal return current, Irt, for the smooth crossing of the light cylinder. The magnetospheric field in region (II) satisfies the relativistic force-free equation; however, unlike region (I), there is freedom in the form of the current as no field lines cross the light cylinder and thus there is no requirement for a return current. However, it is possible that these field lines are twisted. The magnetic field at the interior, regions (III), (IV), and (V) satisfies Eq. (8) corresponding to the barotropic equilibrium. We note that, formally, the cross product of the poloidal magnetic field and the electric current in region (III) exert a force in the ϕ direction, which is required for the spin-down of the star arising from the penetration of the return current into the stellar crust (Karageorgopoulos et al. 2019). In our approach we neglected the impact of this current in the interior and found its equilibrium through Eq. (8) by setting I = 0. Regions (II), (IV), and (V) may contain a poloidal current as a result of the twisting of the magnetic field, Itw. In summary, the subsequent relations describe the field in the five aforementioned regions:

Δ ∗ Ψ = − I rt I rt ′ (I) Δ ∗ Ψ = − I tw I tw ′ (II) Δ ∗ Ψ + R 2 ρ ( r ) S ′ = 0 (III) Δ ∗ Ψ + I tw I tw ′ + R 2 ρ ( r ) S ′ = 0 . (IV), (V) $$ \begin{aligned}&\Delta _*\Psi = -I_{\rm rt} I_{\rm rt}^{\prime }&\text{(I)} \nonumber \\&\Delta _*\Psi = -I_{\rm tw} I_{\rm tw}^{\prime }&\text{(II)} \nonumber \\&\Delta ^*\Psi + R^2 \rho (r)S^{\prime } = 0&\text{(III)} \\&\Delta ^*\Psi + I_{\rm tw}I_{\rm tw}^{\prime } + R^2 \rho (r)S^{\prime } = 0.&\text{(IV),} \text{(V)}\nonumber \end{aligned} $$(16)

The relativistic Grad-Shafranov operator appears explicitly in regions (I) and (II) of the magnetosphere, whereas all of the interior of the star is characterized by the non-relativistic Grad-Shafranov operator. This issue can be resolved even under the assumption that all areas of the star’s magnetic field are solved through the relativistic Grad-Shafranov operator. Given that the stellar surface radius is rns = 0.1RLC, the star’s rotation exerts a minor influence on the internal structure of the magnetic field, as this would correspond to a maximum correction of (rns/RLC)2 = 10−2. We have confirmed that changes in the magnetic flux function, when the above correction is performed, do not surpass 1%. Furthermore, the return current is in general much weaker than the twist current, especially if the latter has a remarkable effect on the magnetosphere, and thus not including it in the calculation of the internal equilibrium does not alter the structure of the internal field.

Regarding the boundary conditions that will be implemented, we proceeded as follows. We assumed that the magnetic flux function, Ψ, is zero along the z-axis:

Ψ ( R = 0 , z ) = 0 . $$ \begin{aligned} \Psi (R= 0, z) = 0. \end{aligned} $$(17)

The value of Ψ at the star’s surface is not defined as that of a dipole; it is determined by solving Eq. (16) within the star. In this problem, the light cylinder is positioned at a distance of 10 stellar radii (RLC = 10rns). The chosen distance of the light cylinder corresponds to an angular velocity of the neutron star of Ω = 3000 rad/sec, indicating that we are formally examining a millisecond pulsar with a spin frequency of approximately 500 Hz.

The magnetic field is symmetric with respect to the northern and southern hemispheres; hence, the subsequent boundary conditions are applied at the equator:

∂ z Ψ ( R < R T , z = 0 ) = 0 , $$ \begin{aligned} \partial _z\Psi (R< R_{\rm T}, z=0) = 0, \end{aligned} $$(18)

where RT denotes the equatorial distance of the outermost closed field line of the magnetosphere. Subsequent to this intersection and along the equatorial plane, the equatorial current sheet of the magnetosphere is characterized by the expression

Ψ ( R ≥ R T , z = 0 ) = Ψ ( R T , 0 ) . $$ \begin{aligned} \Psi (R\ge R_{\rm T}, z=0) = \Psi (R_{\rm T}, 0). \end{aligned} $$(19)

Although most works assume that RT lies exactly on the light cylinder, it may in fact be located between 0.8RLC and 0.9RLC (Contopoulos et al. 2024). When closed magnetic field lines are twisted more than a minimum amount, physically acceptable solutions require RT to be moving towards the star (Ntotsikas et al. 2024), which is accounted for in our models.

To guarantee the smooth crossing of the field lines through the light cylinder, we imposed continuity on the function and its derivative,

Ψ ( R = 1 + , z ) = Ψ ( R = 1 − , z ) , $$ \begin{aligned} \Psi (R = 1^{+},z) = \Psi (R = 1^{-},z), \end{aligned} $$(20)

and imposed the regularization condition

∂ R Ψ ( R = 1 − , z ) = ∂ R Ψ ( R = 1 + , z ) = 1 2 I ( Ψ ) I ′ ( Ψ ) . $$ \begin{aligned} \partial _R\Psi (R=1^{-},z) = \partial _R\Psi (R=1^{+},z) = \frac{1}{2}I(\Psi )I^{\prime }(\Psi ). \end{aligned} $$(21)

The last equation also relates the value of the R derivative of Ψ to I. For the outer region and the open field lines of the magnetosphere, we implemented the split monopole boundary conditions to obtain radial behaviour at the boundary of the integration grid zmax, Rmax. Formally, boundary conditions of this type were applied at an infinite distance from the surface of the star, in which case the magnetic lines eventually acquire a radial shape. However, applying split-monopole conditions to the boundaries of our grid does not change the results compared to a run with a twice larger numerical domain, but it does significantly increase the convergence speed of the computation. Relation (21) enables us to determine the form of I(Ψ) for magnetic flux function values less than Ψ0, where Ψ0 represents the magnetic flux function at the first closed magnetic field line.

We note that a large fraction of the magnetospheric return current enters the star through a current sheet flowing on the separatrix. Such singularities cannot be handled by the finite difference method used here; therefore, we approximated the current sheet using a Gaussian function with a width of 1 × 10−3Ψ0 centred at Ψ = 0.99Ψ0, ensuring that the entire return current flows through the open field lines. Unlike the solutions dealing exclusively with the magnetosphere, in our approach the field emerging from the star is not a pure dipole, and thus the maximum flux function may have a different value depending on the twist imposed. To avoid such discrepancies, we appropriately normalized the magnetic flux function so that the first closed field line for the α = 0 model corresponds to Ψ = 1.23, therefore enabling a direct comparison with the standard solutions (Timokhin 2006), and we used this normalization for the rest of the solutions. Here, α is a parameter that determines the ratio of the electric current distribution to the difference of the flux function, Ψ, minus its value at the first closed magnetic field line.

With regard to the twist current, we focused on two cases. The first is the globally twisted model in which the closed magnetospheric field lines carry poloidal current (region II) and are connected to the star (region IV):

I ( Ψ ) = { I rt , (I) α ( Ψ ( R , z ) − Ψ 0 ) , (II),(IV),(V) 0 , (III) . $$ \begin{aligned} I(\Psi ) = {\left\{ \begin{array}{ll} I_{\rm rt},&\text{(I)} \\ \alpha (\Psi (R,z)-\Psi _0),&\text{(II),(IV),(V)} \\ 0,&\text{(III)}. \end{array}\right.} \end{aligned} $$(22)

Alternatively, in the internally twisted model, any poloidal current corresponding to an internal poloidal field in the star closes within the star, and thus there is no poloidal current in the closed field lines (region II), nor in the field lines within the star that are connected to the closed magnetospheric field lines (region IV):

I ( Ψ ) = { I rt , (I) 0 , (II),(III),(IV) α ( Ψ ( R , z ) − Ψ ssf ) , (V) . $$ \begin{aligned} I(\Psi ) = {\left\{ \begin{array}{ll} I_{\rm rt},&\text{(I)}\\ 0,&\text{(II),} \text{(III),} \text{(IV)} \\ \alpha (\Psi (R,z)-\Psi _{\rm ssf}),&\text{(V)}. \end{array}\right.} \end{aligned} $$(23)

In all cases, Irt is the appropriate return current determined by the smooth light crossing condition. The poloidal twist current in the star is a linear expression given by the flux function, Ψ, reduced by either its value at the last closed field line (Ψ0) or its value at the toroidal loop that touches the stellar surface, for which Ψssf = Ψ(rns, 0), i.e. the value of the flux function on the equator of the star.

Matching of the stellar and magnetospheric solutions was done by requiring continuity on Ψ and I. The solution in the stellar interior satisfies Δ ∗ Ψ + I tw I ′ tw + R 2 ρ S ′ = 0 $ \Delta^{*} \Psi+I_{\mathrm{tw}}I\prime_{\mathrm{tw}}+R^2 \rho S\prime=0 $, while the corresponding one in the magnetosphere is Δ ∗ Ψ + I tw I ′ tw = 0 $ \Delta_{*} \Psi + I_{\mathrm{tw}}I\prime_{\mathrm{tw}}=0 $. The two equations differ in terms of the form of the Grad-Shafranov operator (non-relativistic vs relativistic) and on the term involving the density. We have already commented on the impact of the non-relativistic term of the Grad-Shafranov operator, which leads to differences that scale with the square of the ratio of the neutron star radius to the light cylinder, which in this case is 0.01. The difference due to the density term is minimized as the density formally drops to zero on the surface of the star (see Eq. (10)), and thus the equation is satisfied on both sides of the stellar boundary. We note, however, that due to discretizations and the cylindrical coordinates we have adopted, the density on the surface has some discontinuities, though they do not lead to sharp transitions in the field structure.

In summary, two families of models of magnetic field configurations are developed. The first one (Eq. (22)) corresponds to the following state: regions (II), (IV), and (V) contain a current associated with the internal twist of the magnetic field lines that also leaks out in the closed magnetospheric field lines, while open magnetospheric lines carry the essential return current from the outer magnetosphere to the star’s interior. In the second one (Eq. (23)) there is only a poloidal twist current in the field lines that close within the star, region (V) without crossing the surface of the magnetosphere, while the open field lines, region (I) carry the essential return current. The choice of the linear expression for the twist current is restrictive given the wide parameter space; however, we have chosen to limit ourselves to a simple functional form and explore as a main parameter the ratio of the poloidal current to the flux function. The chosen expressions, I = α(Ψ(R, z)−Ψ0) or I = α(Ψ(R, z)−Ψssf), guarantee that the toroidal field vanishes along this boundary. There is thus no magnetic field discontinuity because of the twist nor an additional current sheet that would have broken the north-south symmetry. A current sheet discontinuity is the one that already exists due to the presence of the return current that closes along the separatrix between the closed and the open magnetic field lines.

4. Simulations and results

We solved Eq. (16) for different values of the parameter α, which determines the current ratio over the magnetic field, applying the simultaneous relaxation method and considering the two cases of global and internal twist. Our simulations were conducted using a FORTRAN 90 framework. As the system is axisymmetric, the functions depend only on R and z. We divided the numerical domain into R and z groups with a typical resolution of 800 × 800. We also ran higher-resolution simulations, 1600 × 1600, to verify the numerical convergence of our results, and a change of less than 0.5% was found. The derivatives were computed numerically using the finite difference method; a central difference scheme was used for the first derivatives, and a three-point stencil for the second derivatives. The magnetosphere code is based on a refined version of the code from Gourgouliatos & Lynden-Bell (2019), which was also used in Contopoulos et al. (2023). The interior code has been developed from the first principles for this particular application and has been benchmarked against Gourgouliatos et al. (2013).

Figure 2 shows magnetic field solutions for the global twist model, for which Itw is applied in the area of the closed magnetic field lines of the magnetosphere, regions (II), (IV), and (V), with coefficients ranging from α = 0.0 to α = 8.0. As the parameter α increases, we find that the current of the closed field lines becomes stronger. Additionally, the first closed field line corresponds to a higher value of Ψ, and the equatorial current sheet must end closer to the star. We present these effects for various values of α in Table 1, which summarizes the main simulation results. Regarding the magnetic field in the interior, we find that the region of the closed field lines shrinks once more twisted is imposed. We increased the value of α up to 10, as beyond this value the twisted region is rather small and under-resolved. Moreover, the innermost point of the equatorial current sheets is twice the radius of the star for the maximum value of αimplemented.

thumbnail Fig. 2.

Magnetic field of a relativistically rotating neutron star and its magnetosphere for the global twist model. The field lines are shown in black, the stellar surface in green, and the light cylinder with a dashed blue line at R = 1; the current sheet is also shown. The red line represents the first closed field line of the untwisted case (panel a) and is shown in subsequent panels for comparison. A poloidal current is applied in the area of the closed magnetic field lines with coefficients α = 0.0,  3.0,  6.0, and 8.0. The colour bar indicates the value of II′ across several parts of the star and the magnetosphere. The maximum value of II′ in the twisted field lines is saturated in colour to allow the depiction of its structure in the rest of the system. The top-right inset is shows a zoomed-in area of the star.

Table 1.

Simulation results for the global twist model.

We also solved for the internally twisted model (see Table 2 and Fig. 3). As the field is only twisted in the interior, the twisted region tends to shrink with increasing α, as in the globally twisted case. No displacement of the innermost point of the current sheet is required here since the twist is only at the interior and thus mainly affects the magnetic flux emanating from the star and its multipolar structure, with higher multipoles appearing for stronger toroidal fields. These effects impact both the flux value of the first closed field line and the structure of the magnetic field in the magnetosphere.

Table 2.

Simulation results for the internal twist model.

thumbnail Fig. 3.

Magnetic field of the internal twist model. Panels a and b correspond to α = 5.0 and 10.0, respectively. The various field lines and colour bars are shown as in Fig. 2.

We note that for both families of solutions, increasing α leads to a larger amount of flux emerging from the star, Ψssf, and a larger maximum value of flux at the interior of the star, Ψmax; the two other quantities in the equation (density profile and S) stay the same. We note, however, that the changes in the maximum flux and the total flux emerging from the surface of the star are small and do not scale with the increase in Ψ0. Thus, the primary cause of an increasing fraction of open field lines is the change in the magnetosphere. Furthermore, we note that while the form of the flux function on the stellar surface is dominated by a dipolar term, higher multipoles co-exist, leading to a surface magnetic field profile that is not a pure dipole.

Comparing the two families of solutions, while the trend is qualitatively the same, we find that the quantitative differences between the two cases are rather pronounced: for the same values of α, the fraction of the open field lines is substantially higher once twist currents flow through the magnetosphere. This has a major impact on various physical properties of the neutron star, ranging from the spin-down power to the twist of the magnetosphere and the polar cap opening angle.

5. Discussion

Here we address the consequences of increasing the twist for the main physical properties of neutron stars. We focused on the spin-down luminosity, the overall twist in the magnetosphere, and the polar cap opening angle.

5.1. Spin-down luminosity

The general trend is that an increase in the twist current leads to a larger fraction of open field lines. For small values of α ≤ 1, the two models give the same result (Fig. 4). However, once α exceeds this value, the global twist model has a much higher value of Ψ0. This is also clearly visible in Figs. 2 and 3, where the first closed field line of the untwisted model (shown in red) is much more drastically displaced from the equator in the globally twisted model than in the internallytwisted one.

thumbnail Fig. 4.

Magnetic flux function of the first closed magnetic field line as a function of α.

The overall spin-down luminosity, which corresponds to the loss of electromagnetic energy from the star’s magnetic field, is expressed by the following relation:

L = 2 ∫ 0 Ψ 0 norm I ( Ψ ) d Ψ . $$ \begin{aligned} L = 2\int _0^{\Psi _0^\mathrm{norm}}I(\Psi )\mathrm{d}\Psi . \end{aligned} $$(24)

Using Eq. (24), we evaluated the star’s spin-down power. In models where the twist current flows in the magnetosphere, the spin-down power can be enhanced by a factor of approximately 16 for α = 10 (Fig. 5), whereas if the twist current closes within the star, for the same value of α the spin-down luminosity is enhanced by a factor of 4. The rapid increase in the spin-down luminosity in the globally twisted model is related to the fact that the current sheet inner edge is closer to the stellar surface; therefore, there is a drastic increase in the number of the open magnetic field lines.

thumbnail Fig. 5.

Spin-down luminosity of the various solutions scaled to the spin-down power of the untwisted solution, as a function of α.

5.2. Magnetospheric twist

Next, we evaluated the twist of a magnetic field line via integration of the following relation:

d R B R = d z B z = R d ϕ B ϕ , $$ \begin{aligned} \frac{\mathrm{d}R}{B_R} = \frac{\mathrm{d}z}{B_z} = \frac{R\mathrm{d}\phi }{B_{\phi }}, \end{aligned} $$(25)

where BR is the radial component of the magnetic field, Bz the axial, and Bϕ the azimuthal. From this we obtained the following:

Δ ϕ = ∫ B ϕ R B R d R . $$ \begin{aligned} \Delta \phi = \int \frac{B_{\phi }}{RB_R}\mathrm{d}R. \end{aligned} $$(26)

The twist can vary in different field lines, i.e. the first closed field lines have zero twist, as defined there by Bϕ = 0, while the field line emerging from the equator of the star has the maximum value of I; however, its extent is formally zero. Therefore, this integration was evaluated on the field line that corresponds to Ψ = 1.1Ψ0norm. We integrated from the surface of the star to the equator; therefore, the results presented in Table 1 represent half the twist values that would be derived by measuring the whole extent of the magnetic field line that emerges from the northern and closes at the southern hemisphere. In agreement with earlier solutions (Lynden-Bell & Boily 1994; Gourgouliatos & Lynden-Bell 2008; Pavan et al. 2009; Akgün et al. 2018a), the twist tends to approach a maximum value of ∼π/2 rather than increasing indefinitely. Therefore, the system can formally approach a split monopole state after a finite amount of twist. We decided to investigate the properties in terms of the parameter α up to a maximum value of α = 10.

Figure 6 shows the twist of the closed magnetic field line region. It is evident that with increasing α, the twist in region (II) becomes higher. However, this region shrinks and is concentrated closer to the star, and the twist at the exterior will eventually vanish.

thumbnail Fig. 6.

Twist of the magnetospheric part of the field line for which Ψ = 1.1Ψ0, as a function of the parameter α.

5.3. Polar cap

We calculated the semi-opening angle of the polar cap by tracing the foot points of the first closed field line corresponding to Ψ0. The location (RΨ0, zΨ0) of its foot points satisfies the following conditions: Ψ(RΨ0, zΨ0) = Ψ0 and RΨ02 + zΨ02 = rns2. The semi-opening angle is given by the following expression:

θ pc = arctan R Ψ 0 z Ψ 0 · $$ \begin{aligned} \theta _{\rm pc} = \arctan \frac{R_{\Psi _0}}{z_{\Psi _0}}\cdot \end{aligned} $$(27)

We plot the semi-opening angle of the polar cap in Fig. 7. We find that as α increases, the semi-opening angle of the polar cap expands for both the internal and global twist. However, there is a difference in magnitude. If the twist current is allowed in the magnetosphere, the polar cap increases very rapidly, and it even becomes a split monopole after some finite twist. In contrast, if the twist current lies entirely within the star and does not populate the closed magnetospheric field lines, changes on the polar cap are minor, increasing from 17.5° to 20°, due to the fact that the magnetic field on the surface is no longer dipolar, but it contains higher multipoles of north-south symmetry.

thumbnail Fig. 7.

Polar cap’s semi-opening angle as a function of α.

6. Applications

The solutions found here illustrate the impact of twist on the global magnetic field structure of a neutron star, which spans from the interior to the exterior. Our models explore the major qualitative difference of the presence of twist in the magnetosphere or exclusively at the interior of the star, and employ an incremental exploration of the parameter space regarding the injection of twist of the field lines. These differences in the magnetic field structure both at the interior and the exterior likely impact the observational behaviour of neutron stars.

We confirm that the inclusion of a toroidal magnetic field strongly affects the spin-down rate. Even if a comparable amount of poloidal magnetic flux crosses the stellar surface, a twisted field will have a more pronounced spin-down rate. We notice major differences in the spin-down rate depending on whether the twist current is confined within the star or whether it populates the magnetosphere. For smaller values of α up to α ≈ 1, the spin-down power scales the same for the two families of models. For values higher than that, we find that there is a drastic difference. In particular, if we allow the twist current to flow in the magnetosphere, for α = 10 the spin-down luminosity is about 16 times higher than in the untwisted case, whereas if the twist current is enclosed in the star it is about 4 times higher compared to the untwisted model. This effect will differentiate the two families of models. This is related to the fact that the equatorial current sheet needs to retract closer to the star if the magnetosphere is twisted, whereas the current sheet can start at the light cylinder in the absence of external twist. This calculation provides further evidence that the twist of the magnetosphere increases the spin-down luminosity of the star (Thompson et al. 2002; Lyutikov 2006; Viganò et al. 2011). It also, quite remarkably, provides evidence that an internal twist leads to a spin-down luminosity even if it only indirectly affects the magnetosphere by altering the lower boundary conditions.

We note that our calculations are constrained by the fact that the light cylinder is only ten times larger than the stellar radius. A smaller star or a more distant light cylinder would allow for a finer exploration of the innermost point of the current sheet and would not stop at this particular value. In the highly twisted case, where the currents flow in the magnetosphere, we further notice that the polar cap becomes larger and the magnetic field is reminiscent of a split monopole even within the light cylinder. Thus, the field lines are much less curved in the poloidal direction, yet there is still curvature due to the toroidal field. This stretching of the magnetic field disfavors the generation of curvature radiation, as a smaller part of the field is curved, as opposed to the previous case.

A critical question is what determines whether the internal twist will also populate the magnetosphere, or if it is only going to be contained in the star, in principle the capacity of the magnetosphere to support twist currents. Twist currents can become much stronger than the spin-down current, which depends on the ratio of the neutron star radius to the light cylinder. They require a larger population of charges in the magnetosphere that is related to the multiplicity of the charges and can untwist within timescales of years (Beloborodov & Thompson 2007), thus providing a transient, but not necessarily explosive, behaviour. As this question is still unresolved, we speculate that magnetospheres with a higher charged particle density are more likely to host such currents, whereas magnetospheres with smaller particle densities are possibly reminiscent of vacuum or minimally twisted magnetospheres. This will have the following consequences: neutron stars with twisted magnetospheres will have a higher spin-down rate, the field lines will have a larger curvature radius, and, as postulated, they will have a higher charge density. This could provide an interpretation of transient magnetar behaviour; since the spin-down rate is higher, the generation of curvature radiation related to radio emission is less likely due to the structure of the field lines, and the presence of charged regions and currents within the magnetosphere may eclipse any radio emission produced near the star (Levin et al. 2019; Lower et al. 2023). We note that such variations on the spin-down power have been noted in simulations of twisted magnetospheres (Parfrey et al. 2012; Ntotsikas et al. 2024) and also in the timing behaviour of pulsars correlated with emission across the electromagnetic spectrum and changes in pulse profile (Urama et al. 2006; Camilo et al. 2007, 2008; Younes et al. 2020a; Lower et al. 2025; Fisher et al. 2024), also accompanied by spectral evolution (Younes et al. 2025).

These models also introduce current sheets that extend within the light cylinder, especially for highly twisted systems. As current sheets are prone to instabilities, they can undergo tearing modes, which leads to dissipation. This effect may be responsible for the release of energy and particle acceleration. Although such effects are likely to occur in pulsars (Comisso et al. 2017; Cerutti & Giacinti 2021; Bransgrove et al. 2023), this model brings them within the light cylinder and very close to the star. The creation of current sheets is intimately related to the twisting of the external field; thus, torque variations due to switches from twisted and untwisted states will be accompanied by a rapid release of energy through explosive events (Archibald et al. 2020). Moreover, the twist, internal or global, affects the size of the polar cap, with larger polar caps associated with a higher twist (Tong 2019), as has been suggested to be the case in sources with flaring episodes (Younes et al. 2020b).

In our simulations, the chosen ratio of the light-cylinder radius over the neutron star radius corresponds to a millisecond pulsar whose frequency is approximately 500 Hz. Although this is much faster than any known magnetar, it may be directly comparable to millisecond magnetars (Dall’Osso & Stella 2022). These sources are newborn magnetars, expected to spin-down rapidly (Çıkıntoğlu et al. 2020). In addition to these sources, an interesting extension would be the scaling of the model to less rapidly rotating sources. In such systems, the light cylinder lies much farther away than the 10 stellar radii we have assumed here, and it is possible that a smaller fraction of the closed field lines will be twisted. However, the inclusion of twist can impact the spin-down properties of lower rotating sources if it extends all the way to the light cylinder. However, in models where the twist current is contained within the star, part of the change in the spin-down power is due to the fact that the field has a multipolar structure on the surface rather than a dipole. This will not affect significantly the field structure near the light cylinder, as higher multipoles decrease rapidly. Nevertheless, even if these effects are mild, a toroidal field affects the overall amount of magnetic flux emerging from the star, thus leading to a higher spin-down rate.

A further remarkable difference between the models in which the twist is contained in the star versus the ones where the magnetosphere is twisted is the volume occupied by the toroidal field within the star. In the former case, the toroidal field is confined in a very small region, within closed loops of magnetic flux; on the contrary, if the twist is allowed to populate the magnetosphere, it reaches much higher latitudes up to approximately 45°, for the models presented here. This expansion of the toroidal field region can have important implications for the magnetic field energy, as the toroidal field can contribute more intensely to the total energy budget, even under equilibrium conditions. On the other hand, in cases where the twisted field was contained in the star, this ratio was rather low (Lander & Jones 2009) unless a special form of the poloidal current was postulated (Ciolfi & Rezzolla 2013). It may also affect the ellipticity of the star, as the toroidal field tends to increase it, and in the fully confined models it can only contribute locally, whereas in this family of models, the toroidal field can be much stronger and extend to a larger region within the star (Haskell et al. 2008).

The question of a hydromagnetic equilibrium of a neutron star with rotation was also posed in the work of Glampedakis et al. (2014). In that work the problem was studied in the non-rotating system, while touching upon the question of rotation, without, however, fully addressing the light cylinder and the relativistic magnetosphere. Our results provide a direct comparison between the two regimes. The obvious difference is in the very structure of the magnetosphere, where the additional physical constraints of the relativistic force-free magnetosphere have been added. In particular, the current flowing along the open field lines is determined by the requirement that the magnetic field crosses the light cylinder smoothly; therefore, solutions where a predetermined current is allowed to flow are not acceptable. The presence of the light cylinder adds an extra length-scale that, in combination with the twist-current, leads to contraction of the closed field line region and the presence of a current sheet closer to the star. Furthermore, the open field lines, especially near the light cylinder and beyond, have a structure of a split monopole: such a field is significantly different from the dipole background assumed there. Notably, split monopole structures can appear in spherical geometry only after a relatively small amount of twist (Lynden-Bell & Boily 1994; Gourgouliatos & Lynden-Bell 2008; Tong 2019) Apart from the expected differences in the magnetosphere, where the actual equilibrium equation is different, there are further, remarkable differences at the interior. In particular, the internal toroidal field is allowed in the equatorial and mid-latitude region, but not near the poles as in the non-rotating solution of Glampedakis et al. (2014). Consequently, the region where the twist current flows, at the interior this time, can maximally be that of the closed field lines, a region determined by the relativistic solution. Overall, the solution with a relativistically rotating magnetosphere is more constrained than the non-rotation inside-out magnetosphere and provides direct estimation of astrophysically relevant quantities for neutron stars including their relativistic magnetosphere.

7. Conclusions

We investigated equilibrium configurations for twisted magnetic fields of neutron stars, taking both the star and the magnetospheres into consideration. Previous studies generally either examined the external magnetosphere without accounting for the stellar interior or presumed a vacuum exterior. Our work highlights a comprehensive method in which the magnetic field, including the stellar and magnetospheric currents, is consistently solved for.

Our results indicate that the inclusion of twist impacts the structure of the neutron star’s magnetic field. Its effect is moderate if the twist is only allowed in the interior of the star and rather drastic if the twist is also in the closed magnetospheric field lines. We have quantified its impact on the spin-down rate, magnetospheric twist, and polar cap opening angle, offering insights into phenomena such as transitions between pulsar and magnetar states or intermittent sources. This work has provided the framework to extract observable parameters and to eventually quantify these effects.

Future research should progress to solutions beyond that of a barotropic equilibrium and simple linear forms for S(Ψ) and the twist current, Itw. Moreover, it is essential to adopt more realistic models in which both the crust and the core are accounted for, as the relevant equilibrium equations are rather different. Furthermore, progress from axisymmetric models to three-dimensional simulations may provide a more precise understanding of magnetospheric processes.

An in-depth examination of the communication between the star’s inside and exterior, particularly over its surface, is crucial. These endeavours will be essential for addressing fundamental inquiries regarding the evolution of magnetic fields in neutron stars and comprehending the relationship between magnetospheric distortions and observable phenomena, such as timing variations in magnetars and even normal sources.

Data availability

The datasets generated during the current study are available from the corresponding author upon reasonable request.

Acknowledgments

We are grateful to an anonymous referee whose constructive comments improved this paper. We are grateful to I. Contopoulos and S. Lander for discussions during the preparation of this manuscript. This work was supported by computational time granted by the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS – under project ID pr017008/simnstar2.

References

  1. Akgün, T., Miralles, J. A., Pons, J. A., & Cerdá-Durán, P. 2016, MNRAS, 462, 1894 [Google Scholar]
  2. Akgün, T., Cerdá-Durán, P., Miralles, J. A., & Pons, J. A. 2018a, MNRAS, 474, 625 [Google Scholar]
  3. Akgün, T., Cerdá-Durán, P., Miralles, J. A., & Pons, J. A. 2018b, MNRAS, 481, 5331 [CrossRef] [Google Scholar]
  4. Archibald, R. F., Kaspi, V. M., Ng, C. Y., et al. 2015, ApJ, 800, 33 [Google Scholar]
  5. Archibald, R. F., Scholz, P., Kaspi, V. M., Tendulkar, S. P., & Beardmore, A. P. 2020, ApJ, 889, 160 [Google Scholar]
  6. Armaza, C., Reisenegger, A., & Valdivia, J. A. 2015, ApJ, 802, 121 [Google Scholar]
  7. Beloborodov, A. M. 2009, ApJ, 703, 1044 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beloborodov, A. M., & Thompson, C. 2007, ApJ, 657, 967 [Google Scholar]
  9. Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bransgrove, A., Beloborodov, A. M., & Levin, Y. 2023, ApJ, 958, L9 [NASA ADS] [CrossRef] [Google Scholar]
  11. Camilo, F., Cognard, I., Ransom, S. M., et al. 2007, ApJ, 663, 497 [Google Scholar]
  12. Camilo, F., Reynolds, J., Johnston, S., Halpern, J. P., & Ransom, S. M. 2008, ApJ, 679, 681 [Google Scholar]
  13. Cerutti, B., & Giacinti, G. 2021, A&A, 656, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Chen, A. Y., & Beloborodov, A. M. 2014, ApJ, 795, L22 [NASA ADS] [CrossRef] [Google Scholar]
  15. Çıkıntoğlu, S., Şaşmaz Muş, S., & Ekşi, K. Y. 2020, MNRAS, 496, 2183 [Google Scholar]
  16. Ciolfi, R., & Rezzolla, L. 2013, MNRAS, 435, L43 [NASA ADS] [CrossRef] [Google Scholar]
  17. Comisso, L., Lingam, M., Huang, Y. M., & Bhattacharjee, A. 2017, ApJ, 850, 142 [Google Scholar]
  18. Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351 [Google Scholar]
  19. Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. N. 2023, MNRAS, 527, L127 [CrossRef] [Google Scholar]
  20. Contopoulos, I., Ntotsikas, D., & Gourgouliatos, K. 2024, MNRAS, 527, L127 [Google Scholar]
  21. Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961 [NASA ADS] [CrossRef] [Google Scholar]
  22. Dall’Osso, S., & Stella, L. 2022, in Millisecond Magnetars, eds. S. Bhattacharyya, A. Papitto, & D. Bhattacharya, Astrophysics and Space Science Library, 465, 245 [Google Scholar]
  23. Dimitropoulos, I., Contopoulos, I., Mpisketzis, V., & Chaniadakis, E. 2024, MNRAS, 528, 3141 [Google Scholar]
  24. Fisher, R., Butterworth, E. M., Rajwade, K. M., et al. 2024, MNRAS, 528, 3833 [Google Scholar]
  25. Fujisawa, K., & Kisaka, S. 2014, MNRAS, 445, 2777 [Google Scholar]
  26. Glampedakis, K., & Lasky, P. D. 2016, MNRAS, 463, 2542 [Google Scholar]
  27. Glampedakis, K., Lander, S., & Andersson, N. 2014, MNRAS, 437, 2 [Google Scholar]
  28. Goodwin, S. P., Mestel, J., Mestel, L., & Wright, G. A. 2004, MNRAS, 349, 213 [Google Scholar]
  29. Gourgouliatos, K. N., & Esposito, P. 2018, in Strongly Magnetized Pulsars: Explosive Events and Evolution, eds. L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, & I. Vidaña, Astrophysics and Space Science Library, 457, 57 [Google Scholar]
  30. Gourgouliatos, K. N., & Hollerbach, R. 2016, MNRAS, 463, 3381 [Google Scholar]
  31. Gourgouliatos, K., & Lynden-Bell, D. 2008, MNRAS, 391, 268 [Google Scholar]
  32. Gourgouliatos, K. N., & Lynden-Bell, D. 2019, MNRAS, 482, 1942 [Google Scholar]
  33. Gourgouliatos, K. N., & Pons, J. A. 2019, Phys. Rev. Res., 1, 032049 [Google Scholar]
  34. Gourgouliatos, K. N., Cumming, A., Reisenegger, A., et al. 2013, MNRAS, 434, 2480 [Google Scholar]
  35. Gruzinov, A. 2005, Phys. Rev. Lett., 94, 021101 [Google Scholar]
  36. Haskell, B., Samuelsson, L., Glampedakis, K., & Andersson, N. 2008, MNRAS, 385, 531 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hu, C. P., & Ng, C. Y. 2019, Astron. Nachr., 340, 340 [Google Scholar]
  38. Hurley, K., Cline, T., Mazets, E., et al. 1999, Nature, 397, 41 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hurley, K., Boggs, S. E., Smith, D. M., et al. 2005, Nature, 434, 1098 [NASA ADS] [CrossRef] [Google Scholar]
  40. Kalapotharakos, C., & Contopoulos, I. 2009, A&A, 496, 495 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44 [NASA ADS] [CrossRef] [Google Scholar]
  42. Karageorgopoulos, V., Gourgouliatos, K. N., & Contopoulos, I. 2019, MNRAS, 487, 3333 [NASA ADS] [CrossRef] [Google Scholar]
  43. Karuppusamy, R., Stappers, B. W., & van Straten, W. 2010, A&A, 515, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Lander, S. K., & Jones, D. I. 2009, MNRAS, 395, 2162 [NASA ADS] [CrossRef] [Google Scholar]
  45. Lattimer, J. M., & Prakash, M. 2001, ApJ, 550, 426 [NASA ADS] [CrossRef] [Google Scholar]
  46. Levin, L., Lyne, A. G., Desvignes, G., et al. 2019, MNRAS, 488, 5251 [Google Scholar]
  47. Lower, M. E., Younes, G., Scholz, P., et al. 2023, ApJ, 945, 153 [Google Scholar]
  48. Lower, M. E., Karastergiou, A., Johnston, S., et al. 2025, MNRAS, 538, 3104 [Google Scholar]
  49. Lynden-Bell, D., & Boily, C. 1994, MNRAS, 267, 146 [Google Scholar]
  50. Lyne, A., Hobbs, G., Kramer, M., Stairs, I., & Stappers, B. 2010, Science, 329, 408 [NASA ADS] [CrossRef] [Google Scholar]
  51. Lyutikov, M. 2006, MNRAS, 367, 1594 [CrossRef] [Google Scholar]
  52. Mazets, E. P., Golentskii, S. V., Ilinskii, V. N., Aptekar, R. L., & Guryan, I. A. 1979, Nature, 282, 587 [NASA ADS] [CrossRef] [Google Scholar]
  53. Ntotsikas, D., Gourgouliatos, K., Contopoulos, I., & Lander, S. 2024, MNRAS, 527, 6691 [Google Scholar]
  54. Parfrey, K., Beloborodov, A. M., & Hui, L. 2012, ApJ, 754, L12 [Google Scholar]
  55. Parfrey, K., Beloborodov, A. M., & Hui, L. 2013, ApJ, 774, 92 [Google Scholar]
  56. Pavan, L., Turolla, R., Zane, S., & Nobili, L. 2009, MNRAS, 395, 753 [NASA ADS] [CrossRef] [Google Scholar]
  57. Philippov, A. A., & Spitkovsky, A. 2014, ApJ, 785, L33 [NASA ADS] [CrossRef] [Google Scholar]
  58. Philippov, A. A., Spitkovsky, A., & Cerutti, B. 2015, ApJS, 801, L19 [Google Scholar]
  59. Philippov, A., Timokhin, A., & Spitkovsky, A. 2020, Phys. Rev. Lett., 124, 245101 [Google Scholar]
  60. Pili, A. G., Bucciantini, N., & Del Zanna, L. 2015, MNRAS, 447, 2821 [NASA ADS] [CrossRef] [Google Scholar]
  61. Pintore, F., Bernardini, F., Mereghetti, S., et al. 2016, MNRAS, 458, 2088 [Google Scholar]
  62. Pons, J. A., Link, B., Miralles, J. A., & Geppert, U. 2007, Phys. Rev. Lett., 98, 071101 [Google Scholar]
  63. Pons, J. A., Viganò, D., & Geppert, U. 2012, A&A, 547, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Rankin, J. M., Wright, G. A. E., & Brown, A. M. 2013, MNRAS, 433, 445 [NASA ADS] [CrossRef] [Google Scholar]
  65. Reisenegger, A. 2009, A&A, 499, 557 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951 [Google Scholar]
  67. Scholz, P., Camilo, F., Sarkissian, J., et al. 2017, ApJ, 841, 126 [Google Scholar]
  68. Soglasnov, V. A., Popov, M. V., Bartel, N., et al. 2004, ApJ, 616, 439 [Google Scholar]
  69. Spitkovsky, A. 2006, ApJ, 648, L51 [Google Scholar]
  70. Stefanou, P., Pons, J. A., & Cerdá-Durán, P. 2022, MNRAS, 518, 6390 [Google Scholar]
  71. Suvorov, A. G., & Glampedakis, K. 2023, Phys. Rev. D, 108, 084006 [Google Scholar]
  72. Thompson, C., Lyutikov, M., & Kulkarni, S. R. 2002, ApJ, 574, 332 [NASA ADS] [CrossRef] [Google Scholar]
  73. Timokhin, A. N. 2006, MNRAS, 368, 1055 [Google Scholar]
  74. Tong, H. 2019, MNRAS, 489, 3769 [Google Scholar]
  75. Tsang, D., & Gourgouliatos, K. N. 2013, ApJ, 773, L17 [NASA ADS] [CrossRef] [Google Scholar]
  76. Urama, J. O., Link, B., & Weisberg, J. M. 2006, MNRAS, 370, L76 [NASA ADS] [Google Scholar]
  77. Uryū, K., Gourgoulhon, E., Markakis, C. M., et al. 2014, Phys. Rev. D, 90, 101501 [CrossRef] [Google Scholar]
  78. Uzdensky, D. A. 2003, ApJ, 598, 446 [NASA ADS] [CrossRef] [Google Scholar]
  79. Viganò, D., Pons, J. A., & Miralles, J. A. 2011, A&A, 533, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Woods, P. M., Kouveliotou, C., Göğüş, E., et al. 2002, ApJ, 576, 381 [Google Scholar]
  81. Younes, G., Güver, T., Kouveliotou, C., et al. 2020a, ApJ, 904, L21 [CrossRef] [Google Scholar]
  82. Younes, G., Baring, M. G., Kouveliotou, C., et al. 2020b, ApJ, 889, L27 [NASA ADS] [CrossRef] [Google Scholar]
  83. Younes, G., Lander, S. K., Baring, M. G., et al. 2025, ApJ, submitted [arXiv:2502.20079] [Google Scholar]

All Tables

Table 1.

Simulation results for the global twist model.

In the text
Table 2.

Simulation results for the internal twist model.

In the text

All Figures

thumbnail Fig. 1.

Various regions of the domain, labelled (I) to (V). The field lines are depicted in black, the last open field line of the magnetosphere in red, the stellar surface in green, and the light cylinder in dashed blue. The top-right inset provides a zoomed-in view of the star.
In the text
thumbnail Fig. 2.

Magnetic field of a relativistically rotating neutron star and its magnetosphere for the global twist model. The field lines are shown in black, the stellar surface in green, and the light cylinder with a dashed blue line at R = 1; the current sheet is also shown. The red line represents the first closed field line of the untwisted case (panel a) and is shown in subsequent panels for comparison. A poloidal current is applied in the area of the closed magnetic field lines with coefficients α = 0.0,  3.0,  6.0, and 8.0. The colour bar indicates the value of II′ across several parts of the star and the magnetosphere. The maximum value of II′ in the twisted field lines is saturated in colour to allow the depiction of its structure in the rest of the system. The top-right inset is shows a zoomed-in area of the star.
In the text
thumbnail Fig. 3.

Magnetic field of the internal twist model. Panels a and b correspond to α = 5.0 and 10.0, respectively. The various field lines and colour bars are shown as in Fig. 2.
In the text
thumbnail Fig. 4.

Magnetic flux function of the first closed magnetic field line as a function of α.
In the text
thumbnail Fig. 5.

Spin-down luminosity of the various solutions scaled to the spin-down power of the untwisted solution, as a function of α.
In the text
thumbnail Fig. 6.

Twist of the magnetospheric part of the field line for which Ψ = 1.1Ψ0, as a function of the parameter α.
In the text
thumbnail Fig. 7.

Polar cap’s semi-opening angle as a function of α.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.