Three-dimensional GRMHD Simulations of Neutron Star Jets

Neutron stars and black holes in X-ray binaries are observed to host strong collimated jets in the hard spectral state. Numerical simulations can act as a valuable tool in understanding the mechanisms behind jet formation and its properties. Although there have been significant efforts in understanding black hole jets from general relativistic magnetohydrodynamic (GRMHD) simulations in recent years, neutron star jets still remain poorly explored. We present the results from three-dimensional GRMHD simulations of accreting neutron stars with oblique magnetospheres for the very first time. The jets in our simulations are produced due to the anchored magnetic field of the rotating star in analogy with the Blandford–Znajek process. We find that for accreting stars, the star–disk magnetic field interaction plays a significant role, and as a result, the jet power becomes directly proportional to Φ2jet , where Φjet is the open flux in the jet. The jet power decreases with increasing stellar magnetic inclination, and finally, for an orthogonal magnetosphere, it reduces by a factor of ≃2.95 compared to the aligned case. We also find that in the strong propeller regime, with a highly oblique magnetosphere, the disk-induced collimation of the open stellar flux preserves parts of the striped wind, resulting in a striped jet.


INTRODUCTION
One of the key features of disk-fed accreting systems is the outflow of material either in the form of collimated jets or conical winds.This occurs in systems hosting various central objects ranging from neutron stars and black-holes in X-ray binaries (XRBs), black holes in active galactic nuclei (AGN) to young stellar objects (YSOs), and white dwarfs.Several mechanisms have been suggested to explain the formation of outflows from the accretion disk (e.g.Begelman et al. 1983;Icke 1980;Vitello & Shlosman 1988) which can all contribute to a larger or lesser extent for a given system (Proga 2007;Pudritz et al. 2007).The mechanisms of relativistic jet formation from black holes and neutron stars are most commonly ascribed to the theoretical models suggested by Blandford & Znajek (1977) and Blandford & Payne (1982).The Blandford-Payne (BP) process suggests that the accretion disk threaded by large poloidal magnetic fields can power astrophysical jets.The Blandford-Znajek (BZ) mechanism on the other hand relies on the frame-dragging effect of a Corresponding author: Pushpita Das p.das2@uva.nlrotating black hole which, if threaded by poloidal magnetic fields, launches a Poynting flux driven outflow along the field lines.For compact stars, a third option exists in which the field lines are rooted in the rotating star, and a jet similar to the BZ process is launched from the magnetically dominated open field lines near the polar cap (Parfrey et al. 2016;Parfrey & Tchekhovskoy 2017;Das et al. 2022).
Jets in X-ray binaries are most commonly studied in the radio band.The radio emission originates close to the central object and is typically associated with non-thermal emission from electrons in the jet (Falcke & Markoff 2000).The evolution of the jet is strongly coupled to the nature of the accretion flow in the vicinity of the compact object.In black-hole X-ray binaries (BHXRBs) a steady radio jet is often observed in the hard X-ray state, during transition to the soft state the jet grows transitional and more powerful, and is finally suppressed during the soft state (Fender et al. 2004).Similar to BHXRBs a strong radio jet is also observed during hard states in neutron star X-ray binaries (NSXRBs).However, as the system evolves towards soft X-ray states, contrary to blackholes, some show signs of jet suppression (Miller-Jones et al. 2010) and some do not (Migliari et al. 2004).Along with the differences in jet behaviour with spectral states, for similar X-ray luminosity, neutron star jets are also observed to be radio-fainter by a factor of ∼ 10 ( Gallo et al. 2018;van den Eijnden et al. 2021) compared to their black hole counterparts and it has been argued that this difference cannot be explained by the additional X-ray emission from the surface (Fender & Kuulkers 2001;Migliari & Fender 2006;Tudor et al. 2017).
There have been substantial studies on the radio vs. Xray correlation in X-ray binaries (Hannikainen et al. 1998;Merloni et al. 2003;Fender et al. 2009;Gallo et al. 2014;Carotenuto et al. 2021;Gallo et al. 2018;van den Eijnden et al. 2021).Besides the systematic radio-quietness of neutron star systems, these studies have shown that NSXRBs are more "varied" than BHXRBs.For example, individual neutron star systems that are otherwise very similar (in terms of Ṁ , mass, spin, and magnetic field strength) show over an order of magnitude difference in radio power (e.g.Russell et al. 2018;Tetarenko et al. 2018).Furthermore, several sources display large uncorrelated jumps in either radio-or X-ray power on timescales of days (e.g.Tudor et al. 2017;van den Eijnden et al. 2020).The differences in black hole and neutron star jet properties are often attributed to a different disk-jet coupling mechanism, however, no consensus has been reached yet.
Numerical simulations of star-disk systems can be used to study jet formation and disk-jet coupling in NSXRBs.Stellar outflows and star-disk magnetic interaction have been extensively studied with non-relativistic MHD simulations (Matt & Pudritz 2005;Zanni & Ferreira 2013;Lovelace et al. 2014;Lii et al. 2014;Romanova et al. 2018;Takasao et al. 2022) and relativistic simulations (Parfrey et al. 2016).Axisymmetric GRMHD simulations of accreting neutron stars have been recently performed by Parfrey & Tchekhovskoy (2017); Das et al. (2022);C ¸ıkıntoglu et al. (2022).Due to their twodimensional nature, the previous simulations were restricted to aligned configurations of the (quadru-) dipolar stellar field with the disk magnetic initial field configuration.In this paper, we show results from the first 3D GRMHD simulations of accreting neutron stars and investigate the role of stellar magnetic inclination in jet formation.

NUMERICAL SETUP
We solve the following ideal MHD equations in a general relativistic framework using The Black Hole Accretion Code (BHAC) (Porth et al. 2017;Olivares et al. 2019), (1) Here, T µν is the energy-momentum tensor of an ideal magneto-fluid; ⋆ F µν is the dual of the Faraday tensor F µν and ρ, u µ , are rest-mass density and fluid 4-velocity respectively.We perform the simulations in Cartesian Schwarzschild coordinates corotating at the stellar angular frequency (Ω star ).The metric components for our tailormade coordinate system are given in Appendix A. We initialize our simulation domain with a standard Fishbone-Moncrief torus (Fishbone & Moncrief 1976), with the inner edge and density maximum located at r in = 45r g (r g = GM/c 2 ) and r max = 65r g respectively.The torus is weakly magnetized with poloidal loops defined as A ϕ ∝ max(ρ/ρ max − 0.1, 0) such that 2p max /b 2 max ≈ 110 (where the subscript max indicates the individual maxima of the quantities as customary in black hole torus setups, Porth et al. 2019).Outside the torus, the magnetosphere is initialized with a stationary dipolar vector potential in Schwarzschild space-time following Wasserman & Shapiro (1983).
Inside the magnetosphere (r < r lc = c/Ω), the initial pressure and density are set such that magnetization (σ t = b 2 /ρ) and plasma beta (β t = 2p/b 2 ) are 70 and 0.014 respectively.Both σ and β transition smoothly to follow an r −6 profile outside the light-cylinder radius.We perturb the pressure with 8% white noise to excite the magnetorotational instability (MRI) inside the torus.The development of MRI results in angular momentum transport in the disk Balbus & Hawley (1991) and leads to accretion onto the neutron star.
In order to maintain force-free like conditions in the magnetosphere, we modify the numerical solution to drive pressure, density, and v ∥ to target values (ρ t = b 2 /σ t , p t = β t b 2 /2 and v ∥,t = 0, see Das et al. (2022) for more details).In order to differentiate the accreting region from the magnetosphere, we introduce a tracer-fluid (τ ).The accretion disk is initialized with τ = 1 and the rest of the domain is initialized with as τ = 0.The tracer is passively advected with the flow by solving the extra equation ∇ µ (τ ρu µ ) = 0. Inside the magnetosphere (regions with τ = 0), we drive the solution to the target values as mentioned above, for 0 < τ < 1 we mix the force-free and MHD solutions by interpolating between the two solutions.Finally, for τ = 1, the solution follows that of the unmodified ideal GRMHD equations.We refer the reader to Parfrey & Tchekhovskoy (2017) for further discussion of this approach to couple the force-free and ideal MHD regimes.For numerical stability, along with driving the solutions to target values, we also maintain global floor and ceiling values for β and σ such that at any point in time β > 0.1β t and σ < 10σ t .
Since the star is part of the computational grid, the stellar boundary conditions require special attention.Inside the stellar radius, we set pressure and density to initial values such that β = 0.014 and σ = 70, and the velocity components (v x , v y , v z ) are overwritten such that the star maintains a solid body rotation with an angular frequency Ω. Thereby, matter that reaches the stellar surface is effectively removed from the simulation, emulating the outflow boundary conditions used in spherical setups (Parfrey & Tchekhovskoy 2017;Das et al. 2022).We update the magnetic field using the constrained transport algorithm of BHAC (Olivares et al. 2019).To fix the stellar magnetic field to the initial oblique dipole, within r/R ⋆ < 0.95 we enforce that the co-rotating electric fields smoothly transition to zero using a sigmoid function.This avoids any kinks of the magnetic field-lines at the stellar boundary.The induction equation is solved in the entire domain, ensuring that no magnetic monopoles are created throughout the evolution.We tested that the method of flux freezing preserves the initial configuration with sufficient accuracy.In particular, the relative difference in the magnetic flux (at r = R ⋆ ) between the initial (t = 0r g /c) and final state (t = 25000r g /c) is ≃ 5 × 10 −3 .In order to test our stellar boundary conditions and driving criterion, we have also performed flat-spacetime isolated neutron star simulations for the aligned (χ star = 0 • ) and the orthogonal case (χ star = 90 • ).We find that the Poynting flux for these two cases closely follow the scaling obtained in the force-free simulations by (Spitkovsky 2006) with, in our case, k 1 ≃ 1.1 and k 2 ≃ 0.9.
For numerical integration, we use the total variation diminishing Lax-Friedrich scheme for fluxes along with Piecewise Parabolic reconstruction scheme and second order modified Euler time-stepper.Our simulation domain extends from −400r g to 400r g in x, y, and z.1 .The stellar radius for all the runs is set to R ⋆ = 4r g .The base resolution of the domain is 128 × 128 × 128 cells with 7 levels of AMR.This allows us to resolve dynamically interesting regions with an effective resolution of 8192 3 cells; in particular, the stellar radius is resolved by over 40 cells.Unless stated otherwise, throughout this letter we adopt geometric units with G = c = 1.

RESULTS
We have carried out five 3D simulations in total.In order to observe the effects of stellar magnetic inclination (χ star ) on the accretion dynamics, we performed four different simulations for χ star ∈ [0 • , 30 • , 60 • , 90 • ] and Ω = 0.03.We have also carried out one simulation to explore dynamics in the propeller case for a star with an inclined magnetosphere (discussed in Section 3.2).
The onset of the simulations is similar to the 2D case: as we start rotating the star, an Alfvén wave is launched from the stellar surface and sweeps across the domain.Behind the Alfvén wave, the solution approaches a force-free equilibrium.Meanwhile, MRI turbulence in the disk causes the inner disk to move inwards until its ram pressure is balanced by the magnetic pressure of the stellar dipole.Fostered by the anti-parallel configuration of the initial star and disk field configuration, the magnetosphere reconnects with the disk field which opens up initially closed dipolar fieldlines and forms accretion columns.Magnetospheric accretion starts at ∼ 3700r g /c.
Figure 1 shows the inner magnetosphere at a quasistationary state for different stellar inclinations with µ = 20 (µ is the magnetic dipole moment) and Ω = 0.03.The first two rows in Figure 1 show density and magnetization profiles in the z-x plane for different inclinations.The third row shows plasma beta in the equatorial (x-y) plane.For the aligned case, the star accretes through two almost symmetrical accretion columns in the upper and lower magnetic hemispheres2 .As we start inclining the stellar magnetic field, similar to the aligned case, the disk opens up the initially closed stellar flux, but now the accretion column further away from the equatorial disk becomes increasingly weaker with higher inclinations and all the matter flows through the column closer to the disk.
In the inclined cases, it is striking that the dominant accretion column is formed below the upper polar cap for x > 0 and above the lower polar cap for x < 0. In fact, as we will discuss further in Section 3.1, the accretion column moves towards the location of anti-parallel magnetic fields.This is most evident in the case of χ star = 90 • (Figure 1 (h)) which shows the column notably below the equatorial region (for x > 0).
Tongues of accreting gas form at the disk-magnetosphere boundary due to Rayleigh-Taylor instabilities (RTI) (Arons & Lea 1976a,b;Scharlemann 1978;Romanova et al. 2003;Kulkarni & Romanova 2008).These tongues of matter result in accretion via secondary columns close to the equatorial plane (most pronounced in the aligned case).The instabilities at the disk-magnetospheric boundary allow matter to sweep between the magnetic fieldlines, resulting in a smaller effective magnetospheric radius compared to the 2D simulations.For χ star ≳ 30 • , RTI plays a lesser role as the matter is channeled directly onto the inclined magnetic poles.This was also observed in previous 3D star-disk simulations by Kulkarni & Romanova (2008).With increasing stellar magnetic inclinations, the magnetic pole comes closer to the equatorial disk, and it becomes energetically favorable to channel matter directly onto the magnetic poles.As a result of this, we see that the oblique magnetosphere becomes more extended in y-direction (the axis perpendicular to the plane of spin and stellar dipole moment) in the (x-y)-plane.

Jet Formation
The accretion disk collimates the open stellar magnetic flux to Poynting-flux dominated jets.Here we define the jet as the magnetically dominated region (σ > 1) around the axis.Figure 2 shows the time series of the mass accretion rate at the stellar surface and jet power extracted at 50r g for different stellar inclinations.The mass accretion rate and jet power are defined as follows, where The jet power decreases consistently with increasing magnetic inclination.However, the mass accretion being governed by the disk does not show any such correspondence.For χ star ∈ [0 • , 30 • , 60 • , 90 • ], the average mass accretion rate is Ṁ = [0.441,0.795, 0.708, 0.394] with a significant variance in each case (c v = SD(⟨ Ṁ⟩) ⟨ Ṁ ⟩ = [0.202,0.13, 0.173, 0.14]).The jets formed in all the cases explored here are quite stable and symmetric in both upper and lower hemispheres (for ∼ 25000r g /c).   the quantities are averaged over t ∈ [15 000, 25 000]r g /c.We also show corresponding quantities for the isolated (nonaccreting) case.For accreting neutron stars with aligned rotation and magnetic axes, previous axisymmetric studies suggest that the disk enhances the isolated open flux, leading to enhanced jet power (Parfrey et al. 2016;Parfrey & Tchekhovskoy 2017;Das et al. 2022).The wind power scaling with stellar magnetic inclination for the isolated case is given by L ≃ L aligned (1 + sin 2 χ star ) (Spitkovsky 2006;Tchekhovskoy et al. 2013).However, for neutron star jets, we find that the variation in power with inclination is opposite to that of the isolated wind power scaling.
There are several effects that come into play regarding the jet power of the oblique rotator.Already the isolated case which was studied in detail by Tchekhovskoy et al. (2016) is quite subtle.In essence, by inclining the dipole axis, the distribution of the open magnetic field becomes a nonuniform function of B r (θ, ϕ) which means, there is no direct correspondence between magnetic flux (surface integral of B r ) and Poynting flux (surface integral of (B r ) 2 ).In Tchekhovskoy et al. (2016), around 40% of the power increase was attributed to an overall increased open flux and 60% to the non-uniform "bunching" of the open magnetic field lines.
In our simulations, inclination-dependent equatorial field line bunching plays only a minor role.This is because nonuniform bunching of isolated MHD winds is sub-dominant against the collimating effect from the pressure of the disk.A quantitative comparison of the distribution of Poynting flux and radial magnetic field across the jet is provided in Appendix B.
Another mechanism that affects the jet power is the diskinduced magnetic flux opening.Similar to previous 2D simulations Das et al. (2022), we initialize our domain with an anti-parallel star-disk magnetic field configuration.In the aligned case, the amount of open magnetic flux is then given by the location of the magnetospheric radius in the equatorial plane (Parfrey et al. 2016;Parfrey & Tchekhovskoy 2017;Das et al. 2022).However, in the 3D case, disk-induced flux opening becomes less efficient with magnetic obliqueness.This is consistent with the increasingly larger closed zone seen also in the middle row of Figure 1.Overall, for the accreting cases presented here, the open magnetic flux varies by a factor ≈ 2 from the aligned case to the orthogonal rotator.
To a good approximation, the jet power in our simulations is related to the open magnetic flux as L = ξΩ 2 Φ 2 open with ξ ∈ [0.34, 0.43] which is indicated by the green shaded area in Figure 3.This suggests that disk-induced flux opening determines the jet power which leads to a decrease in power with increasing obliqueness.
We are hence left with one question: what is the reason for the reduced flux opening efficiency for the oblique case?The answer must be sought in the 3D geometry of the stardisk system.In the aligned configuration, disk-connected field lines are subjected to ballooning and flux opening on both hemispheres (e.g.Uzdensky 2004).The newly opened field lines are subsequently collimated into a (northern-and southern-) jet by the lateral pressure of the accretion disk.In the highly oblique case, the situation changes qualitatively: since the poles point towards the direction of the disk, a larger fraction of already open magnetic field lines are able to reconnect with the disk fields.The sketch in Figure 8 shows the contrasting scenarios for the aligned case and the orthogonal rotator.In the former case, two newly opened field lines are created, one contributing to the northern-and one to the southern-jet.In the latter case, no additional open flux is obtained when the disk reconnects with already open field lines (green boxes in Figure 8), however the reconnected field line can move across the disk to end up in the jet on the other hemisphere.Since the accretion columns are bounded by the last open field line -one of which has just moved across the disk -we can now understand why the accretion columns shown in Figure 1 are tilted towards the anti-parallel region of the star-disk system (top-left and bottom right quadrants).Also in the case of the orthogonal rotator, the closed dipolar stellar field lines are able to reconnect with the disk fields and open up.This situation is highlighted by the two inner purple boxes in panel (d) of Figure 8.In our simulation however, this flux opening is almost exactly compensated by the absorption of initially open flux by the equatorial accretion disk.Thus in the case of the orthogonal rotator, we end up with zero net flux opening (cf. Figure 3).
Efficient reconnection of the anti-parallel fields leads to the situation that the simulated jets are fed by only one pole per hemisphere even for stars with extreme obliqueness of 90 • .To further elucidate this point, we show the footpoints of the open magnetic field in Figure 4.In the left column (a), we show the initial wind-like configuration recovered in our simulations before the onset of MRI-driven accretion.In the right column (b) we illustrate the situation in the late stage of the simulation.The colors represent foot-points of the fieldlines contributing to the upper (blue) and lower (red) jet.In practice, field-lines are traced up to r > 50r g and are tagged as part of the jet as long as they remain in magnetically dominated σ > 1 regions.
For an isolated pulsar, in the aligned case (χ star = 0 • ), for a given rotational hemisphere, the open-flux consists of fieldlines that connect to only one magnetic pole.But with increasing magnetic inclination, more fieldlines of the wind connect to the "other" pole which gives rise to the striped wind (Michel 1971).This case is illustrated in Figure 4a).On the surface of the star, we tag footpoints of field lines that end up in the upper (rotational) hemisphere (z > 0) as blue and as red if they end up in the lower hemisphere (z < 0).With increasing inclination, more field lines with opposite polarity are found in either hemisphere (hence more admixture of red and blue tagged footpoints respectively).After the onset of accretion, the disk absorbs some of the initial open flux and even moves open field lines across the disk (as sketched in Figure 8).Ultimately this leads to a reduction of flux opening efficiency and a single polarity outflow in either hemisphere (Fig 4b).
During the review process, we came across another set of 3D GRMHD accreting neutron star simulations which ex- plores the stellar magnetic inclination parameter space by (Murguia-Berthier et al. 2023).They also report the same trend in jet-power with magnetic inclination for the antiparallel star-disk geometry for their µ = 10 − 20 cases.

Propeller Jets
Along with a parameter sweep on the magnetic inclination angle, we have also performed one run to explore the dynamics in a strong propeller regime (r m ≃ 0.7r lc , ≃ 1.9r co ).We initialize our simulation with a similar setup as discussed in section 2 (µ = 30, Ω = 0.05) and set the stellar magnetic inclination to 60 • .In order to maintain numerical stability in this more challenging scenario, the region inside the lightcylinder is initialized with σ t = 60 and β t = 0.016 and the pressure and densities are floored such that β > 0.003 and σ < 300.
The morphology of the jets for χ star = 60 • is shown in Figure 5 for the accretion regime (a) and propeller regime (b).Similar to the previous cases, the disk collimates the open fieldlines to form a jet.However, there are some differences in the jet structure in these two regimes.First, due to the larger magnetospheric radius, the jet opening angle in the propeller regime is larger compared to the accreting cases.Second, similar to the striped wind from isolated oblique rotators, the misaligned propeller shows current sheets within the jet.The current sheets first occur once the flow crosses the light cylinder but, since the entire flow is collimated by the ambient disk pressure, they are not restricted to the traditional striped wind region located ±χ star around the equa- tor.In our highly oblique accreting scenarios, stripes are also observed at early times, but disappear once the fieldlines in anti-parallel regions reconnect with the disk.At large distances, the jet structure is similar to the aligned case.In the strong propeller regime, reconnection between the open polar fieldlines of the star and the disk magnetic fields (green boxes in Figure 8) is inhibited.As a result, the stripes survive within the jet.This is further illustrated by a field line foot-point map (Figure 6).Unlike the accreting cases, here we can see that at a quasi-stationary state, the jet has two polarity fieldlines in both hemispheres and hence the stripes in the jet.

DISCUSSION & CONCLUSION
Our simulations should be applicable to the hard state of accreting millisecond pulsars.They show that the formation of jets in this parameter regime is quite robust, even for highly oblique cases of the star-disk magnetospheric system.Compared to black hole jets where the entire magnetic field must be sourced outside of the event horizon, the inclination of the internal field of the neutron star offers an additional free parameter.Hence for otherwise identical initial conditions and for the same mass accretion rate, we observe a reduction of jet power by a factor of ≃ 2.95 between the orthogonal rotator and the aligned configuration.This is explained by the less efficient flux opening of the inclined star-disk system.Assuming that neutron star jets are indeed formed by the mechanism studied here, unless all neutron star systems conspire to the same obliqueness angle, we thus expect their jets to be more widely scattered in power than black hole jets.Recent high-sensitivity radio observations lend some support to this prediction (van den Eijnden et al. 2021).At the same time, the jet power in all models presented here is sufficiently strong to explain radio emission from both weak-and strong magnetic field sources such as Cir X-1, Sco X-1 and Swift J0243.6+6124(e.g.Fender et al. 2004;Bradshaw et al. 1999).In particular, scaled to physical units and assuming a polar magnetic field strength of 10 9 G, the χ star = 90 • case (the weakest jet presented here) corresponds to a jet power of ≃ 8 × 10 36 erg s −1 , almost two orders of magnitude above the minimum jet power estimated by Fender et al. (2004) for Cir X-1.We refer the reader to Das et al. (2022) for further discussion of the jet power scaling.
To compare with the earlier 2D results, we have performed one 2D simulation with matched initial conditions to the χ star = 0 • case shown here3 .In 3D simulations, due to Rayleigh-Taylor instabilities developing at the disk-magnetospheric boundary, the magnetospheric radius for χ star = 0 • is much closer to the star compared to the 2D case.However, the mass accretion rate in both cases are quite similar, 0.318 and 0.279 (code units) for 2D and 3D respectively.The variability in the mass accretion rate in the 2D (c v = 0.38) simulations is quite large compared to the 3D (c v = 0.154) simulations.This occurs because in 3D geometry, the accretion flow can always penetrate the magnetosphere in between the magnetic fieldlines leading to a comparably stable surface accretion rate.Similar to the mass accretion rate, the jet power in both these cases is also quite close.The average jet power in 2D and 3D are 0.011 and 0.015 respectively (both in code units).To summarize, for aligned accreting neutron stars, although the accretion dynamics in the inner magnetosphere in 3D change compared to the axisymmetric simulations, these variations are not transferred into the jet.
Another possible factor that might impact the behavior of jets in the 3D simulations is the dependence of flux opening on the initial star-disk magnetic field geometry.In all of our simulations, we have considered anti-parallel star-disk magnetic fields which in principle allows the disk magnetic field to reconnect with the stellar magnetic field which leads to enhanced flux opening (for lower inclinations).Recent 3D GRMHD simulations of the aligned case by Parfrey & Tchekhovskoy (2023) suggest that for r m ≪ r co , the difference between parallel and anti-parallel star-disk geometry is quite prominent, with the parallel geometry resulting in less jet-power compared to the anti-parallel case.However, the disk-polarity dependence is reduced for higher inclinations and when r m ≃ r co such that simulations with both disk-polarities result in comparable jet powers (Parfrey & Tchekhovskoy 2023;Murguia-Berthier et al. 2023).
The simulation of the propeller regime shows that the neutron star jets resulting from an inclined star can contain stripes of magnetic field with opposite polarity.In general, the stripes are possible if the jet is comprised of open field lines that originate from opposite magnetic poles.In the isolated case, the striped wind geometry is well understood and confined to the angular range π/2−χ star < θ < π/2+χ star .Since the disk collimates the jet however, we find that the stripes in the accreting scenario can preside much closer to the polar axis (for χ star = 60 • within 8 • < θ < 20 • and 160 • < θ < 172 • ).
Even though the striped jets only arose in the strong propeller simulation, striped jets might survive also in accreting regimes.As long as opposite polarity fieldlines are collimated into a jet of one hemisphere, current sheets are expected to form within the outflow.This might also happen for the accreting case depending on the reconnection efficiency of the star-disk magnetic field, the stellar magnetic inclination, and the initial magnetic field in the disk.
The possibility of striped jets is interesting since it provide a means for magnetic energy dissipation via magnetic reconnection in the jet.For isolated stars with misaligned spin and magnetic axis, the dissipation in the striped wind has been studied by many authors (Coroniti 1990;Lyubarsky 2003;Komissarov 2013), most recently by Cerutti et al. (2020); Hakobyan et al. (2023) who find that the global structure of the stripped wind follows the split-monopole prediction (Bogovalov 1999) and in the stripped region, plasmoid dominated magnetic reconnection leads to efficient magnetic energy dissipation throughout the wind.In the case of misaligned stars, assuming a split-monopole magnetic wind, the dissipation for r > r lc is given by L(r) = L 0 (1 − β rec ln(r/r lc )).Here, β rec is the reconnection rate, and L 0 is the Poynting flux in a dissipation-less split-monopole magnetosphere (L 0 = 2cB 2 star r 4 star /3r 2 lc ) Cerutti et al. (2020).In our case, for χ star = 60 • due to the flux collimation by the accretion disk, the stripes only survive within 8 • < θ < 20 • and 160 • < θ < 172 • .For a given reconnection speed β rec , we can estimate the dissipation of the current sheets in the jets following Cerutti et al. (2020).With the above quoted extent of the striped region, this leads to L(r) = L 0 (1 − 0.0033β rec ln(r/r lc )).For β rec ∼ 0.1 − 0.2, the resulting dissipation radius is very large (> 10 30 r g ) which indicates that dissipation of stripes is dynamically unimportant.Never the less, the formation of current-sheets in the magnetically dominated jet offers a mechanism for particle acceleration which might power flaring emission in the "striped state".The striped state should emerge when star-disk coupling is moderate: e.g. in the propeller regime when the accretion rate is low and the disk is geometrically thick and able to collimate a jet.In addition, the reconnection of the stardisk system is reduced for parallel field configurations which could give rise to striped jets even for strongly accreting systems.
The 3D simulation presented here mark a first step to understand neutron star jet formation from oblique systems.More work is needed to address the dependence of the jet properties with the parameters of the system, foremost the magnetospheric radius and initial magnetic field geometry.
PD and OP acknowledge funding from the Virtual Institute for Accretion (VIA) within NOVA (Nederlandse Onderzoeksschool voor Astronomie) Network 3 'Astrophysics in extreme conditions'.Simulations have been carried out in part on the Dutch national e-infrastructure with the support of SURF Cooperative (project number NWO-2022.019),Swedish cluster Dardel as a part of PRACE (DECI-17) allocation (project number 17DECI0051) and the HELIOS cluster of the Anton Pannekoek Institute for Astronomy.We thank Anna Watts and Nathalie Degenaar for valuable comments.The line element of the Schwarzschild coordinates in standard spherical coordinates (t, r, θ, ϕ) can be written as, With the following transformations, the covariant components of the cartesian corotating Schwarzschild metric can be written as, Here Ω is the angular corotation frequency.In our simulations, this is always set to Ω star .

Figure 1 .
Figure 1.Inner magnetosphere for different rotation and magnetic axis inclinations at t = 12000rg/c.The upper panel shows log10ρ and the middle panel shows log 10 σ in the z-x plane.The lower panel shows log 10 β in the equatorial plane (x-y).The gray dashed lines show the corotation radius.The black lines in the lower panel shows the β = 1 contour.
Fig 3 (a) shows variation of jet power (P jet ) and panel (b) shows variation in open flux (Φ open ) with inclination.All

Figure 2 .
Figure 2. (a) Mass accretion rate at the stellar surface for different stellar magnetic inclinations.(b) Evolution of jet power with time for different magnetic inclinations.The jet is defined as σ = 1 contour.The mass accretion rate and the jet power is extracted at 4.2rg and 50rg respectively.

Figure 3 .
Figure 3. Variation of Jet power (a) and open flux (b) with magnetic axis inclination (χstar).The vertical axis shows the mean quantities normalized by the aligned isolated cases for t ∈ [15000,25000]rg/c.The red dashed lines represent the variation of power and open flux in the isolated cases.The values of isolated open flux as a function of χstar are adopted from Tchekhovskoy et al. (2013).

Figure 4 .
Figure 4. Surface maps of the open magnetic fieldlines powering the jet for different stellar magnetic axis inclinations in the accreting regime.The blue and red colors represent the fieldlines ending up in the jet in the northern (z > 0) and southern hemispheres (z < 0) respectively.The first column (a) represents the magnetic foot-points of the open fieldlines at an initial time (mimicking the isolated case) and the second column (b) shows the foot-points of the open field lines in the jet at a later time (after reaching the quasi-stationary state).The green arrow indicates the direction of the magnetic moment.For a dipolar geometry, the polarity of the magnetic field flips at the magnetic equator (outlined by the green circle), b r > 0 (outgoing fieldlines) in the hemisphere along the magnetic moment and b r < 0 (in-going fieldlines) in the rest of the hemisphere.

Figure 6 .
Figure 6.Surface maps of the open magnetic fieldlines powering the jet in the upper (blue) and lower (red) hemispheres.Here χstar = 60 • , µ = 30 and Ω = 0.05.The two panels show the surface maps at two different times (a) Early (b) Quasi-stationary state.All the labels are the same as in Fig 4.

Figure 7 .
Figure 7. θ profiles of different quantities for χstar = 0 • (blue) and 90 • (red).Panel (a) shows the time evolution of the quantities at the stellar surface.Panel (b) and (c) show the time evolution at 50rg in the regions with σ > 1.The dotted, dot-dashed lines in (a) show ⟨B r 2 ⟩ ϕ at t ∈ [0, 20000]rg/c.The solid and dashed lines in panel (b) and (c) show ⟨B r 2 ⟩ ϕ at t = t0 and ⟨T r t,EM⟩ϕ respectively.Here t0 = 200rg/c (b) and t0 = 20000rg/c (c).

Figure 8 .
Figure 8. Stellar and disk magnetic field interaction for χstar = 0 • (left) and χstar = 90 • (right).The upper panel shows the initial magnetic field and the lower panel shows the interaction at an intermediate time.The blue and red lines in the lower panel show the fieldlines in the upper and lower jet respectively.The magnetic reconnection of fieldlines with opposing polarity is indicated by the purple and green boxes.The purple box represents the reconnection of the initially closed stellar zone, while the green box represents reconnection of the initially open stellar fieldlines.The dotted lines in panel (c) represents the newly open fieldlines as a result of reconnection of the stellar closed zone.The dashed lines in panel (d) show how the initially open fieldlines in the lower right and upper left get re-arranged by the disk due to magnetic reconnection.The accretion column has a shifted position as a consequence of the re-arrangement of polar field lines.