Twisted magnetar magnetospheres

Magnetar magnetospheres are strongly twisted, and are able to power sudden energetic events through the rapid release of stored electromagnetic energy. In this paper, we investigate twisted relativistic force-free axisymmetric magnetospheres of rotating neutron stars. We obtain numerical solutions of such configurations using the method of simultaneous relaxation for the magnetic field inside and outside the light-cylinder. We introduce a toroidal magnetic field in the region of closed field-lines that is associated with a poloidal electric current distribution in that region, and explore various mathematical expressions for that distribution. We find that, by increasing the twist, a larger fraction of magnetic field-lines crosses the light-cylinder and opens up to infinity, thus increasing the size of the polar caps and enhancing the spin-down rate. We also find that, for moderately to strongly twisted magnetospheres, the region of closed field-lines ends at some distance inside the light-cylinder. We discuss the implications of these solutions on the variation of magnetar spin-down rates, moding and nulling of pulsars, the relation between the angular shear and the twist and the overall shape of the magnetosphere.


INTRODUCTION
A major part of pulsar study is the modeling of their magnetosphere, that is the neutron star atmosphere with its non-trivial distribution of magnetic field-lines, charge density, and electric currents.In a broader astrophysical context, early studies of magnetospheres focused initially on vacuum fields (Deutsch 1955).The work of Goldreich & Julian (1969) established the framework for a force-free magnetosphere of a relativistically rotating pulsar, initiating the effort for analytical and numerical solutions.The first solution of an aligned, axisymmetric and rotating neutron star with a smooth transition of the magnetic field-lines from the inner magnetosphere to the outer magnetosphere through the light-cylinder (Contopoulos et al. 1999) led to a series of studies of relativistic pulsar magnetospheres in axisymmetry (Uzdensky 2003;Goodwin et al. 2004;Gruzinov 2005;Timokhin 2006).This was followed by timedependent simulations of dipolar axisymmetric magnetospheres of neutron stars carried out within the framework of relativistic magnetohydrodynamics (Komissarov 2006), three-dimensional models describing the magnetosphere of oblique and orthogonal rotators (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009) and resistive force-free electrodynamics (FFE) (Kalapotharakos et al. 2012).
A key element of the above models is that field-lines crossing the light-cylinder are swept back, they are thus twisted.This ensures that the transition through the light-cylinder is smooth, thus solving one of the major hurdles of early modelling efforts (Mestel   ⋆ Email: d.ntotsikas@upnet.com† Email: kngourg@upatras.gr‡ Email: icontop@academyofathens.gr et al. 1979, 1985).Yet, these models do not include any twist along magnetic field-lines that do not cross the light cylinder.Twisted magnetospheres may be relevant to magnetars which are neutron stars with anomalous X-ray emission and Soft Gamma-Ray repeating events (Turolla et al. 2015;Kaspi & Beloborodov 2017).Their strong magnetic field powers these events (Duncan & Thompson 1992;Thompson & Duncan 1995, 1996;Turolla et al. 2015;Kaspi & Beloborodov 2017;Gourgouliatos & Esposito 2018;Esposito et al. 2021).These events release electromagnetic energy through bursts, outbursts and flares (Mazets et al. 1979;Hurley et al. 1999Hurley et al. , 2005;;Woods & Thompson 2004;Coti Zelati et al. 2018), and are often also accompanied by a variation of their 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).
As both the radiative profile and spin-down rate of neutron stars are determined by the magnetic field structure, variations of the magnetic field will affect them both.Moreover, the amount of energy stored within the magnetosphere which can be rapidly released in powering the aforementioned energetic events is altered.This can be achieved by twisting the magnetic field-lines, through shearing of their foot-points.Models of such twisted and sheared magnetic fields have been presented in the astrophysical literature (Aly 1994;Lynden-Bell & Boily 1994;Wolfson 1995;Gourgouliatos 2008) with particular focus on magnetar flares and outbursts (Thompson et al. 2002;Beloborodov & Thompson 2007;Pavan et al. 2009;Glampedakis et al. 2014;Akgün et al. 2016;Kojima 2017;Akgün et al. 2018;Carrasco et al. 2019).The basic idea of these models is that the magnetic field can accommodate a maximum amount of twist, beyond which it becomes unstable and releases its energy in an explosive manner, possibly through ideal or resistive instabilities (Lyutikov 2003(Lyutikov , 2006;;Beloborodov 2012;Tong 2019;Mahlmann et al. 2019).These mechanisms are typically studied for non-rotating magnetospheres.This is justified if most of the activity remains close to the stellar surface where the effects of rotation are minimal.However, as it has been shown in numerical modelling (Parfrey et al. 2012;Parfrey et al. 2013), the relativistic effects cannot be neglected if the twist is strong enough, since the poloidal current may occupy a large fraction of the magnetosphere that extends up to the lightcylinder.
A magnetar field can be twisted either through shearing of the stellar crust that is driven by internal magnetic stresses in the star (Lander & Gourgouliatos 2019;Gourgouliatos & Lander 2021), or through currents that flow from the crust to the magnetosphere (Fujisawa & Kisaka 2014;Mahlmann et al. 2019).The flow of poloidal currents within an axisymmetric magnetosphere, implies the presence of a toroidal field.If a closed magnetic field-line, i.e. a fieldline whose both footpoints are anchored on the crust, is differentially sheared a toroidal field will develop.On the contrary, an open fieldline, i.e. a field-line whose one footpoint is anchored on the crust and extends to infinity, cannot be twisted permanently from action originating at the star as the twist will relax beyond the light-cylinder, leading only to transient behaviour (Bransgrove et al. 2020).In this paper we impose a toroidal field that arises from the azimuthal shearing of the closed field-lines.Open field-lines contain only the poloidal current that is required for their smooth crossing of the light-cylinder, which cannot be externally imposed.
The plan of the paper is as follows.In section 2, we present the mathematical setup that describes a twisted magnetosphere under the force-free condition.We present the results of the numerical solution in section 3, and discuss them in section 4. Their astrophysical applications are presented in section 5, and we conclude with section 6.

MATHEMATICAL FORMULATION
In this section we provide the mathematical equations that describe axisymmetric ideal force-free magnetospheres.Assuming that the electromagnetic forces dominate gravitational and inertial ones, as well as pressure gradients, a system in equilibrium satisfies the condition of zero Lorentz force, namely where E is the electric field, B the magnetic field, ρ q the electric charge density, J the electric current density and c the speed of light.
In what follows we will work in cylindrical coordinates (R, ϕ, z), and will consider an axisymmetric configuration.Thus, the magnetic field can be expressed in terms of two scalar functions Ψ(R, z) and I(R, z): where, ∇ϕ ≡ φ/R.The function Ψ(R, z) is the poloidal flux and I(R, z) is proportional to the poloidal electric current contained within a circle of radius R whose centre lies on the axis of symmetry of the system at distance z from the horizontal plane.
The velocity is given by the expression: where Ω is the angular frequency.Assuming ideal MHD, we obtain the electric field by Ohm's law: Taking the divergence of the electric field we obtain the electric charge density.Then we substitute into equation ( 1), and by appropriately normalising lengths so that the light-cylinder radius is at R = R LC ≡ c/Ω = 1, we obtain the following equation: which is the axisymmetric pulsar equation derived by Scharlemann & Wagoner (1973).
In the standard solution of the pulsar equation (Contopoulos et al. 1999), the magnetosphere comprises of two parts: a region where both ends of the field-lines are anchored onto the surface of the star (the region of closed field-lines), where it is postulated that I(Ψ) = 0 and a region containing the field-lines that cross the light-cylinder and open up to infinity (the region of open field-lines), where in general I(Ψ) 0. The functional form of I(Ψ) in that region is prescribed by the requirement that field-lines cross the light-cylinder smoothly.
The condition that I(Ψ) = 0 for the closed field-lines is based on the assumption that the field of the pulsar does not have any intrinsic twist.While this assumption corresponds to a magnetic field of minimal complexity, the idea of twisted magnetic fields is commonplace in magnetar models.Solutions with non-zero I(Ψ) in the closed field-lines could be relevant and applicable to strongly magnetised systems.Thus, we study systematically systems with nonzero I(Ψ) within the closed region.
We integrate numerically equation ( 5) within an orthogonal domain of the northern half-space of the system R × z ∈ [0, R max ] × [0, z max ] excluding the star, which is assumed to have a radius of r NS , thus R 2 + z 2 > r 2 NS .In that domain, we apply the standard boundary conditions of the axisymmetric pulsar equation.For R = 0, that is the z-axis, the following relation holds: while on the surface of the star r NS = √ R 2 + z 2 , the magnetic flux function is that of a dipole: In the numerical setup, we place the inner boundary at r NS = 0.1R LC , which, in physical units, implies a rotation angular velocity Ω = 3000rad s −1 and a spin frequency of 500Hz.Because of North-South symmetry, the vertical derivative of the magnetic flux at the equator for the closed field-lines is equal to zero, namely Here, R Y is the radial extent of the region of closed field-lines, namely the cylindrical radius where the last closed field-line crosses the equator.That position is often called the Y-point.Beyond the last closed field-line, the system forms an equatorial current sheet, thus the boundary condition there is given by the following expression: In most studies, R Y is typically the light-cylinder radius, yet it has been recently argued that it is likely to be located at around 0.9 of the light-cylinder for the untwisted magnetosphere (Contopoulos et al. 2023).In the twisted case, enforcing the inner edge of the current sheet at the light-cylinder is too restrictive for some of the cases studied and has to be relaxed to allow for physically acceptable solutions.
As the pulsar equation has a critical point at the light-cylinder, we demand that the magnetic field is non-singular and continuous there, thus the magnetic field-lines cross the light-cylinder smoothly.This leads to Dirichlet and Neumann boundary conditions: where a prime denotes differentiation with respect to Ψ.The above relations allow us to determine the functional form of I(Ψ) for Ψ < Ψ 0 , where Ψ 0 is the flux function corresponding to the last open magnetic field-line.In the region of the closed magnetic field-lines (Ψ > Ψ 0 ) the electric current function is set to the following form: where α, n are parameters that determine the details of the electric current distribution.This permits the inclusion of twist in the closed magnetic field-lines.
In principle, one could have chosen a functional form for I(Ψ) where the toroidal field does not vanish at the first closed field line located at Ψ = Ψ 0 .The reason we opted for I(Ψ 0 ) = 0 is the following: In the standard untwisted pulsar magnetosphere, a discontinuity on the toroidal field appears across the separatrix between the open and closed field-lines.This discontinuity is accompanied by a poloidal current sheet, through which most of the current that flows through the polar cap and crosses the light-cylinder returns to the star.Therefore, most of the poloidal electric current that is contained at the open field-line region closes along the separatrix between the open and closed field-lines.As the toroidal magnetic field is equal to zero there, the modulus of the discontinuity of the poloidal current across the separatrix is equal in the north and the south hemisphere and the system maintains a symmetry across the equator.
Once we include a toroidal field in the closed-field line region, this could, in principle, interact asymmetrically with the separatrix.More precisely, while the toroidal field of the open field-lines just outside the separatrix has opposite directions in the northern and southern hemisphere, the toroidal field that we impose in the closed field-lines has the same polarity.This north-south asymmetry of the toroidal field in the closed-line region would, in general, cause a north-south asymmetry in the poloidal field configuration if the toroidal field just inside the separatrix is non-zero.The discontinuity of the toroidal field determines the intensity of the poloidal current sheet flowing along the separatrix.Thus, for the hemisphere where the toroidal field has the same direction across the separatrix, the discontinuity will be milder, whereas in the other one it will be steeper.Consequently, the poloidal return current sheet will be different in the northern and southern separatrices.To avoid such complications, we set the toroidal field at the edge of the closed field-line region (Ψ = Ψ 0 ) equal to zero, as it is evident from equation ( 12).This, ensures that the north-south symmetry of the poloidal magnetic field configuration is not broken.Otherwise, the toroidal field in the closed-line region is not north-south (anti)symmetric.
At the external boundary of the domain (z = z max and R = R max ) we implement the split monopole boundary conditions, as follows: Thus, the open magnetic lines become radial.These conditions hold at an infinite distance from the surface of the star.However, we apply them at the outer boundary of the integration domain.We have found that they have minimal impact on our calculation, when we increased the integration domain by a factor of two, the solution changed by less than 1%.

Magnetic field structure
We have integrated equation ( 5) for various choices of the parameters n and α appearing in equation ( 12) that determine the electric current flowing along the closed magnetic field-lines.We have used a typical (R, z) resolution of 400 × 1200, and we have also doubled the resolution in several cases to ensure the convergence of our calculations.We have solved the above equation applying the simultaneous relaxation method in an appropriately modified version of the numerical code used in Gourgouliatos & Lynden-Bell (2019).In the open field-line region, the value of I(Ψ) is determined by enforcing the conditions described in equations ( 10) and ( 11).At the closed field-lines we input the prescribed function for the electric current, given by equation ( 12), and we apply the elliptic solver.This process is repeated until the system converges to the solution of the system.The current flowing in the region of the open field-lines should formally return through a current sheet on the equatorial plane and the separatrix between the open and closed field-lines.As it is not possible to handle numerically an infinitesimal return current sheet in a finite difference code, we approximate it using a Gaussian function centred at Ψ = 0.975Ψ 0 and a width of 0.01Ψ 0 .This amplitude of the Gaussian is chosen so that the entire current closes through the star.Thus, in this narrow layer we do not apply the conditions of equations ( 10) and (11) but the above predetermined expression.
First, we have tested our numerical code for a system without any poloidal current at the closed magnetic field-line region (α = 0) to ensure it converges to the standard pulsar solution, figure 1, panel a.We find that the last open field-line occurs at Ψ 0 = 1.23 in agreement with the results of Timokhin (2006).Next, we have explored a wide variety of combinations of the parameters α and n shown in Tables 1-3.
The linear current relationship corresponds to n = 1, and we have solved for α = 0.1 to α = 10.In figure 1 we show solutions for the parameters starting from α = 0, the untwisted case, to α = 5.0, panels a to d.As α increases closed magnetic lines of the untwisted system open up gradually, which is reflected on the increase of Ψ 0 , that is found through the solution of the pulsar equation.Increasing α above 2.5, while enforcing the innermost point of the current sheet to remain at R Y = 1 leads to magnetospheric structures that contain magnetic field-lines that are disconnected from the star.While such solutions satisfy the mathematical form of the equations and the boundary conditions, they are physically unacceptable, as they include field-lines emerging from the equator without being linked to the star.To prevent such unphysical solutions we move the innermost point of the current sheet within the light-cylinder and closer to the surface of the star.This way, the current sheet starts earlier and all magnetic field-lines remain connected to the star.While we set manually the position of R Y , this is not done arbitrarily, but in a way to impose the minimum displacement of R Y to avoid the appearance of disconnected magnetic field lines.For α = 3.0, the current sheet starts at R Y = 0.95.The next shift of the current sheet occurs for α = 4.0 where the current sheet inner edge must be shifted to R Y = 0.9.We have used increments of 0.05 for the displacement of R Y , in R LC units.Should we have allowed for a finer change, the movement of R Y would be continuous with α.In table 1 we can see the increase in the value of the magnetic flux function at the last open magnetic  field-line with respect to the coefficient α and the displacement of the current sheet.Such a behavior is expected, as moving R Y closer to the star leads to larger values of the flux function.This, in combination with the increase of α leads to a nonlinear increase of the magnetic flux function in the corresponding last open magnetic line.The same behaviour is found for the nonlinear dependence of the current distribution to the magnetic flux function.In particular, we have studied the magnetosphere for n = 1.5 and α ranging from 0.1 to 10, with some solutions shown in Figure 2. We can obtain solutions with the current sheet inner edge being at R Y = 1, while α ⩽ 1.5.Once α reaches 2.0 a drastic shift of the inner edge of the current sheet is required to R Y = 0.85.Here also, the red line denotes the magnetic field-line for which Ψ 0 = 1.23.The increase of the coefficient α is followed by a stronger shift of the current sheet towards the star reaching R Y = 0.40 for α = 5.0.Similar to the linear case, the magnetic energy of the region of the closed magnetic lines increases initially with α.The change in the position of the current sheet also leads to a change in the volume of the occupied by the closed magnetic field-lines and eventually a decrease of the magnetic energy.
Qualitatively similar results can be drawn studying the non-linear current distribution for n = 2.0 and α ranging from 0.1 to 10.In this case, the inner edge of the current sheet has to be displaced even more drastically compared to the other two cases.Figure 3 shows the magnetic field structure for a poloidal current distribution with n = 2.0 for various values of α.The value of α = 1.0 corresponds to the maximum achievable value of the current function with n = 2.0, for which no displacement of the current sheet from the position R Y = 1 is required.For the highest value α = 10 we have studied, we found that R Y = 0.13 and the system is very close to the structure of a split monopole.

Spin-down luminosity and twist of the magnetosphere
Introducing a poloidal current in the closed field-line region, and increasing its intensity, results in more magnetic field-lines becoming open, Figure 4.As the spin-down rate is related to the torque exerted by the open magnetic field-lines, such an increase leads to a more efficient loss of the rotational energy of the star.At a given rotation rate, the spin-down luminosity L scales approximately as follows: We note that a more accurate approximation of energy losses can be obtained from the relation (Timokhin 2006): (the total spin-down luminosity is equal to the electromagnetic energy loss from both poles of the star).We have found that the difference of the results found by the two formulae are minimal and in the results we quote the ones that scale as Ψ 2 0 .From this we are able to obtain an estimate of the energy loss of the star under the influence of the increasing current in the region of the closed magnetic lines.Increasing the value of the current in the region of the closed magnetic lines is followed by an increase spin-down rate which is more effective for higher n.It is interesting to see how the star loses part of its rotational energy in the case where the region of closed magnetic lines in its magnetosphere is driven by a non-linear current distribution, with n = 2.0.In this case, we can observe a significantly steeper increase in the energy loss as the value of α increases compared to the previous two cases.It is also possible to discern sudden increases in the spin down power which also coincide with the respective displacement of the current sheet towards the inner magnetosphere of the star.If this displacement is significant, equation ( 15) may substantially underestimate the true spin-down rate; see section 4.3.
The equation of a field-line can be found by integrating the following relation: where B R , B z , B ϕ are the components of the magnetic field in cylindrical coordinates.We can evaluate the twist of a field-line, through integration of the equation involving the R and the ϕ components of the magnetic field.To estimate the maximum possible twist that the closed magnetic field-lines of the magnetosphere can undergo, the calculations were performed within the region of the closed magnetic field-lines at a position corresponding to a magnetic line having a magnetic flux Ψ = 1.1Ψ 0 .We performed the integration on various field-lines within the closed region and we have found for the particular choices of magnetic configuration the above expression provides the maximum twist.Using the first equality of equation ( 17) and integrating along a constant Ψ we obtain the following relation: We integrate from a point on the surface of the star in the northern hemisphere, where the field line emanates, to the equator, thus the twists quoted in the tables correspond to half of the field-line's twist.We plot the results in Figure 5.We notice that for the system where the current sheet is closest to the integration domain inner edge and α has the highest value quoted, (n = 1 α = 10; n = 1.5, α = 10 and n = 2, α = 10), bf the twist does not increase indefinitely but rather tends to approach a maximum value ∼ π/2 which is consistent with previous solutions (Lynden-Bell & Boily 1994; Gourgouliatos & Lynden-Bell 2008;Pavan et al. 2009;Akgün et al. 2018).Thus, after a finite amount of twist, the system can formally reach a state of a split monopole.In our calculation we have chosen to explore the characteristics in terms of the parameter α going up to a maximum value of α = 10.Thus while two simulations have the same value of α, different n corresponds to different amounts of twist, with the increase in the twist being more drastic for larger n that also require pushing the Y point closer to the star.
To calculate the opening angle of the polar cap, the boundary condition describing the magnetic flux function on the surface of the star was used, equation ( 7), which can be written using the polar angle as follows: and consequently, the angle of the polar cap will be given by the relation: In Figure 6 we plot the polar cap angle θ pc versus α.We verify that by increasing the twist the polar cap expands.A further conclusion we can draw is that the position of the current sheet plays a crucial role in shaping the polar cap, which is consistent with its impact on the magnetosphere in general.This claim can become more evident by studying the polar cap in the case where the magnetosphere is described by the non-linear current distribution α(Ψ − Ψ 0 ) 2.0 .In this case, in addition to the very large values obtained for the current      function, compared to the other two cases, the current sheet is subject to the strongest currents compared to the other cases.In these magnetospheres we can also observe the largest increases of the polar cap opening angle.

DISCUSSION
The inclusion of a twisted region in the relativistic pulsar magnetosphere has significant qualitative implications for its structure and main characteristics.We have identified the following changes: displacement of the inner edge of the equatorial current sheet; enhancement of the spin-down rate related to the twist and an increase to the size of the polar cap.

Location of the current sheet
An essential part of the axisymmetric pulsar magnetosphere is the equatorial current sheet.The field-lines that cross the light-cylinder cannot cross the equatorial plane, as this would imply that the particles that are frozen onto them should move faster than the speed of light.Therefore, these field-lines form two discontinuous streams separated by the equator.While the above requirement enforces the presence of a current sheet beyond the light-cylinder, there is no physical reason to prevent the inner edge of the current sheet to be closer to the star and within the light-cylinder.Such solutions are mathematically and physically viable.In most pulsar studies, the location of the current sheet inner edge is set to the light-cylinder, where a Y−point appears.Despite that, it has been shown that a family of smooth solutions can be obtained for current sheets starting within the light-cylinder (Timokhin 2006).In these solutions, there is no physical requirement to demand that the current sheet starts within the light-cylinder.However, as a twist is introduced in the magnetic field-lines, the displacement of the current sheet within the light-cylinder is no longer just an alternative solution, but a necessity in order to obtain a physically acceptable configuration.This is because in the untwisted pulsar magnetosphere, there is only one characteristic length-scale, the light-cylinder.In the twisted magnetosphere however, an additional length scale is introduced, which is inversely proportional to α, that determines the ratio of the magnetic field to the electric current flowing along a given magnetic field-line.This becomes evident from the solution of an axisymmetric forcefree field in spherical geometry, the so-called spheromak solution (Chandrasekhar & Kendall 1957), that has a characteristic radius at which the field-lines close.We further note that such constraints appear in relativistic generalisations of force-free solutions (Gourgouliatos & Vlahakis 2010; Barkov et al. 2022).Thus, in the solutions studied in this paper, we see the interplay between these two lengthscales: the light-cylinder and the force-free twisted field.As long as the twist is small, and consequently α −1 large, the force-free length scale is sufficiently long, the location where the twisted field-lines must formally close lies beyond the light-cylinder.Therefore, the inner edge of the current sheet is not affected.As α increases there will be a point where the length scale of the twisted field will be- come comparable to the light-cylinder.Once this stage is reached, we obtain solutions where the current sheet inner edge is located inside the light-cylinder.We further note that a similar effect has been noted in studies where additional length scales are introduced, such as the the Alfvèn radius (Metzger et al. 2011;Lander & Jones 2020).Similarly, dynamical studies of magnetospheres lead to transient configurations where the current sheet starts well within the light-cylinder (Barkov et al. 2022;Sharma et al. 2023).Moreover, in solutions obtained through force-free electrodynamics (Spitkovsky 2006), the magnetosphere very quickly opens up at R = 0.6R LC and it took it more than 20 rotations to gradually reach the light cylinder.

Polar cap
The current density throughout the magnetosphere of a pulsar is fuelled by the creation of positron-electron pairs, possibly in the polar cap region, filling the magnetosphere with charged particles (Ruderman & Sutherland 1975;Scharlemann et al. 1978;Muslimov & Tsygan 1992).As the region of the polar cap has been associated with radiation mechanisms and particle acceleration its physical size may be relevant to pulse shape.Enhancing the twist of the magnetosphere, leads to an expansion of the polar cap region.This is related to the fact that more field-lines become open.Since the opening angle of the polar cap depends on the region where the last open field-line crosses the star, having a larger fraction of open field-lines implies a larger polar cap.We notice significant increases of the polar cap once it is required to displace the current sheet closer to the surface of the star in magnetospheres with a current sheet near the surface of the star.

Twist and spin-down rate
In general, twisted magnetospheres are expected to have more efficient spin-down rates due to the fact that field-lines tend to become radial and approach the structure of a split monopole (Michel 1982;Thompson et al. 2002;Lyutikov 2006).An additional factor contributing to more efficient spin-down of twisted magnetospheres is the displacement of the current sheet closer to the star.In figure 7 we plot the spin-down rate as a function of the displacement of the current sheet.It is evident that it increases as the current sheet moves inwards in the magnetosphere.We note however, that while we may have the same displacement of the Y-point, the spin-down rate does not become identical in these cases, as the details depend on the functional form of I = I(Ψ) in the closed field-line region, leading to some deviations.
In fact, our results may underestimate the spindown rate for these cases, since we assume the luminosity scales with Ψ 2 0 only.In practice, there may be an additional contribution from the qualitative change to the magnetospheric structure: because both the polar cap region is larger and R Y closer to the star, there is a larger volume of the exterior where electromagnetic waves carry away angular momentum from the star.It is interesting to compare our results with the case of adding a neutrino-driven wind, which also results in a larger volume of open field-lines.In the semi-analytical analysis of Lander & Jones (2020) (specifically, combining their equations 2,6, and 12), it may be shown that the spindown luminosity has an additional scaling with (R L /R Y ) 2 due to the changes in magnetospheric structure.A similar contribution is likely to exist in our case, and would result in lines of steeper gradient in fig. 7.

APPLICATIONS
While the discussion has been primarily numerical and focusing on the generic features of twisted magnetospheres we can still consider some applications to specific astrophysical systems.The effects described in this work may be related to the switching of neutron stars between pulsar and magnetar states, as this has been observed in PSR J1119-6127 (Majid et al. 2017).This source has been exhibiting the following behaviour: Two glitches occurred in 2004 and 2007.The glitch events were followed by intermittent and RRATtype emission; the pulsar recovered from the glitch in an interval between 20 and 210 days.The spin-down rate recovered to a value smaller than the pre-glitch one (Weltevrede et al. 2011).Another event occurred in 2016 with a glitch event accompanied by radiative changes (Dai et al. 2018).While our model is a series of static axisymmetric equilibrium states and cannot capture this type of event in its totality, some of the characteristics can be described in this framework.A glitch can provide a sudden twist in the magnetosphere, part of which may be channeled to the open field-lines, but a fraction of the twist may be trapped in the closed field-lines (Bransgrove et al. 2020).Following that, a twisted magnetosphere will lead to a higher spin-down rate compared to the pre-glitch state.Should the magnetosphere have some resistivity, the spin-down rate will gradually drop, with a minimum value corresponding to the completely untwisted one.Furthermore, the process of twisting has three effects: an increase in the current flowing within the closed field-line region, the expansion of the polar cap and the displacement of the current sheet within the light-cylinder.The first one may increase the amount of plasma in the closed field-lines and shadow radio pulses, which is actually the case, while particles accelerated in this region may bombard the crust and lead to X-ray emission (Beloborodov 2009).The displacement of the current sheet and the expansion of the polar cap affects the interface between the closed and open fieldlines, which is related to the gaps where the creation of pairs is likely.Thus, it is likely that the presence of a twisted region of the fieldlines can lead to magnetar behaviour to an otherwise normal pulsar, even if it is not highly-magnetised.
We further remark that some pulsars exhibit moding and nulling behaviour, and changes in the subpulses, (Deich et al. 1986;Esamdin et al. 2005;Bhattacharyya et al. 2010;Rankin et al. 2013).These systems include the intermittent pulsars: B1931+24 (Rea et al. 2008), J1832+0029 (Lorimer et al. 2012), J1107-5907 (Meyers et al. 2018), J1832+0029 and J1841-0500 (Wang et al. 2020).These sources switch between states with different pulse shape, or they even completely turn off, while at the same time their spin-down frequency changes.These effects have been explained in the context of changes in the size of the corotating region of the magnetosphere and transitions between metastable configurations (Timokhin 2010).While we can construct solutions with different spindowns, we cannot at this point offer a clear mechanism that would show how the system can switch between these two configurations.However, qualitatively, we may expect the generation of such states: in a twisted magnetosphere the location of the current sheet is displaced, strongly affecting the closed field-line area and consequently the spin-down rate.Therefore, the system may transition between states with different spin-down rates.In particular, the intermittent pulsar B1931+24 alternates between an on and an off state, when the radio emission switches on and off respectively, with the on state lasting about 5-10 days, and the off state 25-35.The spin down rate alternates between: von = (−16.30± 0.04) × 10 −15 s −2 and vof f = (−10.80± 0.02) × 10 −15 s −2 .Considering the runs with n = 1, we notice that the states with α = 2 and α = 3.5 have spin-down powers whose ratio is 2.13/1.32= 1.61.This ratio is close to the ratio of the spin-down powers of the two states, with the highly twisted one being the on phase and the less twisted one being the off phase.

CONCLUSIONS
The main focus of this paper is the structure of a twisted axisymmetric magnetosphere of a relativistically rotating neutron star, by numerically solving the relativistic force-free equation.The increase of the poloidal current in the closed magnetic lines, or equivalently the introduction of a twist in the closed field-lines, affects the global structure of the magnetosphere and it is found to be qualitatively different from the standard (untwisted) pulsar magnetosphere.We find that part of the magnetic field-lines becomes open.Moreover, we find that it is not possible to have a simultaneous increase of the polar current and a physically consistent structure of the magnetosphere without a displacement of the current sheet towards the surface of the star.We conclude that the shift of the current sheet is unavoidable if the current in the magnetosphere increases above a certain value.During the movement of the current sheet towards the surface of the star, the volume of closed field-lines decreases.Thus, following an initial increase in the energy of the closed field-line, it eventually drops.Furthermore, the maximum value of the twist is approximately ∆ϕ max = π/2, which is in agreement with previous studies, thus after finite twist, the magnetosphere will open up completely adopting the structure of a split monopole.
It is possible that variations in the timing behavior of magnetars before and after outbursts can be associated with the magnetospheric twist.We conclude that by twisting the closed magnetic field-lines the field becomes more radial, a larger fraction of magnetic flux crosses the light-cylinder and the spin-down rate becomes higher.As the twist cannot be indefinitely increased we anticipate that after some twisting episode the field will return to the untwisted state either through a sudden event, i.e. reconnection of the open field-lines, or through gradual dissipation of the current.
We note that these results can be further explored in a 3-D study since an axisymmetric system cannot simulate a pulsating source.Even in axial symmetry the inclusion of quadrupolar, or in general multipolar fields, breaks the north-south symmetry.The lack of a north-south symmetry leads to a current sheet that does not coincide with the equator, and would require a drastically different treatment from the finite difference code we are employing here.
While these results are drawn from an idealised case of an axisymmetric rotator, they can provide useful conclusions related to the behaviour of magnetars and neutron stars.For instance, changes in the twist affects the shape of the closed region, the spin down rate and the polar cap.More interestingly, we notice that a pulsar can switch between twisted and untwisted states, without drastically changing the energy content of the closed region.This could be possibly related to intermittent sources, and sources switching from pulsar to magnetar behaviour and vice versa.

Figure 1 .
Figure 1.Pulsar magnetospheres with a linear poloidal current distribution for various values of α, ranging from α = 0 to α = 5.The contours of constant Ψ (poloidal field-lines) are shown in black, and the poloidal current I is shown in colour.The red line corresponds to the last open field-line of the standard pulsar solution Ψ 0 = 1.23.

Figure 2 .
Figure 2. Magnetospheres including a non linear poloidal current distribution (n = 1.5) for values of α ranging from α = 2.0 to α = 5.The contours of constant Ψ (poloidal field-lines) are shown in black and in color I.The red line corresponds to the last open field-line of the standard pulsar solution Ψ 0 = 1.23.

Figure 3 .
Figure 3. Magnetospheres under the influence of a non linear polar current distribution (n = 2.0) for values of α ranging from α = 1.5 to α = 4.0.The contours of constant Ψ (poloidal field-lines) are shown in black and in color I.The red line corresponds to the last open field-line of the standard pulsar solution Ψ 0 = 1.23.

Figure 4 .
Figure 4. Magnetic flux function of the last closed magnetic line for a variety of current distributions.

Figure 5 .
Figure5.An approximation of how the magnetosphere in the region of the closed magnetic field-line twists, as the value of the current function increases in the closed magnetic lines of its magnetosphere and as the current sheet shifts towards its surface.

Figure 6 .
Figure6.An approximation of how the polar cap increases , as the value of the current function increases in the closed magnetic lines of its magnetosphere and as the current sheet shifts towards its surface.

Figure 7 .
Figure7.The spin-down power of a twisted magnetosphere scaled to the spin-down power of an untwisted one against α for the various values of n.

Table 2 .
Simulation results for a range of values of the parameter α and n = 1.5.

Table 3 .
Simulation results for a range of values of the parameter α and n = 2.0.