Radio Emission of Pulsars. I. Slow Tearing of a Quantizing Magnetic Field

The pulsed radio emission of rotating neutron stars is connected to slow resistive instabilities feeding off an inhomogeneous twist profile within the open circuit. This paper considers the stability of a weakly sheared, quantizing magnetic field in which the current is supported by a relativistic particle flow. The electromagnetic field is almost perfectly force-free, and the particles are confined to the lowest Landau state, experiencing no appreciable curvature drift. In a charge-neutral plasma, we find multiple branches of slowly growing tearing modes, relativistic analogs of the double tearing mode, with peak growth rate $s \gtrsim 4\pi \widetilde k_y J_z/B_z$. Here, $B_z$ is the strong (nearly potential) guide magnetic field, $J_z$ the field-aligned current density, and $\widetilde k_y$ is the mode wavenumber normalized by the current gradient scale. These modes are overstable when the plasma carries net charge, with real frequency $\omega \sim s\cdot |n_0^+ - n_0^-|/(n_0^+ + n_0^-)$ proportional to the imbalance in the densities of positive and negative charges. An isolated current sheet thinner than the skin depth supports localized tearing modes with growth rate scaling as (sheet thickness/skin depth)$^{-1/2}$. In a pulsar, the peak growth rate is comparable to the angular frequency of rotation, $s \gtrsim 2\widetilde k_y \Omega$, slow compared with the longitudinal oscillations of particles and fields in a polar gap. The resistive modes experience azimuthal drift reminiscent of sub-pulse drift and are a promising driver of pulse-to-pulse flux variations. A companion paper demonstrates a Cerenkov-like instability of current-carrying Alfv\'en waves in thin current sheets with relativistic particle flow, and proposes coherent curvature emission by these waves as a source of pulsar radio emission.


INTRODUCTION
An outstanding question about stellar magnetic fields centers on the role that internal resistive instabilities play in redistributing magnetic twist and the associated electric current. This phenomenon has long been encountered in magnetic plasma confinement devices (tokamaks), where an elaborate theory of these instabilities has been developed (White 2013). High-order resistive instabilities involving multiple magnetic tearing surfaces are known to play a central role in establishing the large-scale twist profile within the plasma torus. This process can be viewed, in some respects, as a relaxation to a lower-energy state while conserving magnetic helicity (Taylor 1986).
As regards stellar magnetospheres, most studies of magnetic tearing have focused on the effect of divergences in the magnetic shear that is imposed at the surface of the star, or on collisions between discrete magnetic structures (Parker 1972;Priest & Forbes 2000;Nalewajko et al. 2016). A few investigations have been made of the microscopic underpinnings of turbulent relaxation in stellar magnetospheres (e.g. Browning et al. 1986;Boozer 2020), but none directed toward the relativistic magnetosphere of a neutron star.
The case of a neutron star presents unique challenges and interesting opportunities. The main thesis of this paper is that the twist profile imposed on a pulsar magnetosphere by the global charge flow is susceptible to internal tearing instabilities. The open magnetic flux bundle is very narrow near the star; high-order instabilities are less sensitive to the global connectivity of the magnetic field, and in principle will arise even in zones that do not sustain electron-positron pair creation. We find that internal tearing in the inner pulsar magnetosphere is slower than the plasmoid instability seen near the magnetic separatrix in numerical simulations of the outer magnetosphere and wind zone (Cerutti & Philippov 2017;Philippov et al. 2019). The characteristic scale is small -around the plasma skin depth of ∼ 0.1 − 100 cm -and the growth time comparable to the rotation period of the star.
A companion paper (Thompson 2021) identifies a secondary, Alfvénic instability of small-scale current sheets, which shows promise in explaining the charge clumping that is needed to sustain coherent curvature emission in the radio band.

Tearing in a Quantizing Magnetic Field
This paper has a few interconnected goals. The first is to extend the theory of magnetic reconnection in a collisonless, relativistic plasma to the quantizing regime, where a strong guide magnetic field is present, the magnetization may approach ∼ 10 20 , and the effects of particle gyromotion and curvature drift are essentially negligible. We adopt a kinetic approach, as this problem cannot reliably be approached by a relativistic form of resistive magnetohydrodynamics. Previous kinetic investigations of non-relativistic (Galeev & Zelenyǐ 1976;Drake & Lee 1977;Basu & Coppi 1981;Ottaviani & Porcelli 1993) and relativistic (Zelenyi & Krasnoselskikh 1979;Zenitani & Hoshino 2008;Sironi & Spitkovsky 2011;Cerutti et al. 2013;Sironi & Spitkovsky 2014;Hoshino 2021) reconnection have focused on plasma configurations where curvature drift is not negligible, either because particles gain relatively high energies, or because the guide magnetic field is not very strong.
In the relativistic regime with a strong guide field, we find that a localized tearing mode can be excited at an isolated current sheet if the sheet thickness is comparable to or smaller than the skin depth. In addition, we allow the plasma to have a finite charge density, finding that the tearing instability becomes an overstability: magnetic islands experience a general drift along the tearing surface.
A second goal is to build theoretical tools that will be useful in studying the spontaneous relaxation of a smoothly twisted magnetic field in a stellar magnetosphere. This process is, to a limited extent, constrained by the tying of magnetic field lines at the stellar surface; the endpoint is generally not a simple linear force-free magnetic field, with current density and magnetic flux density connected by a spatially uniform constant (e.g. Russell et al. 2015). A smoothly twisted field is, nonetheless, generally susceptible to some degree of relaxation. In a tokamak, where the relaxation process is less well described by Taylor theory than in a reversed-field pinch, configurations with magnetic twist peaking off the toroidal axis will quickly relax to a more centrally peaked configuration. Internal tearing modes feeding off multiple interacting current sheets appear to play a significant role. They are a major focus of our study.
The plasma in the magnetosphere of a neutron star is so dilute that electron inertia can play a stronger role than Coulomb collisions in mediating the breakdown of magnetic flux freezing. We also focus on the case where ions are absent, corresponding to a negative corotation charge in the magnetosphere. The ordering of the electron gyro-scale and Debye scale is reversed in the relativistic plasma, as compared with the more familiar non-relativistic case. We may therefore contrast our approach with Galeev & Zelenyǐ (1976), who in their treatment of non-relativistic reconnection considered tearing modes that are strongly influence by electron gyromotion, even while the magnetic gradients are restricted to scales much exceeding the Debye length. In a more complete treatment of non-relativistic reconnection (Drake & Lee 1977;Basu & Coppi 1981), both the electron gyroscale and Debye scales are resolved, but the tearing eigenfunction incorporates an outer ideal zone where ions provide a neutralizing charge. In contrast with the relativistic case, the electrostatic potential contributes minimally to the parallel electric field. Ottaviani & Porcelli (1993) give a simplified treatment of interacting tearing surfaces that resolves the Debye scale and neglects collisions, but also incorporates heavy, neutralizing ions.
In the case of a relativistic and dilute plasma, where collisions are slow and ions are absent, we find that the outer tearing mode never fully reaches the ideal regime. Growing modes are obtained in the presence of multiple interacting (periodic) tearing surfaces. Similarly to the collisionless, non-relativistic case (e.g. Ottaviani & Porcelli 1993), we find that when the separation between tearing surfaces much exceeds the Debye scale, the approximation of a uniform magnetic flux perturbation breaks down near each tearing surface; in a related effect, the mode displacement is strongly concentrated in the bulk. However, we also find a tower of overtone modes with finite wavenumber transverse to the current sheet and growth rates nearly as large as the fundamental mode. These long range oscillations in the perturbed magnetic field are ascribable to the absence of heavy ions and breakdown of screening.

Pulsars and Magnetars
Our third goal is to investigate the role that the spontaneous relaxation of twist may play in sourcing radiative emissions from a neutron star. Twist is injected into the magnetosphere of a radio pulsar by an escaping flow of corotation charge (Goldreich & Julian 1969;Spitkovsky 2006); in a magnetar, by an internal magnetoelastic instability (Thompson et al. 2002;Chen & Beloborodov 2017).
The radio emission of rotating neutron stars exhibits (i) very short-timescale variations (microstructure), implying strong angular variations within a single pulse (Graham-Smith 2003); and (ii) much slower pulse-to-pulse variations in the envelope of a radio pulse, including a secular drift in the phase distinct sub-pulses (e.g. Deshpande & Rankin 2001). The tearing mode analysis presented here makes some suggestive connections with these well-known properties of pulsars. First, tearing is found to be concentrated on small lengthscales, comparable to the skin depth of ∼ 0.1 − 100 cm in the open pulsar circuit. Second, the growth of the tearing instabilities is constrained by the strong background dipolar magnetic field, being limited to a rate comparable to the angular frequency of rotation. In addition, high-order tearing modes in a nearly azimuthally symmetric plasma can form a discrete ring of current 'knots' (Bierwage et al. 2005;Wang et al. 2007). The tearing process turns out to be an overstability in the presence of a finite charge density, as is expected in the pulsar magnetosphere. Angular drift is therefore an intrinsic property of such structures.
The main implication of this work for pulsars is that the process of radio emission is generally top-down, involving internal tearing of a large-scale magnetic field. It has long been suspected that the emission of Crab-like pulsars arises from the separatrix between open and closed magnetic field lines, or the equatorial current sheet in the wind (e.g. Gruzinov 2007;Cerutti & Philippov 2017). Internal tearing modes are slowly growing in the open circuit, but still potentially play a potent role in adjusting the current profile over a single rotation.
The magnetosphere of a magnetar appears to sustain a relatively stronger current than is encountered in a rotationpowered pulsar. The instabilities found here also provide a starting point for investigating the the redistribution of currents within the magnetosphere, which is mediated by a combination of internal tearing (Thompson 2008) and ohmic dissipation (Beloborodov 2009), the relative importance of which is not understood.

Plan of the Paper
The plan of this paper is as follows. Section 2 sets the stage by reviewing the plasma properties expected in the magnetosphere of a radio pulsar or magnetar, and some basics about internal tearing in a magnetized plasma. Next, in Section 3, we clarify the connection between charge conservation along relativistic magnetic flux tubes and vortical deformations of a nearly force-free and quantizing magnetic field. Our kinetic approach to linear perturbations is described in Section 4. Growth rates are derived in Section 5 for linear tearing modes in a neutral plasma supporting a periodic flux component with a wavelength somewhat larger than the skin depth. This result is generalized in Section 6 to a single current sheet. A shift from instability to overstability in a plasma with net charge density is demonstrated in Section 7. The effect of line tying on the tearing mode and the plastic deformation of the upper neutron star crust are discussed in Section 8. A summary of our results is presented in Section 9, along with some implications for radio pulsars and magnetars.
Throughout this paper, we adopt the shorthand X = X n × 10 n to describe the normalization of quantity X in c.g.s. units.

COLLISIONLESS PLASMA IN A QUANTIZING MAGNETIC FIELD
In this paper, we investigate the collective dynamics of charged particles and a quantizing magnetic field (B 10 8 G in radio pulsars and 10 14 G in magnetars). The particles relax rapidly to the lowest Landau state and, since B is deformed only slowly, are guided along this field. In other words, one may locally define a frame where the electric field E B and E 2 B 2 , although globally this is not possible. The particle flow is driven by a gentle transverse shearing of the magnetic field. The self-inductance of the system is enormous: the non-potential component of the magnetic field may carry several orders of magnitude more energy than the charges. The particle distribution function is locally one-dimensional and highly nonthermal.
The ratio of plasma frequency to gyrofrequency is extremely small. In a frame where the plasma has vanishing net momentum, where σ = B 2 /4π(nγmc 2 ) is the magnetization, γ is a characteristic Lorentz factor, and q, m, n are the particle charge, mass and number density. We will therefore neglect curvature drift and other effects of gyromotion.
We consider two cases: (i) a charge-balanced plasma (charges ±q are present in equal numbers) and (ii) a plasma composed of charges of a single sign, with net charge density ±|q|n (as appropriate to some parts of the pulsar magnetosphere where e ± pairs may be relatively rare).
Our analysis of tearing instability focuses on the simplest case, where m and |q| are the same for all plasma particles. After deriving the basic kinetic equations, we first consider the simplest case of a charge-balanced plasma (Sections 5 and 6). Next, deviations from charge neutrality are treated as a perturbation to the charge-neutral case (Section 7).

Pulsars
Our approach is motivated by the application to radio pulsars and magnetars. When a neutron star rotates with frequency Ω, its corotating magnetosphere supports a net charge density ρ c ∼ ∓ΩB/2πc. Allowing for M ± pairs per corotation charge (and for the possibility that the seed charges are protons rather than electrons) gives where P = 2π/Ω is the spin period and −e and m e are the electron charge and mass. The plasma skin depth is some 7 to 9 orders of magnitude larger than the gyroscale, and is given by where g(m, γ, γ ± , M ± ) = [(m e /m) γ −3 + 2M ± γ −3 ± ] 1/2 and ... denotes an average over the one-dimensional momentum distribution. The parts of the polar magnetic flux bundle where k −1 p ∼ 0.1 − 10 2 cm are of greatest interest for radio emission.
Models of nebular synchrotron radio emission around young, rapidly spinning pulsars provide constraints on the pair multiplicity M ± , in at least some parts of the pulsar circuit. In the case of the Vela pulsar, M ± ∼ 10 5 is inferred (de Jager 2007; Bucciantini et al. 2011). Theoretical models of pair creation in the polar cap of a neutron star depend on the self-consistent solution for E · B, and imply a range of multiplicities (Hibschman & Arons 2001;Timokhin & Harding 2019). The outcome of these models is influenced by secondary effects such as relativistic frame dragging (Muslimov & Tsygan 1992;Philippov et al. 2015), and will also be influenced by changes in the local current density driven by the instability described here.
The role of Coulomb collisions in magnetic tearing in the pulsar circuit depends in interesting ways on M ± , as well as on the strength of the magnetic field relative to the QED field B QED = m 2 e c 3 /e = 4.4 × 10 13 G. Collisions that are trans-relativistic in the center-of-momentum frame have a total cross section comparable to the unmagnetized value when B B QED . Comparing the collision rate s col with a mode growth rate s ∼ 4π|J|/B, one finds s col /s ∼ α em (B/B QED )M ± when J ∼ ρ c c (and the corotation charge is supplied by electrons). Thus, collisions will have a negligible effect in the parts of the pulsar circuit where the e ± multiplicity is modest, and especially where the particle kinetic energy remain below threshold for triggering a e ± pair cascade. In Paper II, we show that matching the drift rate of internal tearing modes with the angular drift rate of radio sub-pulses suggests low or vanishing M ± within the drifting magnetic structures.
A new channel for e ± backscattering opens up when B 5 B QED : a single gamma ray produced by e ± annihilation will quickly re-convert back to a pair, e + + e − → γ → e + + e − (Thompson & Kostenko 2020). This is mediated by the resonant quantum-electrodynamic s channel, with a cross section ∼ 10 2 times larger than for the Coulomb-like t channel. This suggests that collisions will have a stronger effect on internal tearing in the magnetosphere of a magnetar than in ordinary radio pulsars.

Magnetars
So far, two complementary models have been developed for the plasma state surrounding a magnetar. It is possible that a counterstreaming flow of electrons and positrons is sustained in a double-layer structure by the resonant scattering of blackbody keV-energy photons flowing from the heated stellar surface (Beloborodov & Thompson 2007). In this case, the plasma is nearly collisionless and the e ± flow is close to what is required to sustain a large-scale (e.g. hemispheric) magnetic twist, n ± |∇ × B|/4πe. Hence, the magnetization σ ∼ eB B /γm e c 2 ∼ 10 16 B 15 B,7 /γ 3 , High-order tearing modes associated with overlapping rational surfaces are triggered when the twist peaks off the main toroidal axis. The analysis in this paper focused on the stability of high-wavenumber current gradients that are superimposed on the mean helical magnetic field. We consider a local cartesian patch, with guide magnetic field oriented in theẑ direction (corresponding approximately to directionφ in the torus andr in the open circuit of a radio pulsar). In a magnetic confinement device, these modes drive the redistribution of twist across the small axis of the torus; in a radio pulsar, we propose that they mediate a release of non-potential magnetic field energy and are the underlying instability driving coherent radio emission.. The linear stability analysis presented here is combined with a conjecture that the redistribution of twist across magnetic flux surfaces is driven by small-scale current variations, extending down to scales comparable to the skin depth.
where B = B/|∇ × B| is the magnetic twist length and γ ∼ 0.03 eB/m e ckT bb ∼ 10 3 B 15 (kT bb /keV) −1 is the Lorentz factor of the resonantly interacting e ± . A recent re-examination of QED processes in ultrastrong magnetic fields has revealed the possibility of a collisional plasma state with a high resistivity, unlike the most common picture of a pulsar magnetosphere (Thompson & Kostenko 2020). This collisional state would be sustained by frequent pair annihilation and reconversion, e + +e − → γ → e + +e − . The required current density could be sourced by yielding of the magnetar crurst in compact fault-like zones (as suggested by ab initio calculations of magnetic diffusion: Gourgouliatos et al. 2016;Thompson et al. 2017). The state of the e ± plasma is trans-relativistic, quasi-thermal and stable, with density n ± ∼ 15 |∇ × B|/4πe and magnetization about 10 2 times higher than in the relativisitic double-layer model. One also obtains a direct explanation for the presence of a hard X-ray component of the magnetar spectrum: the annihilation of trans-relativisitic e + − e − pairs in magnetic fields B 10 14 G produces a broad, bremsstrahlung-like spectrum of X-rays (Thompson & Kostenko 2020), similar to that observed. The higher density in this state also provides a promising context for collective plasma emission in the IR-optical band, which is observed from quiescent magnetars at a rate far exceeding the expected surface blackbody flux (Kaspi & Beloborodov 2017).

Internal Tearing in Cylindrical and Toroidal Geometry
In this section, we first review some basic aspects of internal tearing in a magnetized plasma with toroidal or cylindrical symmetry, which may be useful for an astrophysical audience. Then we make contact with the nonpotential magnetic field configuration encountered in the polar regions of a neutron star magnetosphere. Our focus in this paper is on the interaction of multiple, closely spaced tearing surfaces, which can be studied using a (locally) cartesian toy model. The analysis in Sections 3-7 is based on such a cartesian model.
A simple point of notational confusion may arise here ( Figure 1). Distinct coordinate systems are adopted in studies of pulsar electrodynamics, where the active magnetic field lines are concentrated in a narrow bundle near the magnetic dipole axis, and in the literature dealing with the breakdown of confinement in toroidal plasmas (tokamaks). In the first case, standard spherical polar coordinates are frequently adopted. In the latter case, the toroidal magnetic field B φ winds around the large radius of the torus, and the poloidal field B θ around the small radius.
In toroidal geometry, the unperturbed magnetic field can be written where r is a locally cylindrical radial coordinate ( Figure 1). This expression manifestly satisfies ∇ · B 0 = 0, and is useful for analyzing current-driven instabilities. The flux functions Ψ φ , Ψ θ are both constant along the magnetic field: In the case of a slender torus, the poloidal field B θ0 does not depend significantly on the large radius R c of the torus, which enters this equation via ∇φ =φ/R c . It is therefore sometimes useful to write fields and currents in terms of the linear poloidal flux, The magnetic twist is commonly described in terms of a 'safety' factor: A toroidally symmetric equilibrium corresponds to q = q(r). Large q (strong 'safety') corresponds to slow rotations of the magnetic field lines in the θ-direction. This regime of weak magnetic twist has a direct analog in the open pulsar circuit.
Enhanced tearing is encountered in configurations in which q(r) has extrema at r > 0, and can be triggered even when 1/q is small (the twist is weak). Important for us is the ability of high-wavenumber oscillations in q(r) to induce tearing, even when superposed on a smooth q profile that dominates the current density. This is typically demonstrated by identifying the mean helical flux surfaceΨ θ , and then considering a linear mode with wavevector k perpendicular to this surface. The remaining flux component ∆Ψ θ sources a more strongly variable field, ∆B θ = R −1 c (d∆Ψ θ /dr)θ, which may exhibit multiple null surfaces ∆B θ = 0. An essential feature of the tearing process is that the perturbation equations directly involve the background current density only through its gradient (Biskamp 2005). This property holds also in the nearly force-free regime, as will become clear in Section 4. As long as the gradient scale ofΨ θ is much longer than that of ∆Ψ θ , the mean current does not significantly modify the growth of the high-wavenumber tearing mode. This is seen directly in our kinetic approach, which follows the perturbed distribution function of the current-carrying charges, whose overall drift speed is fixed by the mean magnetic twist.

High-wavenumber Tearing in Cylindrical Symmetry
With a goal of investigating internal tearing in the pulsar emission zone, consider now a very slender torus and focus on a small segment of length ∆z = R c ∆φ R c . A locally cylindrical coordinate system (r, θ, z) is adopted. The torus may experience long-wavelength deformations that change R c and break toroidal symmetry. This complicates low-wavenumber tearing modes but can have a weaker impact on high-wavenumber modes that feed off small-scale gradients in the magnetic field.
For now, we maintain toroidal symmetry in the background plasma state (invariance under rotations in φ) and look for linear modes depending on the angular variable with wavevector This is chosen to be orthogonal to the mean helical magnetic field, k ·B 0 = 0. The unperturbed magnetic field decomposes as where B z0 = B φ0 varies weakly across the flux bundle in the force-free regime. We will take m and k z both to be large, given that ∆B 0 varies on a small scale. The flux variable decomposes similarly as The unperturbed current profile is then where ∇ 2 ⊥ is the Laplacian in the (r, θ) plane. The first term in J z , which is sourced by the mean helical field, varies smoothly and largely factors out of the tearing mode analysis. The θ-component of J 0 is needed to maintain force-free equilibrium (see Section 3).

Neutron Star Magnetosphere
Now consider the polar region of a neutron star magnetosphere (bottom panel of Figure 1), through which a current is driven either by the corotation charge flow (Goldreich & Julian 1969), or by a subsurface instability, as in the case of magnetars (Thompson et al. 2002). A change of coordinates and notation is now required. The guide magnetic field is identified not with the toroidal component (as in the torus geometry), but with a dipole field centered on the star. We adopt spherical polar coordinates (r , θ , φ ). Near the magnetic dipole axis, this poloidal field is approximately radial, with a small flaring B θ 1 2 θ B r . Within the corotating magnetosphere, a flux surface anchored at polar angle θ s on the star will reach a maximum radius r max = r s / sin 2 θ s (here r s is the stellar radius); its curvature weakens with increasing radius, R c (r ) ∼ (r r max ) 1/2 . Magnetic field lines extending beyond the corotating magnetosphere are confined to an approximately cylindrical bundle of opening angle θ open (r ) (r Ω/c) 1/2 , with a radius of curvature R c ∼ (r c/Ω) 1/2 . The curvature of the open field lines experiences mild oscillations during a single rotation of the star.
The weak magnetic field driving reconnection is now the (approximately) toroidal component. Consider, first, the corotating part of the magnetosphere, within which a flux surface may be twisted through an angle ∆φ between the north and south polar regions. This twist angle generally depends on the footpoint colatitude θ s . The toroidal magnetic field B φ 1 2 ∆φ θ 2 B θ 1 4 ∆φ θ 3 B r is supported by the current density J r = (cθ 2 B r /4πr )∆φ . On the open-field bundle, the toroidal field is sustained by the corotation charge flow J r ρ co c, where ρ co −ΩB r /2πc. This outward flow is compensated by a return current flowing through an annular sheath surrounding the polar cap. In contrast with the current supporting a closed-field magnetic twist, the rotation-driven current J r has the same sign in each hemisphere.
Consider the twist profile as a function of transverse radius = θ r in the simplest case of an axisymmetric current flow. Writing In the simplest case where a central core of the open flux bundle maintains a flat current density, one sees that the profile of qR c is also flat. The tearing instability considered here therefore feeds off gradients in f ( ). Otherwise, the profile of qR c peaks near = 0. Three-dimensional, force-free models of the pulsar magnetosphere reveal strong departures from uniform current density near the center of the open flux bundle, as well as significant departures from axial symmetry (Timokhin & Arons 2013;Gralla et al. 2017). The azimuthally averaged twist profile generally shows an off-axis maximum, which in a tokamak is known to trigger a rapid tearing instability that is driven by the interaction between rational surfaces (White 2013).
In the case of magnetars, global models of yielding in the solid crust suggest that twist is injected in localized, fault-like zones (Gourgouliatos et al. 2016;Thompson et al. 2017). The inhomogeneities in the excited twist are then even strong than those that are imposed by spindown in the open circuit. Tearing modes similar to those investigated here are implicated in the broader spreading of the current through the closed magnetosphere (Thompson 2008).

NEARLY FORCE FREE DYNAMICS
In this section, we summarize the connection between charge conservation and force-free equilibrium. The closed magnetosphere of a magnetar and, to a lesser extent, the open circuit of a rotating neutron star, has an enormous selfinductance: the energy of the charges supporting the background current is miniscule compared with the energy stored in the non-potential magnetic field. The magnetic field is assumed to be weakly sheared, meaning that a strong 'guide' field is present, to be identified with the mean helical field encountered in Section 2.3. The force balance separates into a transverse component, which is dominated by the electromagnetic field, and the longitudinal dynamics of the charges, which is not force free. Our focus here is on the slow drift of charges perpendicular to the guide field.
Here and in the remainder of this paper, we work in cartesian geometry. The guide magnetic field flows in the z direction; hence where A is the vector potential. The symbol ⊥ represents a projection onto the x − y plane. In this section, our focus is on the motion of charges perpendicular to the guide field, and so we only need the perpendicular electric field E ⊥ . The magnetic field is assumed to deform slowly (e.g. E 2 B 2 ), as will be justified after the fact for the linear modes uncovered in this paper.
To an excellent approximation, the charges (of sign ±q) are guided along the magnetic field, with a finite mean speed β ± β ± z that may depend on the sign of q. The transverse drift speed β ± ⊥ is easily obtained from the Lorentz force equation, when Taking the perpendicular divergence of this, combining with the equation of charge conservation, Here ρ ± is the density of positive (negative) charge, The quantity explicitly conserved here is the charge flow guided along a magnetic flux tube, as may be seen by multiplying Equation (16) by a small constant flux δΦ δA ⊥ · B z , and defining the linear charge density δλ = δA ⊥ ρ and current δI z = J z δA ⊥ : The tiny effects of particle inertia can be included, in a manner analogous to the treatment of non-relativistic magnetofluids, by taking the transverse curl of the Euler equation (now with ∂/∂z → 0), Here, n = n + + n − = ρ + /q − ρ − /q and the transverse momentum density is (γmc)nβ which is suppressed compared with the left-hand side by a factor ∼ 4πγmnc 2 /B 2 z . In the case of a pulsar, the suppression is almost complete (∼ 10 −15 − 10 −20 ).
The leading term on the right-hand side of Equation (19) will be familiar from the standard treatment of tearing in a non-relativistic, incompressible magnetofluid (e.g. Biskamp 2005). In the non-relativistic regime, the transverse velocity potential is proportional to the electrostatic potential that describes E × B drift: E = −∇φ and v E×B = z × ∇ imply = cφ/B z and w = ∇ 2 = (c/B z )∇ 2 φ. In the regime of extreme magnetization, the electromagnetic field plays a central role; but we will see that a kinetic description of linear perturbations still is written most succinctly in terms of an effective Langrangian displacement field (Section 4).

KINETIC DESCRIPTION OF LINEAR PERTURBATIONS
The unperturbed state contains a non-potential magnetic field that varies in cartesian direction x but not in y or z, The second equation describes force-free equilibrium, corresponding to (d/dx)(B 2 y0 + B 2 z0 ) = 0. The example investigated in detail below is a periodic shear with a wavelength 2π/k x that can be adjusted with respect to the skin depth: The perturbed system remains uniform in z. The electromagnetic field is the sum of the background (20) and a time-dependent perturbation (labelled 1) that depends on x, y and t, The guide field B z is taken effectively to be a constant, B z B 0 B ⊥0 . In a quantizing magnetic field, the unperturbed particle distribution function f 0 may depend on fewer than three momentum coordinates. For example, if the current sheet is translationally symmetric in one coordinate (y), as is assumed here, then there is no E × B drift out of the plane of the sheet in the background state. In this situation, the perturbed Boltzman equation may involve derivatives with respect to all three spatial coordinates, 1 and, perturbing Equation (15), An electric field is excited in the plane of the current sheet in the perturbed state, so that β x,1 = 0 and one must include a term β ⊥1 · ∇f 0 = β x,1 ∂f 0 /∂x in the Boltzman equation to represent the advection of charges across the magnetic field. The distribution function is taken to be a narrow top hat centered at momentump 0 ≡γ 0β0 mc for q > 0 (and −p 0 for q < 0). When there is a mixture of a positive (negative) charges, each with space density n + 0 (n − 0 = n 0 − n + 0 ), where p 0± =p 0 ± ∆p 0 /2, and Θ is the Heaviside function. 2 We also assume that 1 −β 0 (B y0 /B z0 ) 2 , corresponding toγ 0 B z0 /B ⊥0 . The background particle density is and the skin depth The strong momentum dependence of k p0 is a consequence of the narrow momentum distribution assumed; more generally k p0 ∝ γ −3 0 1/2 . The background is invariant under translations in y, and so each mode is fourier decomposed as Our goal is to calculate the mode growth rate s and real frequency ω as functions of k p0 /k x , k y /k x . Although ω vanishes in the charge-symmetric state (counter-streaming e + and e − with vanishing net charge density, n + 0 = n − 0 ), we find that growing modes are overstable in a charge-asymmetric state. This will be be important in the application to pulsars, in particular with regard to the phenomenon of sub-pulse drift.
The perturbed Boltzmann equation reads where ± labels positive and negative charges. Substituting Equation (15) and integrating the sum consistent with Equation (16). Here the perturbed charge and current densities ρ 1 , J z1 are In regions where the current distribution evolves slowly, one may take ∂/∂t → 0 and Fourier transform in y to get a modified Grad-Shafronov equation, The electrostatic term on the left-hand side generalizes the equation describing the ideal flux perturbation in the standard tearing mode analysis.
Returning to the fully time-dependent problem, we adopt a useful shorthand Then Equation (29) reads where Integrating the quantity q(s + f + 1 − s − f − 1 ) over p gives a modified conservation equation for the longitudinal current, This reduces to Equation (32) as s − iω → 0. The perturbed current is conserved along each magnetic flux tube only in the absence of transverse gradients in ρ 0 or J x0 .

Perturbed Charge and Current Density
Here we evaluate the perturbed charge and longitudinal current densities, as given by Equation (31), using Equation (34) for f ± 1 . When evaluating the transverse gradient ∂f ± 0 /∂x, we choose n ± 0 to be constant (or, more precisely, assume a gradient scale much exceeding 2π/k x ). The p-integral of the term proportional to ∂f ± 0 /∂x is straightforward due to the narrowness of the momentum distribution, giving where In the case of a neutral plasma, the mode is purely growing or decaying (ω = 0) and this reduces to In what follows, we takeβ 0 to be independent of x, that is, we assume that the high-wavenumber magnetic flux component driving reconnection is superposed on a smoother component that contributes a greater fraction of the total current (Section 2.3). These equations for ρ 1 and J z1 can be further simplified by a change of variables. The presence of a strong guide field allows the definition a hydromagnetic displacement field ξ ⊥ = (ξ x , ξ y ) and a local 'rest' frame in which The first term in the expression for A z1 represents the effect of a hydromagnetic displacement. This is the magnetic perturbation associated with the second term in E in Equation (22). This displacement scales as ξ x ∝ x near x = 0, whereas the magnetic potential perturbation A z1 (0) remains finite in a tearing instability. Therefore, the non-ideal displacement ξ rec vanishes only in the ideal regime.
To evaluate the mode spectrum, we must iteratively compute the profiles of A z1 and φ 1 , or equivalently of ξ x and ξ rec x . Coulomb's law implies Substituting Equation (41) gives 3 Similarly, Ampere's law These equations simplify further when the current profile is harmonic in coordinate x. One sees that the coupling between the electrostatic and vector potential perturbations is driven by the non-ideal part of the displacement field.

NEUTRAL PAIR PLASMA WITH A STRONG GUIDE MAGNETIC FIELD
This section is devoted to an evaluation of the mode spectrum, in the case of a sinusoidal current profile (Equation (21)). This toy problem encapsulates the effects of multiple interacting tearing surfaces, where the non-potential magnetic field vanishes, The resulting eigenvalue problem is tractable. We first consider the simplest case of a charge-neutral plasma, n + 0 = n − 0 , with positive and negative charges counterstreaming, β − z0 = −β + z0 . Then the tearing modes are purely growing or decaying (ω = 0); we find multiple branches at k p0 /k x 1. These results are generalized in Section 7 to a non-neutral plasma (where ω = 0) and in Section 6 to a single tearing surface.

Growth Rates
A wide spectrum of modes is present when k p0 k x . Figures 2 and 3 respectively show the dependence of growth rate on (i) the ratio k p0 /k x of current sheet spacing to skin depth and (ii) on the mode wavenumber k y . The peak growth rate is, when k p0 k x , In the open pulsar circuit, where J z0 ∼ ρ co c = ΩB z0 /2π, the peak growth rate is therefore comparable to the angular frequency of rotation of the star, Faster growth is seen when the separation between tearing surfaces is smaller than the skin depth, k x > k p0 . Near a rotating neutron star, the tearing modes uncovered here grow slowly compared with the Alfvén frequency c/R NS , where R NS is the stellar radius. Furthermore, because s/ck is never a large number, no part of the eigenmode can be understand as evolving in the ideal regime. The non-ideal displacement ξ rec x dominates the hydromagnetic displacement ξ x not only near the surfaces B y0 = 0 (where A z1 is finite and so ξ rec x ∝ B −1 y0 ∝ x −1 ), but everywhere within the bulk. In other words, these modes cannot be divided into a narrow layer where flux freezing breaks down, surrounded by a more extended ideal magnetohydrodynamic zone, as in the standard tearing mode analysis of a collisionless plasma (Drake & Lee 1977;Basu & Coppi 1981;Ottaviani & Porcelli 1993). Right panel: Gradient scale l ξ of the electrostatic perturbation near x = 0, defined by ξx = (x/l ξ )Az1(0)/B ⊥0 . New branches of the dispersion curve appear as the tearing surfaces move apart relative to the skin depth (kp0/kx increases). These higher-order modes contain additional nodes in ξ rec x ; their growth rates converge on s ∼ 1, suggesting that a flux variation of low kx triggers a broad spectrum of secondary variations with wavenumber ranging from ∼ kx to ∼ kp0. Colors indicate the number Nrec of nodes in ξ rec x , ranging from Nrec = 0 (black) to Nrec 15 (dark red). Dispersion curves show a sudden shift in Nrec at high kp0/kx, marking a singularity in the electric field at x = 0 (where |l ξ | −1 → ∞ and Ex1(0) flips sign). In the right panel, square points denote l ξ > 0 and crosses l ξ < 0.  In all panels, the mode wavenumber ky = 0.3kx and dotted curves show perturbations with a reversed sign, e.g. Az1 < 0 in the top-left panel. As kp0/kx increases, the mode becomes strongly concentrated in the bulk, away from the tearing surfaces. Top-right panel: Az1 but now normalized to its peak value at kxx = ±π/2. Red-dotted lines shown the analytic approximation (55), which applies when kp0 kx. Red-dashed curve shows for comparison the profile (57) that would obtain if the bulk zone were fully in the ideal regime, with s → 0. Bottom-left panel: hydromagnetic x-displacement ξx, which is proportional to the scalar potential φ1 (Equation (41)). Bottom-right panel: electric field component Ex1 = −∂xφ1 transverse to the By = 0 surfaces; labelling of curves follows the other panels. Note that Az1 and Ex are symmetric about x = 0, whereas ξx and φ1 are anti-symmetric.  Figure 4. ξx always vanishes at the tearing surfaces By0 = 0 (x = 0, ±nπ) and ξy vanishes at the midpoint between these surfaces. As the system grows in size compared with the skin depth (kp0/kx increases) the displacement field is increasingly aligned with the tearing surfaces (|ξy| |ξx|, excepting at the midpoint. Several branches are apparent in the mode spectrum. Each branch is characterized by the number N rec nodes 4 in the profile of the non-ideal displacement field ξ rec x , as defined in Equation (41), over the interval 0 < x < π/2k x . The mode with the highest growth rate s has a profile with the smallest number of nodes (0). An increment in N rec is seen on each branch at a large, discrete value of k p0 . Here, the gradient of ξ x diverges near x = 0, implying a singularity in the transverse electric field -see the right panels of Figures 2 and 3.

Eigenmodes
The profile of the magnetic potential perturbation A z1 is shown in the top two panels of Figure 4, for several values of k p0 /k x (ranging from 10 −1 to 10 1.5 ) and mode wavenumber k y = 0.3k x . The main trend is a growing concentration of the mode displacement away from the surfaces B y0 = 0 as the system size increases with respect to the skin depth. For example, when k p0 /k x = 10 1.5 , the magnetic potential perturbation is ∼ 1 × 10 6 times larger at the mid-point (k x x = π/2) than it is at x = 0, ±π.
This result is easily understood quantitatively. When k p0 k x and |ξ x |, |ξ rec x | are large compared with A z1 (0)/B ⊥0 , Equations (48) and (49) combine to give where s = s ·β 0 ck y (B ⊥0 /B z0 ). Writing A z1 = Ce g(x) and approximating s,β 0 1, one has to leading order This formula (shown as dotted curves in the top-right panel of Figure 4) provides an excellent fit to the full solution for k p0 /k x = 10, 10 1.5 , and still a good fit for k p0 /k x = 10 0.5 . A useful comparison can be made with internal tearing in a collisionless but sub-relativistic plasma (Ottaviani & Porcelli 1993). Here the current system is nearly quasi-static throughout the bulk, corresponding to s 1. Then Equation (49) reduces to the equilibrium given by Equation (32), which may be written The solution which is symmetric about x = π/2k x is where κ = (k 2 x − k 2 y ) 1/2 . This is plotted as the dashed red line in Figure 4, showing that the mode is concentrated in the bulk, but not as sharply as the non-ideal modes obtained above.
We also obtain a simple derivation of the growth rate s 1 when k p0 /k x > 1, To see this, note that A z1 grows exponentially away from the surface B y0 = 0 but still satisfies the boundary condition at the mid-point, and Equation (58) follows by setting sin 2 (k x x) = 1 in Equation (54).
Other details of the eigenmodes are shown in the bottom two panels of Figure 4 and Figure 5: the x-and y−displacements ξ x,y and the transverse electric field. This provides an alternative view of how the displacement is concentrated away from the surfaces where B y0 = 0 when these surfaces are well separated compared with the skin depth.

ISOLATED CURRENT SHEET
We now shift to consider tearing at an isolated current sheet in a quantizing magnetic field. This configuration supports a more restricted set of tearing modes than does the harmonic current distribution with multiple tearing surfaces that we analyzed in Section 5. Growing modes are found only when the current sheet is relatively narrow; the overtone modes found previously for k x k p0 are absent when the mode is constrained to have a finite energy. We also note that our kinetic analysis does not uncover a separation between a non-ideal sublayer and a more extended, nearly ideal displacement of the magnetic field -in contrast with a previous analysis based on a resistive formulation Figure 6. Growth rate of tearing mode at an isolated current sheet, whose magnetic field profile is given in Equation (59), for three values of the ratio kp0/kx of current sheet thickness to skin depth. Mode has a finite energy, with Az1 and φ1 (or equivalently ξx and ξ rec x ) both decaying exponentially in the exterior zone with uniform By0. The growth rate drops toward the minimum smin given by Equation (61) as kp0/kx approaches a maximum value 1.86. Red dotted curve shows smin for kp0 = 1.5kx; the growth rate (black curve) drops to this value at ky 0.54kx and growth shuts off at larger ky. Modes with s > 0 also disappear when ky > kp0/β0; this is plotted as the blue dotted line in the case kp0 = 10 −1/2 kx. of force-free electrodynamics (Lyutikov 2003). In particular, the standard approximation of a nearly uniform magnetic flux perturbation in a resistive sublayer cannot be made in this situation. The flux perturbation retains a significant non-ideal component outside the current sheet; instability to tearing cannot be expressed in a simple way in terms of the jump in ∂ log A z1 /∂x across a resistive sublayer (e.g. Biskamp 2005). We have checked that attempting to enforce constant A z1 nearl x = 0 does not generate self-consistent, finite-energy tearing modes.
The isolated current sheet is taken to have the same current profile as previously (Equation (21)), but now restricted to a single reversal of B y0 , i.e., The perturbation equations (48) and (49) still apply in the current-free zones after taking k x → 0. Then the eigenfunctions for ξ x , ξ rec x can be written in the form Ae ikexx + Be −ikexx , where and s is defined after Equation (54). The behavior may be oscillatory or exponential, but the mode energy converges only when the perturbation is exponentially decaying at |x| > π/2k x . This constraint implies a minimum growth rate, which is obtained by setting k ex = 0 in Equation (60), Figure 7. Dependence of the growth rate of a tearing mode localized near a single current sheet on the ratio kp0/kx of current sheet thickness to skin depth. Curves correspond to the regime of small mode wavenumber ky, with growth possible for kp0 >β0ky. Growth rate diverges for small current sheet thickness, approximately as s ∝ (kx/kp0) 1/2 for kp0/kx < 0.1.
where ε k ≡ k 2 y /k 2 p0 . One finds real s min only for k y > k p0 /β 0 . The growth rates of these finite-energy modes are plotted in Figure 6 for a few ratios of current sheet thickness to skin depth. One finds that as k p0 /k x approaches a limiting value close to 2 ( 1.86 for k y /k x = 10 −2 ), then s → s min and k ex → 0. These results are obtained, as before, by iterating on s and dφ 1 /dx at x = 0, but now requiring that d ln A z1 /dx = d ln φ 1 /dx = −k ex sgn(x) at |x| = π/2k x .
The dependence of the growth rate on the thickness of the current sheet is shown in Figure 7. One sees that s diverges slowly, approximate as s ∝ k 1/2 x , as the current sheet grows thinner in comparison with the skin depth.
To summarize, we have demonstrated the existence of finite-energy and purely growing tearing perturbations, localized around an isolated current sheet in a relativistic, quantizing magnetic field, that have similar growth rates to the lowest-order modes found previously in the multiple-sheet configuration. These growing modes localized around a single current sheet exist only for a more restricted range of k p0 /k x and k y /k p0 than do the multiple-sheet modes. The current sheet thickness k −1 x must be smaller than or comparable to the skin depth, as must the reversal scale k −1 y of the excited mode along the tearing surface. The increase in growth rate with decreasing gradient scale k −1 x supports the formation of strong, localized current sheets during the non-linear development of the tearing mode, along the lines of the Syrovatskiǐ (1971) model.

NON-NEUTRAL PLASMA
We now turn to the more general case of a non-neutral plasma, writing The parameter ε ρ is a measure of the charge imbalance, since The quantity appearing in Equation (39) is now a complex number, meaning that we must take the frequency to have a real component, s → s − iω, and consider both A z1 and φ 1 to be complex numbers.
The problem simplifies when we consider ε ρ to be a small parameter and treat the effect of the charge imbalance as a perturbation to the charge-balanced solutions obtained in Section 5. Then the perturbations to A z1 and φ 1 are purely imaginary, δA z1 = iδA I z1 ; Figure 8. Real frequency induced in the dispersion relation of Figure 2 by a small charge imbalance in the plasma, as measured by the parameter ερ (Equation (62)). Frequency is negative for ερ > 0, and positive for ερ < 0. Growth rate s = s0 is unperturbed to first order in ερ. Results are obtained by iterating on solutions to Equations (69) and (70) so as to maintain the boundary conditions (71). Colors mark the same tower of modes as in Figure 2, with ξ rec x node number Nrec increasing from 0 (black) to 15 (dark red). and the real component of the growth rate perturbation vanishes to leading order, s 0 → s 0 − iω. To first order in ω and ε ρ , where is proportional to the small parameters ω/s 0 , ε ρ . We will also need where For the background current configuration, we return to the case of a sinusoidal profile (Equation (21)) and solve the mode Equations (44) and (45) expressed in terms of the displacement fields (41). Writing ξ x → ξ x + iδξ I x and ξ rec x → ξ rec x + iδξ rec,I x , and subtracting out the unperturbed eigenmode equations, we have and Here ξ rec x is obtained from the charge-neutral solution to Equations (48) and (49). Only the non-ideal part of the electromagnetic perturbation is responsible for driving the oscillation, just as it played a key role in driving the instability in the charge-neutral case.
We search for solutions to Equations (69) and (70) that have δA I z1 symmetric about x = 0, δξ x antisymmetric, and all perturbations vanishing at x = 0: The result is shown in Figures 8 and 9.

LINE TYING AND SURFACE PLASTIC FLOW
The tearing modes described in this paper involve a nearly incompressible and force-free rearrangement of a strong guide magnetic field. The transverse scale ∼ k −1 p,ex ∼ 0.1 − 100 cm of the mode is tiny compared with the length r c/Ω ∼ 5 × 10 8 P −1 cm of the sheared magnetic field. This field is embedded in the outer layers of a neutron star, whose shear strength varies greatly with depth. The transition from force-free to magnetoelastic equilibrium is spread over a distance somewhat greater than the magnetospheric skin depth.
Line tying can reduce the growth rate of a tearing mode (Huang & Zweibel 2009). We find that its effect is limited in the open circuit of a pulsar, given the extended length of the magnetic field lines. More significant is its effect on suppressing tearing modes in PIC simulations of pulsar magnetospheres, especially simulations focusing on the polar cap region. To avoid such an artificial suppression, a numerical model must simultaneously cover the entire magnetosphere as well as small-scale, tranverse structure in the open magnetic field bundle.
To understand the effect of line tying on reconnection in a force-free, relativistic plasma, consider first the case of low magnetization, σ 1 (Alfvén speed c). In the usual formulation of the problem, the sheared magnetic field is tied at both ends, with a length . The growth rate is reduced when k z ∼ −1 k ∼ k y B ⊥0 /B z0 ; then the mode becomes nearly force free and s ∝ (Huang & Zweibel 2009).
In the case of a quantizing magnetic field, as examined here, inertial effects are generally negligible and a growing tearing mode passes through a sequence of force-free equilibria. The magnetic field in the open pulsar circuit is tied only at one end, its length measured out to the light cylinder being comparable to c/Ω. The growth rate for infinite is s ∼ ck ∼ 2Ω(k y /k x ) (Equation (53)). Line tying will therefore have an effect on the growth rate when k y k x /2; the reduction in growth rate s → s · k for k < 1 then implies s → 4Ω(k y /k x ) 2 . High-order tearing of the slender open magnetic flux rope is insensitive to large-scale curvature of the rope (Section 2). The growth rate is also insensitive to magnetic flaring: when evaluated using a locally cartesian model, s is proportional to J z0 /B z0 , which is independent of radius inside the speed-of-light cylinder. Any curent irregularities formed at small radius will propagate along the open field bundle and be advected out into the pulsar wind.
The force-free structure of the tearing mode is modified by solid stresses only in a narrow layer below the stellar surface. The footpoints of the magnetic field are not fixed in position; instead the outer layers of the star flow plastically with them (Li & Beloborodov 2015;Thompson et al. 2017), down to a critical depth that we now derive. We start by balancing one component of the Maxwell stress with the yield stress in the solid, where ε y ∼ 10 −2 − 10 −1 is a temperature dependent yield strain (Chugunov & Horowitz 2010) and the shear modulus is µ 0.12(Ze) 2 (4πρ/3Am n ) 4/3 in a Coulomb solid composed of nuclei of charge Ze and mass Am n (Strohmayer et al. 1991). Then at the yielding depth, the crustal mass density is Here, we have set the background current to a fraction f J of ρ co c, that is, k x B y0 ∼ 2f J ΩB z0 /c, and normalized k x to the plasma scale using k 2 p = 4π(1 + 2M ± )|eρ co |/m e c 2 . The depth of the yielding layer can be compared to the magnetospheric skin depth. Using the hydrostatic relation between depth and Fermi energy, E F,e = (A/Z)m n g|z|, where g is surface gravity, approximating the degenerate electron gas as non-relativistic, and choosing a 26 Fe composition, one finds that the yielding depth is somewhat larger than the magnetospheric skin depth, k p |z| y = 2.8 (1 + 2M ± ) 1/4 B 5/4 z,12 There is a moderate enhancement in this ratio when M ± 1.

SUMMARY AND APPLICATION TO PULSARS AND MAGNETARS
We have identified and characterized internal tearing modes in a quantizing and nearly force-free magnetic field. The electric current profile contains both a smooth component and a high-wavenumber component that drives internal tearing. The embedded plasma is collisionless and relativistic, experiencing negligible curvature drift. Electron inertia plays a dominant role in the breakdown of magnetic flux freezing, although collisions are in practice expected to have non-negligible effects in magnetar magnetospheres. Our results form the basis for a 'top-down' approach to pulsar radio emission, in which slow tearing operating away from the magnetic separatrix are the ultimate driver of the emissionas opposed to fast plasma oscillations that are excited in a quasi-electrostatic gap (Levinson et al. 2005;Philippov et al. 2020).
Our main findings are as follows. 1. A robust but slow instability is uncovered in the presence of multiple, interacting tearing surfaces. The growth rate is suppressed by the inverse of the strong guide magnetic field running perpendicular to the plane of reconnection. The dispersion relation is derived for a wide range of separations between tearing surfaces, measured relative to the skin depth. When this separation is large, the peak growth rate is very nearly s 0 ≡ 4π(k y /k x )J z0 /B z0 , where J z0 is the high-wavenumber component of the current density and B z0 the guide field ( Figure 2 and Equation (52)). Faster growth is obtained when the background current is irregular on a scale smaller than the skin depth k −1 p0 ∼ γ −3 0 −1/2 (mc 2 /4πn 0 q 2 ) 1/2 . We find that current irregularities are susceptible to tearing over a wide range of scales. 2. Overstable modes are uncovered when the plasma carries net charge, as is expected in the magnetosphere of a radio pulsar. The real frequency ω is comparable to the growth rate s times the relative imbalance |n + 0 − n − 0 |/|n + 0 + n − 0 | in the densities of positive and negative charges (Figure 8 and 9).
3. An isolated current sheet whose thickness is smaller than or comparable to the skin depth (k p0 /k x 2) is shown to sustain localized tearing modes (Figures 6 and 7). Mode growth remains slow but scales with the ratio k p0 /k x of current sheet thickness to skin depth as s ∼ (k x /k p0 ) 1/2 s 0 for large k x /k p0 . This opens up a mechanism for runaway collapse of the current sheet even in the presence of a strong guide field. 4. Overtone tearing modes, showing one or more nodes in the non-ideal displacement field ξ rec x (x), can be excited when neighboring tearing surfaces are separated by a distance exceeding the skin depth, k x < k p0 (Figure 2). The excitation of these overtone modes opens up a mechanism for creating small-scale structure in the current surrounding a smooth extremum in the twist profile, as is needed to redistribute twist across a magnetic flux bundle. 5. Line tying has a limited effect on the growth of small-scale tearing modes in the open circuit of a rotating neutron star. In significant part, that is because the growth length c/s is comparable to the length of the open-field bundle, but also because inertial effects have a negligible effect on mode growth even in the absence of line tying. We infer a smooth transition from the force-free state of the magnetosphere to a weak magnetoelastic deformation of the upper solid crust. This transformation is concentrated at a depth somewhat larger than the transverse size of the magnetospheric current irregularities (the magnetospheric skin depth; Equation (74)).
6. It nonetheless should be emphasized that line tying will have the practical effect of suppressing tearing mode growth in a PIC simulation of the open pulsar circuit that covers only a limited range ∆r of radius near the surface of the star: growth rates will be suppressed by a factor ∼ Ω∆r/c (see Huang & Zweibel 2009 for the case of non-relativistic reconnection). 9.1. Peculiarities of Tearing Modes in a Quantizing Magnetic Field 1. The dynamics can be understood in terms of charge flow along separate magnetic flux tubes, which interact via the Coulomb and Lorentz forces. Charge density and current density perturbations are of comparable importance in driving instability. The vorticity field that is a central actor in non-relativistic reconnection is replaced by the charge density field (Equations (16) and (17)).
2. The modes uncovered here are a relativistic generalization of the double tearing mode, but with the distinction that the non-ideal and ideal displacements are comparable in magnitude throughout the bulk. There is no separation between a narrow non-ideal layer and a larger, hydromagnetic deformation of the magnetic field. Growth rates are obtained for separations of the tearing surfaces up to ∼ 30 times the skin depth. In the perturbation equations, the non-ideal displacement (see Equation (41)) is seen to be the driver of mode growth (Equations (44) and (45)). Similarly to the double tearing mode (Bierwage et al. 2005;Wang et al. 2007), the displacement becomes very strongly peaked in the bulk, as the surfaces B y0 = 0 become separated by more than a few skin depths (Figure 4).
3. Also related to the non-ideal nature of the tearing instability is the failure of the 'constant-Ψ' approximation, 5 as applied to a single tearing surface in quantizing magnetic field. In this approximation, the energy of the perturbation diverges with distance from the current sheet. This divergence is traced to the appearance of coordinated current and charge oscillations.
4. The tearing mode with fastest growth has the fewest number of nodes N rec in the profile of ξ rec x transverse to the tearing surface. When k p0 k x (the tearing surfaces are separated by a distance much greater than the skin depth), a large tower of independent branches appears in the dispersion relation, labeled by N rec > 0. The wavenumber of the x-oscillation lies in the range k x < k < k p0 . The growth rates of these overtone modes are nearly as large as the growth rate of the fundamental mode with zero nodes (Figure 2). 5. Nonlinear effects associated with magnetic island saturation (Rutherford 1973) are expected to have a limited effect on mode growth, given that the displacement is concentrated away from the tearing surface when k p0 k x .

Future Directions
The nonlinear development of tearing in a quantizing magnetic field is left to future work. Existing time-dependent calculations of tearing in a non-relativistic plasma show strong nonlinear growth when current sheets of opposing signs are in close proximity (Ishii et al. 2002;Bierwage et al. 2005;Zhang & Ma 2011). This growth is driven by the expansion and eventual interaction of distinct populations of magnetic islands. Thus, the energy released by the instability described here may greatly exceed the simplest estimate extracted from linear theory. Current irregularities forming near the skin depth also fit within a hierarchy of magnetic structures. For example, the energy carried by a tearing mode feeding off a current irregularity of wavenumber k x will scale as k −1 x B 2 y0 ∝ k −3 x , given that the seed current irregularity J z0 ∼ ck x B y0 /4π is normalized by the mean current density ∼ ρ co c in the pulsar circuit. The fastgrowing overtone modes that are found for k x k p0 give an indication of how longer-wavelength current irregularities may spawn smaller-scale magnetic structures.
Another avenue for further work involves generalizing the tearing dispersion relation derived here (for a narrow top-hat particle distribution function, with counterstreaming positive and negative charges) to particle distributions with (i) a broader spread in momentum and (ii) a collective flow of positive and negative charges (as expected in the parts of the pulsar polar cap experiencing a pair cascade).
Heavy ions will also have an effect on the mode structure. A space density ∼ ρ co /|Ze| is expected when the corotation charge is positive. In the absence of pairs, the characteristic scale of the tearing eigenmodes is increased by a factor ∼ (Am p /m e ) 1/2 (here Z and A are the ion charge and mass in atomic units) but the growth rate and drift rates are otherwise independent of the mass of the particle carrying the corotation charge. In the presence of pairs, the oscillatory behavior of the overtone tearing eigenmodes will be modified compared with the case of negative corotation charge.
9.3. Application to Radio Pulsars 1. Slow mode growth and pulsar flux variability. A role for the instability uncovered here in pulsar radiation emission is motivated by the presence of a strongly inhomogeneous current profile in the open magnetic flux bundle of a rotating neutron star (Timokhin & Arons 2013;Gralla et al. 2017). In many cases, the bundle shows non-axisymmetric features and an inhomogeneous distribution of magnetic twist that is known to trigger fast resistive instabilities in tokamaks (e.g. the double tearing mode) and to induce a rapid relaxation in the twist profile toward a local maximum (White 2013). The vigor of such an instability is limited in the pulsar context by the presence of a strong guide field: the peak growth rate is ∼ Ω, the angular rotation frequency of the star, which is much slower than the characteristic frequencies of particle and field oscillations in a polar gap. This is slow enough to be relevant to the dramatic variations that are observed in radio emission over successive rotations (Graham-Smith 2003).
Establishing a connection between magnetospheric tearing and radio emission depends on identifying a secondary instability that can induce charge clumping in current sheets. Paper II describes an overstability of propagating charge perturbations that are localized near Debye-scale current sheets and are excited by a Cerenkov-like process.
2. Sub-pulse drift. When the plasma carries net charge, as it must in the pulsar magnetosphere, we uncover overstable tearing modes with a finite real frequency. In a nearly axisymmetric configuration, as expected in a pulsar with small inclination between the magnetic dipole and rotation axes, this overstability corresponds to a general drift of the tearing modes in the azimuthal angle around the polar cap. There is a promising connection here with the secular 'sub-pulse' drift that is observed in many pulsars (Deshpande & Rankin 2001). Indeed, calculations of high-order (double, triple...) tearing modes in non-relativistic plasmas driven by closely spaced tearing surfaces (B y0 = 0) show such an angular structure (Bierwage et al. 2005).
The magnitude of the drift rate depends on the particle density. It is comparable to observed drift rates in parts of the magnetosphere that have a low, or vanishing, e ± pair density in the absence of tearing (Paper II).
3. Tearing at high and low current density. The slow tearing instability described here does not require that the background current density J z exceed the limiting corotation charge flow ρ co c. Self-consistent models of this particle flow indicate that the charges begin to stream relativistically as J z approaches ρ co c (Timokhin & Arons 2013). The nonlinear development of slow tearing will increase the current density locally and, at a given rotation rate, expand the part of the polar cap experiencing e + -e − pair cascades. In the absence of pair creation, the pulsar magnetosphere collapses from the configuration postulated by Goldreich & Julian (1969) to a static, charge-separated structure without relativistic particle flow (Krause-Polstorff & Michel 1985;Chen & Beloborodov 2014;Philippov & Spitkovsky 2014).
Small-scale tearing may therefore supplement general-relativistic frame dragging (Muslimov & Tsygan 1992;Philippov et al. 2015) as a driver of charge starvation in the polar cap and of activity as a pulsing radio source.

Application to Magnetars
Neutron stars with magnetic fields strong enough to cause significant crustal yielding -magnetars -are frequently observed as intense sources of non-thermal X-ray, with a luminosity far exceeding what can be supplied by the measured loss of rotational energy (Kaspi & Beloborodov 2017). One infers the presence of strong currents flowing through the closed magnetosphere and, from the simultaneous measurement of thermal X-ray hotspots, a strong concentration of these currents on the magnetar surface.
The excitation of strong, localized magnetic shear during fault-like displacements of the crust has been proposed as a trigger for runaway pair creation in magnetar X-ray flares (Thompson & Duncan 2001;Thompson et al. 2017). A path toward understanding how thermalization might occur in the stellar magnetosphere is provided by the tearing instability described here, in combination with the secondary Cerenkov instability of shear Alfvén waves investigated in Paper II. High-frequency Alfvén waves that experience rapid, in situ damping may be excited directly in the magnetosphere, as opposed to the indirect excitation by crustal elastic modes (Blaes et al. 1989). In the case of the giant magnetar X-ray flares, crustal elastic modes are an inefficient source of magnetospheric dissipation (e.g. Thompson & Duncan 2001;Link 2014). Magnetar outbursts sometimes involve a brief non-thermal X-ray transient, followed by a slower relaxation lasting months (e.g. Woods et al. 2001Woods et al. , 2004. The slow component has been modeled in terms of ohmic decay of a localized magnetospheric current (Beloborodov 2009); alternatively, a global magneto-thermal-elastic model of the crust reveals continued creep at crustal faults following an outburst (Thompson et al. 2017), from which one infers a role for internal tearing in powering delayed persistent X-ray emission. The instabilities described here and in Paper II may also play a role in the microphysical development of a global reconnection event, involving the expulsion of a magnetic loop by the closed magnetosphere, during a giant magnetar flare (Lyutikov 2006).