Monster Radiative Shocks in the Perturbed Magnetospheres of Neutron Stars

Magnetospheres of neutron stars can be perturbed by star quakes, interaction in a binary system, or sudden collapse of the star. The perturbations are typically in the kilohertz band and excite magnetohydrodynamic waves. We show that compressive magnetospheric waves steepen into monster shocks, possibly the strongest shocks in the Universe. The shocks are radiative, i.e., the plasma energy is radiated before it crosses the shock. As the kilohertz wave with the radiative shock expands through the magnetosphere, it produces a bright X-ray burst. Then, it launches an approximately adiabatic blast wave, which will expand far from the neutron star. These results suggest a new mechanism for X-ray bursts from magnetars and support the connection of magnetar X-ray activity with fast radio bursts. Similar shocks may occur in magnetized neutron-star binaries before they merge, generating an X-ray precursor of the merger. Powerful radiative shocks are also predicted in the magnetosphere of a neutron star when it collapses into a black hole, producing a bright X-ray transient.

1. INTRODUCTION Magnetospheres of neutron stars are formed by plasma immersed in a strong magnetic field B bg .They have a high magnetization parameter, where ρ bg is the plasma mass density and c is the speed of light; subscript "bg" stands for "background" for waves investigated in this paper.The closed magnetosphere is approximately dipole at radii r much greater than the star radius R ⋆ .It ends and a magnetized wind begins near the light-cylinder R LC = c/Ω, where Ω is the star rotation rate.This basic picture is confirmed by extensive studies of pulsars (Philippov & Kramer 2022).

Perturbations of neutron star magnetospheres
In some cases, the magnetosphere becomes significantly perturbed.In particular, quakes in magnetars launch magnetospheric waves (Blaes et al. 1989;Thompson & Duncan 1996;Bransgrove et al. 2020).Some mechanism quickly dissipates the waves, generating Xray bursts that are observed as the main form of magnetar activity (Kaspi & Beloborodov 2017).Significant perturbations are also expected in tight neutron-star binaries, where the two magnetospheres interact with each amb@phys.columbia.eduother.A huge disturbance of a neutron star magnetosphere occurs when the star collapses into a black hole (Lehner et al. 2012).In all these cases, the excited waves are typically in the kHz band.
Such perturbations are well described by magnetohydrodynamics (MHD), which supports waves of the perturbed magnetic field and electric field E ⊥ B bg .These MHD modes can have wavevectors k in any direction and can be of two types: (1) shear Alfvén waves with k • E ̸ = 0, and (2) compressive waves with E ⊥ k, socalled "fast magnetosonic modes."Both modes propagate with an ultrarelativistic group speed, Phase speed v ph = ω/k of the magnetosonic waves is also close to c while Alfvén waves have v ph ≈ c cos α, where α is the angle between k and B bg .A detailed discussion of waves in e ± plasma in a strong background magnetic field B bg , including two-fluid and single-fluid (MHD) descriptions, is found in Arons & Barnard (1986).

MHD waves from quakes
Substantial work was devoted to the relativistic Alfvén waves excited by magnetar quakes (e.g.Thompson & Blaes 1998;Troischt & Thompson 2004;Yuan et al. 2020;Nättilä & Beloborodov 2022), and little attention was given to the magnetosonic waves.In fact, quakes can excite both modes.Neutron star quakes involve shear oscillations of the crust with horizontal displacements δr.Shear waves inside the neutron star crust propagate with speed v sh = ω/k sh ≈ 10 8 cm/s (Blaes et al. 1989).Their characteristic lowest frequency is ω ∼ v sh /h ∼ 10 4 rad/s, where h is the hydrostatic scale height.A large fraction of the quake energy may be at much higher frequencies, e.g.ω ∼ 10 5 rad/s.The crustal waves leak to the magnetosphere with a transmission coefficient T ∼ 0.1 ω 3/5 5 B 2/5 15 (Blaes et al. 1989;Bransgrove et al. 2020).
To quickly see that quakes can emit both Alfvén and magnetosonic modes, one can consider crustal oscillations with ω ≫ c/R ⋆ and a wavevector that is exactly radial in a region of the stellar surface (Figure 1).Such oscillations satisfy ∂ ϕ = ∂ θ = 0 in spherical coordinates r, θ, ϕ, and this symmetry is preserved during wave transmission, so the emitted magnetospheric wave also has a radial k.This simple setup is reduced to a locally plane-parallel transmission problem; it was used by Blaes et al. (1989) to study how quakes with δr ∥ k × B bg excite Alfvén waves in the magnetosphere.It is easy to see that quakes with δr ∦ k × B bg emit magnetosonic modes.Polarization of the emitted wave is set by the electric field in the crustal oscillation, E = B bg × v/c (where v = −iωδr), and its type is determined by the orientation of E relative to the vector n ≡ k × B bg .It is a pure Alfvén mode if E ⊥ n and a pure magnetosonic wave if E ∥ n.The condition for magnetosonic excitation E • n ̸ = 0 may be stated as The transmitted quake power is partitioned between the Alfvén and magnetosonic modes as E 2 A /E 2 ms = tan 2 ψ where ψ is the angle between the quake motion v and the horizontal component of the magnetic field B h bg .
A complete picture of wave emission is complicated by several additional effects: (1) Besides the vertical (radial) wavevector k r , crustal waves have a horizontal component k h (Bransgrove et al. 2020) and excite non-radial magnetospheric waves.As an example, Yuan et al. (2020) examined wave generation in a uniform oblique B bg attached to a horizontal conducting surface perturbed by axisymmetric horizontal shear.They showed that the shear excites a mixture of Alfvén and magnetosonic waves.These calculations should be extended in the future to include a realistic density profile of the crust.
(2) Waves of lower frequencies ω < ∼ c/R ⋆ develop at larger radii R w ∼ R ⋆ + c/ω.In the region R ⋆ < r < R w the magnetosphere adjusts to the surface oscillation in a quasi-static manner and, effectively, the magnetospheric footpoints in the wave are relocated from the stellar surface to the sphere of radius r ≈ R w .The quasi-static deformation of the magnetosphere at r < ∼ R w differs from the crustal deformation and needs further investigation of its compressive component that could drive a smallamplitude magnetosonic wave.For example, a strongly twisted magnetosphere in axisymmetric equilibrium responds to additional surface shear ∂ θ v ϕ ̸ = 0 by inflating at r < ∼ R w (Parfrey et al. 2013).
(3) The twisted magnetospheres near magnetars have significant spatial variations of the toroidal magnetic field component B ϕ , so vector n = k × B bg can change its direction on a scale comparable to r.Then, any attempt to launch a pure Alfvén mode at r ∼ R w ∼ c/ω (for ω < ∼ c/R ⋆ ) will inevitably generate a mixture of the Alfvén and magnetosonic waves because the two linear modes have different refraction indices c/v ph and will be unable to adiabatically track the changing local n until propagating to r ≫ c/ω.
Pure Alfvénic excitations occur in the simple case of axisymmetric quakes with azimuthal v and an untwisted background magnetosphere (then k ϕ = 0 and B ϕ bg = 0, so n = k × B bg ∥ v).Even in this case, the emitted Alfvén waves can convert to magnetosonic waves because of nonlinear effects.Alfvén waves have v gr ∥ B bg , so they are ducted along the closed magnetic loops, and the nonlinear conversion peaks when the wave reaches the top of the loop (Yuan et al. 2021).In addition, at later times, the bouncing Alfvén waves trapped in the loop develop a nonlinear turbulent cascade which emits magnetosonic waves (e.g.Li et al. 2019).
Besides sudden quakes, the crust may flow plastically and slowly twist a magnetospheric flux bundle to an instability threshold.Three-dimensional simulations of this process demonstrate that relaxation of the unstable flux bundle generates fast magnetosonic waves in the magnetosphere (Mahlmann et al. 2023).

Shocks
In the limit of high σ bg , the emitted magnetosonic waves are equivalent to electromagnetic radio waves propagating without coupling to the magnetized plasma: the oscillating E ⊥ B bg drives a tiny electric current, negligible compared to the displacement current ∂ t E/4π.Therefore, compressive perturbations of the magnetosphere are usually assumed to propagate as vacuum electromagnetic waves superimposed linearly with the background field B bg , and freely escape. 1 The linear propagation would imply no shock formation.
However, this simple picture is safe only for lowamplitude waves, |E| ≪ B bg , i.e. near the neutron star where B bg is strong.As the wave expands to larger radii r, the dipole background field decreases as B bg ∝ r −3 while the wave amplitude E 0 decreases as r −1 .Their ratio E 0 /B bg grows as r 2 and eventually approaches unity.It is easy to see that the linear propagation of an oscillating wave becomes impossible when where α is the angle between B bg and the wavevector k.Indeed, the linear superposition of B bg with the vacuum electromagnetic wave (E w = E and where we used E 2 w = B 2 w .MHD description breaks if B 2 < E 2 .This occurs when the condition (4) is met.
For instance, consider the equatorial plane of a dipole magnetosphere.Here, a spherical wave front propagates with k ⊥ B bg (α = π/2) and the magnetic field perturbation B w oscillates along B bg (Figure 2).The linear superposition wave+background gives minimum Note that the plasma oscillates with the drift speed v/c = E × B/B 2 , and E 2 → B 2 corresponds to |v| → c.This implies a runaway growth of the plasma kinetic energy, and the MHD wave can no longer be approximated as a vacuum electromagnetic wave.Such energy conversion is a general effect of E 2 → B 2 (it also happens in Alfvén waves (Li et al. 2019(Li et al. , 2021)); see Levinson (2022) for a recent discussion).In magnetosonic waves, energy conversion prevents reaching E 2 = B 2 by steepening the wave into a shock.The appearance of shocks at E 0 = B bg /2 was previously noted in the context of waves in pulsar winds (Lyubarsky 2003).A fully kinetic (local-box) particlein-cell simulation of a magnetosonic wave propagating through a decreasing B bg has recently been performed by Chen et al. (2022).Their simulation demonstrates the sudden steepening of the wave into a shock when Waves consisting of half oscillation with B w • B bg > 0 and no part with B w • B bg < 0 would never face the E 2 = B 2 limit.One may think that in this case the wave avoids shock formation.However, as explained below, such half-waves also steepen into shocks, although this occurs gradually, at a larger radius.This gradual steepening creates a forward shock, launching an accelerating blast wave that expands far beyond the magnetosphere.

This paper
Since shocks appear to be a generic outcome of magnetosonic perturbations, a few questions arise: How strong are the shocks?What fraction of the original wave en-ergy gets dissipated by the shock?What fraction of the dissipated energy is radiated?
These questions can be answered by solving a welldefined MHD problem with a simple initial condition: launch a spherically expanding magnetosonic wave with an initial amplitude E 0 /B bg ≪ 1 so that initially (at small radii) it behaves as a vacuum electromagnetic wave.Then, track its expansion through the dipole magnetosphere and examine how it steepens into a shock and propagates afterward.This problem is solved in the present paper.We will show that the plasma Lorentz factor γ in the magnetospheric shocks caused by E 2 → B 2 reaches huge values γ ∝ σ bg , likely making them the strongest shocks in the universe.Therefore, we call them "monster shocks."They differ from normal collisionless plasma shocks because they are highly radiative: we will show that the plasma approaching the monster shock radiates its kinetic energy before forming the downstream flow.
Tracking the evolution of magnetospheric waves with shock formation presents an interesting technical challenge.The magnetization parameter σ bg at radii of main interest can exceed 10 10 (Section 2.1).The large σ bg is usually replaced by σ bg → ∞, which corresponds to taking the force-free electrodynamics (FFE) limit of MHD.However, FFE cannot describe shock formation.FFE neglects plasma inertia, so its wave modes propagate with exactly speed of light and cannot steepen.
Since FFE does not provide a suitable framework, one has to examine the problem using the full MHD equations.The equations state conservation of energy and momentum and could be solved numerically with customary discretization methods.However, in practice such methods fail at high magnetization σ (most existing MHD codes have to keep σ < 100 to avoid numerical issues).We take a different approach: we solve the MHD equations along characteristics.It provides an efficient method for both finding and understanding the solution at arbitrarily high σ bg .The solution is easiest to find for short waves, with wavelength λ ≪ r.
Calculations in this paper are performed for axisymmetric wave packets (∂ ϕ = 0).Setting ∂ ϕ = 0 should also be a good approximation more generally for waves far from their source, where wavevectors k are radial.The wave evolution along the radial ray is controlled by the local B bg encountered along the ray; we consider the simple case of a dipole B bg .
Formulation of the problem and our method for solving it are described in Section 2. Section 3 explains how we track shocks in MHD waves and then Section 4 presents the full numerical simulation, performed for an equatorial wave in a dipole magnetosphere.Remarkably, the problem also admits a complete analytical description, which is given in Section 5. Emission from the monster shock is discussed in Section 6.The results are discussed in Section 7.
2. NONLINEAR MAGNETOSONIC WAVES MHD fluid is described by the plasma mass density ρ, velocity v = cβ, magnetic field B, and electric field E. We wish to investigate magnetosonic waves launched in a dipole magnetosphere and find their nonlinear evolution.Term "nonlinear" here means that (a) the wave oscillation is not necessarily small compared to the background field B bg , and (b) the wave can be strongly deformed during its evolution and form shocks.

Unperturbed background magnetosphere
The wave will be shown to experience strong evolution at radii r much greater than the neutron star radius R ⋆ ∼ 10 6 cm, but well inside the light cylinder R LC (all known magnetars in our galaxy have R LC > ∼ 10 10 cm).Our calculations below use the background at radii R ⋆ ≪ r ≪ R LC , the magnetosphere is approximated as dipole, and its rotation is neglected.Thus, for simplicity, we neglect any twists of the outer magnetosphere, i.e. assume B ϕ bg ≈ 0 at radii of interest.The electromagnetic field of a static dipole magnetosphere is where (e r , e θ , e ϕ ) is the normalized basis in spherical coordinates r, θ, ϕ with the polar axis along the magnetic dipole moment µ.Magnetars have µ ∼ 10 33 G cm 3 .Plasma in the unperturbed static magnetosphere will be approximated by where the dimensionless parameter N ≡ n bg r 3 is approximately constant with r (Beloborodov 2020); m is the particle (e ± ) mass, and n bg = ρ bg /m.The lowest N ∼ 10 31 µ 33 (Ω/rad s −1 ) is set by the Goldreich-Julian density.The more typical N for magnetars is much higher, N ∼ 10 37 , due to strong electric currents near the star accompanied by e ± pair creation with a high multiplicity (Beloborodov 2021a). 2 The magnetization parameter of the background plasma The magnetosphere at r > 10 7 cm is populated with mildly relativistic e ± , as they are decelerated by drag 2 The created pairs flow along the magnetic field lines with a speed controlled by the magnetar's radiation field (Beloborodov 2013).
In the outer closed magnetosphere, the plasma accumulates and annihilates at θ ≈ π/2, in a layer of high e ± density controlled by annihilation balance and the thickness of the annihilation layer.The e ± outflow terminated by the annihilation sink at θ = π/2 may have N ∼ 10 37 at all θ except the thin equatorial layer.
We neglect these background motions; thus, the plasma energy density (including rest mass) in the unperturbed background ahead of the wave is approximated as ρ bg c 2 .

Axisymmetric waves
We will examine the radial expansion of an axisymmetric magnetosonic wave.Its electric field E is perpendicular to both B bg and the (radial) propagation direction, so vector E oscillates along ϕ.It is convenient to define scalar E by The wave magnetic field is perpendicular to E and oscillates in the poloidal plane, Any axisymmetric wave can be described by a toroidal electromagnetic potential, This class of fields includes the static dipole background as a special case, with Charge density remains zero in magnetosonic waves, 4πρ e = ∇ • E = 0. Electric current density j satisfies the Maxwell equation 4πj = −∂ t E + c ∇ × B. In axisymmetric waves the current is toroidal: Note also that E are B are related by induction equation As a concrete example, we will consider a wave launched with an initial sine profile One full oscillation corresponds to 0 < ξ < λ/c = 2π/ω.

Short waves (λ ≪ r)
The condition λ ≪ r considerably simplifies the wave propagation problem.In particular, it leads to a simple expression for current j, which will be used in Section 2.4 to describe the wave-plasma interaction.
Let us define the wave potential A w ≡ (A − A bg )e ϕ .It determines the wave fields and j: For short waves, it is helpful to view fields as functions of t, ξ, θ instead of t, r, θ.The fast oscillation then becomes isolated in the coordinate ξ.Differential equations can be rewritten in variables t, ξ, θ using Using the θ component of induction equation ( 14), we rewrite Equation (17) as In short waves, the oscillation of A w with ξ is much faster than its variation with t or θ at fixed ξ.Hence, as B r w does not contain the large derivative ∂ ξ (rA w ).Note also that rE = ∂ ξ (rA w ) t + ∂ t (rA w ) ξ , and Therefore, the expression for j simplifies to There is also a useful relation between plasma density n = ρ/m and velocity v in short waves.In general, n and v satisfy the continuity equation ∂ t n + ∇ • (nv) = 0.In the short-wave limit, it simplifies to F ξ is the particle flux through the surface of ξ = const.
It is uniform across the wave packet, and evolves as the packet propagates to larger radii: F ξ = cn bg ∝ r −3 .

Wave-plasma interaction
We wish to find the evolution of the wave profile E(ξ).It evolves because the electromagnetic field exchanges energy and momentum with the plasma oscillating in the wave.The plasma motion is described by four-velocity, where β = v/c.The equation of motion is with the derivative d/dt taken along the fluid streamline: d/dt = ∂ t +v•∇.Taking the scalar product of both sides with β, and using E + β × B = 0, one obtains This equation expresses energy conservation.Note that using ρc 2 = nmc 2 in Equation ( 26) assumes a negligible contribution of internal (thermal) energy to the plasma inertial mass.As explained below, this is a reasonable approximation, because significant wave-plasma interaction happens where the plasma has a low temperature and a high γ.
Substitution of Equation ( 23) into Equation (27) gives where we used dξ = dt − dr/c = (1 − β r )dt along the fluid streamline and ρ(1 − β r ) = ρ bg (Equation 24).The derivative dγ/dξ written out in coordinates x α = (t, ξ, θ, ϕ) takes the form, One might expect the terms with ∂ t γ and ∂ θ γ to be small compared to ∂ ξ γ by the factor of ∼ λ/r ≪ 1.However, nonlinear evolution can make the term with ∂ t γ important.As shown below, this happens when γ 3 > ∼ σ bg .Now we can write the energy equation in its final form: where derivative ∂ t is taken at fixed ξ, θ, ϕ, i.e. along the radial ray r = ct + const.This equation connects the evolution of E(t, ξ, θ) and γ(t, ξ, θ).Equation ( 30) has a simple intuitive interpretation.Multiplication of both sides by an infinitesimal δξ gives Here, Ṅξ ≡ 4πr 2 F ξ is the isotropic equivalent of the particle flux F ξ , and δγ = (dγ/dξ)δξ is the change of γ along the fluid streamline as it crosses δξ; it determines the energy gain of the plasma crossing δξ per unit time, Ṅξ mc 2 δγ.The quantity δE = c r 2 E 2 δξ is the electromagnetic wave energy (isotropic equivalent) contained in the infinitesimal part δξ of the wave profile.This interpretation of δE is consistent with the energy density for the electromagnetic wave in the short-wave limit.
Thus, Equation (30) merely states that the wave profile r 2 E 2 (ξ) evolves by exchanging energy with the plasma.The plasma Lorentz factor γ(ξ) oscillates, and E(ξ) becomes deformed because the plasma passing through the wave receives energy at some ξ 1 and returns it to the electromagnetic field at another ξ 2 .E(ξ) is systematically reduced where dγ/dξ > 0 and increased where dγ/dξ < 0, leading to the steepening of E(ξ).
The FFE limit of MHD would correspond to setting ρ bg → 0, so that the r.h.s. of Equation ( 30) vanishes.Hence, the FFE solution is rE = const at ξ = const.This wave excites no electric current j, as can be verified using Equation ( 13), and so the FFE limit gives a simple vacuum wave superimposed on B bg .
Equation ( 30) is sufficient to find the evolution of MHD waves if they drive pure radial plasma motions, as happens in the equatorial waves described below.Then, energy conservation (or the radial component of momentum equation ( 26)) contains the complete dynamical information.Wave propagation outside the equatorial plane (θ ̸ = π/2) involves additional θ-motions governed by the θ component of Equation ( 26).

Equatorial waves
Main features of shock development will be shown by tracking waves at θ = π/2 that are symmetric about the equatorial plane.Symmetry implies u θ = 0 and B r = 0 at θ = π/2.The MHD fluid then oscillates with velocity We will use the following notation: These definitions imply while E, B, and β may be positive or negative.Note that Equation ( 30) for the equatorial wave becomes Recall that the partial derivatives ∂ t and ∂ ξ are defined in variables (t, ξ) and can act on any function E(t, ξ), γ(t, ξ), or r(t, ξ) = c(t − ξ) entering Equation ( 35).
The equatorial wave problem can be reduced to one unknown function of t, ξ.Indeed, any short wave packet propagating into an initially unperturbed plasma has a useful feature: fluid compression n/n bg at any point in the packet is related to the local speed β (Equation 24).The magnetic field is frozen in the fluid and compressed by the same factor, so The electric field E = βB can be expressed as Thus, all MHD quantities in a short wave are functions of β, and Equation ( 35) can be recast so that it contains only derivatives of γ.This requires expressing ∂ t (rE) in terms of ∂ t γ.Differentiating Equation (37), we find3 where we used ∂ t (rB bg ) ξ = c d(rB bg )/dr = −2cB bg and dγ = γ 3 β dβ.Then, substituting Equation (38) into Equation ( 35), we find Note that the term (1 − β) −1 in the bracket (brought by ∂ t γ on the r.h.s. of Equation ( 35) where ρ and B are measured in the fluid rest frame, and 39) rewritten in terms of κ becomes Here, in the coefficient of ∂ t κ we neglected κ 2 /2 ≪ 2σ bg κ 3 , using The equatorial wave evolution will be found if we solve Equation (41) for κ(t, ξ).This first-order partial differential equation is linear in the derivatives ∂ t κ and ∂ ξ κ, and can be solved using the method of characteristics.

Characteristics
One can rewrite Equation (41) as where the derivative dκ/dt is taken along curves C + (characteristics) defined by The characteristics ξ + (t) can also be described by their radial speed, In an accompanying paper we give an alternative derivation of the wave evolution equation generalized to relativistically hot plasmas, and discuss two families of MHD characteristics C ± propagating with radial speeds 4 Here, β s has the meaning of wave speed relative to the plasma (the "magnetosonic speed").In a magnetically dominated plasma β s is close to unity.In particular, for a cold plasma, The FFE limit corresponds to σ bg → ∞ and β s → 1.In this limit, Equations ( 43) and ( 44) simplify to dκ/d ln r = κ − κ −1 and dr + /dt = c (⇔ dξ + /dt = 0), and give the solution where

Shock formation
Significant MHD corrections to FFE arise if σ = κσ bg drops, i.e. if the plasma experiences strong expansion in the wave, κ ≪ 1.This occurs where the plasma is accelerated as E 2 approaches B 2 .Note that strong expansion also implies strong adiabatic cooling.Thus, shock formation is expected in the part of the wave with large γ, small κ, and a reduced temperature.The process of MHD shock formation is illustrated in Figure 3.The shock occurs where the characteristics cross, bringing different values of β to the same location and thus creating a discontinuity.
We consider waves emitted at sufficiently small radii where B bg far exceeds the wave electric field E, and the plasma oscillates with small |β| ≪ 1, which implies a negligible change in plasma density, |κ − 1| ≪ 1.At the small radii, characteristics propagate with speed β + = 1 − O(σ −1 bg ) and FFE is an excellent approximation.In particular, the wave initially propagates with negligible 4 C − did not explicitly appear in our derivation above.For short waves, integration of MHD equations along C − is just another way to get the relation between n, B, and β (Equation 36), and the wave evolution problem is reduced to integration along C + .distortion: each C + characteristic satisfies dξ + /dt = 0, i.e. keeps a constant coordinate ξ = t − r/c = ξ i , and one can define an initial (undeformed) profile of E(ξ i ).
More precisely, the profile of rE(ξ) is static at small r while the normalization of E(ξ) is decreasing, E ∝ r −1 .This initial profile is the only parameter of the problem besides the magnetospheric parameters µ = r 3 B bg and D = r 3 σ bg .It is conveniently described by We will calculate the wave evolution by tracking the C + characteristics, each described by its initial position ξ i and K(ξ i ).Note that K(ξ i ) may be positive or negative.
Characteristics with K < 0 will develop κ ≪ 1, leading to monster shock formation.
It is easy to see that even an arbitrarily high σ bg → ∞ does not save the wave from breaking.Indeed, consider a wave with E oscillating between ±E 0 .As the wave expands from small radii, where E 0 (r) ≪ B bg (r), the ratio E 0 /B bg ∝ r 2 grows and eventually approaches 1/2 at some radius R × .At this moment, the minimum An arbitrarily high σ bg does not prevent the decrease of σ = κσ bg , which reduces the speed of C + characteristics, dr + /dt < c.The characteristics become bent and eventually collide, forming a shock (Figure 3).A formal proof of shock formation is provided by the solution of the coupled Equations ( 43) and ( 44) for κ(t) and ξ + (t) or r + (t) = c(t − ξ + (t)) (Appendix A).Numerical examples will be given in Section 4.

METHOD FOR TRACKING SHOCKS 3.1. Shock strength
The plasma speed β is discontinuous at the MHD shock because the upstream and downstream characteristics bring to the shock different values of β: β u ̸ = β d (subscripts "u" and "d" refer to the immediate upstream and downstream of the shock).The corresponding jump in the proper density characterizes the shock strength, It is related to the Lorentz factor of the upstream plasma relative to the downstream, Γ u|d = γ u γ d (1 − β u β d ).A similar expression for Γ u|d holds in the shock rest frame K ′ .To distinguish between different frames, dynamical quantities measured in the shock frame will be denoted with a prime, and quantities measured in the drift frame K (in which Ẽ = 0) are denoted with tilde.Using γ ′ u , γ ′ d ≫ 1, one can write the continuity of mass flux as ρu γ ′ u ≈ ρd γ ′ d , and find

Shock speed
We can track the evolution of waves with shocks if we know the shock Lorentz factor relative to the downstream plasma γ sh|d .The standard jump conditions for perpendicular relativistic magnetized shocks give γ sh|d ≈ √ σ u ≫ 1.The shock is collisionless and mediated by Larmor rotation, so the jump normally occurs on the Larmor scale, as verified by detailed kinetic simulations (Sironi et al. 2021).Dissipation of the upstream kinetic energy occurs through gyration and thermalization enabled by the instability of coherent gyration in the downstream.54).
However, the standard picture of Larmor-mediated shocks can fail for the monster shocks, where particles experience extremely fast radiative losses.The ultrarelativistic plasma flow entering the shock is subject to radiation reaction in response to curvature in the particle trajectory.The upstream flow may experience a strong radiative drag before it completes one Larmor rotation and joins the downstream.
The resulting jump conditions can be evaluated by going beyond the MHD description and examining microphysics of the shock transition on scales smaller than the particle Larmor radius.This is done in Appendix B, where we calculate the structure of the flow across the jump numerically and also derive an approximate solution analytically.The result is shown in Figure 4, and the jump condition may be stated as which smoothly matches the results at χ ≪ 1 and χ ≫ 1.The parameter χ is defined in Equation (B29) in terms of the upstream Lorentz factor γ ′ u and magnetic field where Bu = κ u B bg and σ u = κ u σ bg .Equations ( 54) and ( 55) can be solved for γ sh|d and χ, which gives This expression is found for χ ≫ 1, and the solution for χ < 1 is not needed as its value makes no difference for The shock Lorentz factor in the lab frame is where we used 1 + The shock motion in the ξ-coordinate is described by dξ sh /dt = 1 − β sh ≈ (2γ 2 sh ) −1 .Thus, we find There is one caveat in the derivation of this result: we neglected that the shock can generate an electromagnetic precursor, which interacts with the upstream flow and may reduce its Lorentz factor γ u (Section 5.5).Note however that dξ sh /dt is independent of κ u ≈ (2γ u ) −1 when χ ≫ 1.Furthermore, the shock speed will be dynamically regulated to a value independent of κ u and κ d when κ d ≪ 1, as explained in Section 5.4.

Tracking waves with shocks
The MHD wave evolution is controlled by the flow of C + characteristics in coordinate ξ = t−r/c.This flow is described by Equations ( 43) and ( 44).They determine both the shape of characteristics and the values of κ, β, and E on each C + .
After caustic formation, the born shock separates the C + flow into two regions: upstream and downstream (ahead and behind the shock).The characteristics propagate with different speeds in these regions and collide at the shock.The colliding characteristics C + u and C + d are terminated and disappear from the wave evolution problem.The location of C + u -C + d collision moves with speed β sh determined by Equation ( 58).The shock always propagates faster than C + u and slower than C + d : 43), (44), and (58) give a closed description for MHD waves with shocks.The cold approximation γ 2 s = σ used in Equations ( 44) turns out sufficient for the kHz waves studied in this paper.The accurate γ s < ∞ is important in the pre-shock region where κ drops so much that the C + flow experiences significant deformation dξ + /dt (Figure 3).This region is also coldest (due to adiabatic cooling accompanying the drop in κ), and so γ 2 s = σ.In the post-shock region ξ > ξ sh , fast radiative losses also allow one to use γ 2 s ≈ σ.Here, the evolution is anyway simple, not sensitive to the precise γ s : the post-shock γ s is so high that C + remain practically static in the ξ coordinate during the main shock dissipation phase r < ∼ 3R × .

Numerical implementation
The described method for calculating the wave evolution is easily implemented in a numerical simulation.When launching a wave, we set up an initially uniform grid in ξ i of size N + , and then use the N + characteristics to track the wave evolution in the ξ coordinate.In the simulation presented below we used N + = 105 .At each timestep dt, the displacement dξ + of each characteristic and the change of its κ are determined by Equations ( 44) and ( 43).After each timestep, the code examines the updated positions or the characteristics and checks for their crossing to detect shock formation.Once the shock is born, the code begins to track its motion according to Equation ( 58) and also check at each timestep which characteristics terminate at the shock.
The simulation thus follows the entire evolution of the wave, from its initial deformation-free propagation at r < R × to shock formation at r ≈ R × to subsequent evolution with the embedded shock.We use an adaptive timestep to resolve fast changes in MHD quantities that occur near R × .Note also that the density of characteristics drops ahead of the shock, where κ is lowest and dξ + /dt is highest.To maintain sufficient spatial resolution we use adaptive mesh refinement in ξ i without changing the total number N + of active (not terminated) characteristics.This is achieved by launching new characteristics in the region of low resolution while discarding the characteristics terminated at the shock.
The calculated ξ + (t) and κ(t) along each C + determine the wave profile κ(ξ) at any time t.Thus, we find the evolution of the wave profile.

SAMPLE NUMERICAL MODELS
As a concrete example, consider a wave with the initial profile E(ξ i ) = E 0 sin(ωξ i ).Recall that E ∝ r −1 along each characteristic C + until it develops an extremely low compression factor κ where B 2 − E 2 approaches zero.The initial profile of E sets the parameter K = rE/µ on each C + , The wave trough Radius R × is also related to the isotropic equivalent of the wave power L ≈ cr 2 E 2 0 /2 and the magnetic dipole moment of the star µ: The background magnetization parameter at R × is determined by Equation ( 8), Our sample models assume the plasma density parameter N = 10 39 and 10 37 (Section 2.1).In both models, the magnetosphere has the dipole moment µ = 10 33 G cm 3 .The wave has frequency ν = ω/2π = 10 kHz and initial power L = 10 43 erg/s.The simulation results are shown in Figure 5 and may be summarized as follows.
The wave travels with no distortion until it comes to R × where E 2 nearly reaches B 2 .The plasma at this point develops an ultra-relativistic drift with u ≈ −(1/2) B bg /(B − |E|), leading to immediate shock formation.In particular, in the model with N = 10 39 the plasma four-velocity upstream of the shock reaches u u ∼ −10 5 , and the shock compression factor q = κ d /κ u exceeds 10 4 .The shock has a large radiative parameter χ ≫ 1(Figure 6).An analytical derivation reproducing the numerical results will be given in Section 5. We will show there that the Lorentz factor of the accelerated plasma scales as γ ∝ L/N ν, and N ∼ 10 37 gives so high γ that the MHD description breaks and the plasma dynamics needs a two-fluid description.
As one can see in Figure 5, the wave profile E(ξ) develops a plateau of E ≈ −B bg /2, which corresponds to E 2 ≈ B 2 .The small variation δE along the plateau implies that the plateau is formed by stretching a small interval δξ of the initial profile E(ξ).The stretching is clearly demonstrated by the C + flow.The part of the wave profile where E 2 approaches B 2 develops a low κ 2 = (B − |E|)/B bg , and here dξ + /dt steeply increases (Equation 44).As a result, characteristics with E 2 ≈ B 2 swiftly "fall" toward the shock, with acceleration, like a waterfall.This effect protects the MHD wave from breaking the condition E 2 < B 2 .Where κ becomes so low that κ 3 σ bg ≪ 1, the C + motion in ξ saturates at dξ + /dt ≈ 2. The fast, large displacement of As the wave propagates to larger radii and B bg ∝ r −3 decreases, the plateau level E = −B bg /2 moves up and away from the original minimum of the sine wave (E = −E 0 ), and hence its width grows, so the plateau is forced to occupy an increasing part of the wave profile.At radii r ≫ R × , where B bg /2E 0 = R 2 × /r 2 ≪ 1, the plateau approaches E ≈ 0, i.e. the part of the wave with E < 0 becomes erased.As a result, the wave loses about half of its original energy. 5Waves with multiple oscillations lose even more than half: one can see that the shock formed in the first oscillation at ξ ≈ 3π/2ω Five snapshots are shown, when the wave packet reaches r/R× = 0.9 (black), 1 (red), 1.1 (green), 1.5 (blue), and 3.7 (magenta).Electric field E is normalized to E0 that would be the wave amplitude if it propagated in vacuum.The plasma Lorentz factor γ is related to the proper compression κ by γ = (1 + κ 2 )/2κ.Black dotted curves show the analytical result for κ(ξ) (Equation 78).The simulation neglected the shock precursor effect, which can reduce γ in the interval 3π/2ω < ξ < ξ sh (Section 5.5).Right: Same wave but now launched into the magnetosphere with N = 10 37 .Here, κ becomes extremely small, breaking the MHD description and transitioning to a two-fluid regime.We argue in Section 5.3 that the two-fluid calculation will likely give the same evolution of E(ξ) and κ(ξ) as found in the MHD model.eventually enters the second oscillation ξ > 2π/ω before stalling there, so the final plateau of E = 0 occupies slightly more than half of the oscillation.
The plasma speed β = E/B = E/(B bg +E) is positive where E > 0 and reaches maximum at the wave crest.This maximum β approaches unity at r ≫ R × .Here, the plasma speed relative to the wave, 1 − β ≈ B bg /E, becomes small, increasing the time it takes the plasma to cross the first half of the wave oscillation, The short-wave approximation t cross ≪ r/c holds throughout our simulation.The simulation did not follow how the plasma eventually becomes trapped in the wave, t cross > r/c; this occurs later, when the wave propagates to larger radii.Besides the monster shock launched at ξ ≈ 3π/2ω at radius R × , the simulation shows gradual steepening of the wave at the leading edge ξ = 0 (and at ξ ≈ 2π/ω).At radius R F ≈ 10R × , the wave launches a forward shock as shown in Section 5.6.

ANALYTICAL DESCRIPTION
The wave evolution demonstrated by the numerical simulation may also be derived analytically.Two key dimensionless parameters of the problem are where B × ≡ B bg (R × ) and ω × ≡ eB × /mc.

Shock formation: caustic in the C + flow
Each C + starts at small radii with ξ + = ξ i , and one can think of ξ i as a Lagrangian coordinate in the flow Figure 6.Evolution of shock parameters in the simulation with N = 10 39 : upstream Lorentz factor γu, shock compression factor q ≡ κ d /κu, and radiative parameter χ.Red dotted curves show the analytical results for γu(r), q(r), χ(r) derived in Section 5 (Equations ( 75), ( 92), ( 94)).The simulation neglected the precursor effect (Section 5.5).
of characteristics ξ i → ξ + (ξ i , t) with initial condition ξ + (ξ i , 0) = ξ i .The C + flow satisfies the relation (A11) found from Equations ( 43) and (44): where D = σ bg r 3 = σ × R 3 × .The caustic in the C + flow (birth of the shock) appears where (∂ξ + /∂ξ i ) t vanishes.It means that characteristics carrying different values of K begin to cross, creating a discontinuity.This happens at some ξ i = ξ c i and time t = t c , which may be derived analytically (Appendix A).We also find the plasma compression factor at the caustic, κ c .
In the leading order of the small parameters ζ and (R × ω/c) −1 we obtain Powerful kHz waves have √ ζ ≪ c/R × ω.In this regime, C + characteristics turn back toward the star (β + < 0) before forming the caustic.The caustic occurs at ξ c near 3π/2ω:

Plasma motion upstream of the shock
As one can see in Figure 5, the plasma develops a very low compression factor κ ahead of the shock, a result of huge plasma acceleration toward the star.The acceleration occurs along the plateau of E ≈ −B bg /2, and the four-velocity u of the plasma drift in the wave is controlled by how close B 2 − E 2 approaches zero: The plasma acceleration is accompanied by its expansion by the factor of κ −1 ≈ 2γ, as follows from the relation It is the drop of κ that boosts the motion of C + characteristics in the ξ-coordinate, deforming the wave profile.
The relation between κ and the displacement of characteristics ξ + − ξ i is given by Equation (66).Substituting K = K 0 sin(ωξ i ) and One can use this relation to evaluate κ u on the characteristic C + u that reaches the shock at time t: ξ u + = ξ sh .We are interested in waves that develop κ u ≪ 1, and then Equation (72) gives where ∆ξ ≡ ξ u + − ξ u i is the displacement of C + u .For the nascent shock at the caustic we find ω∆ξ ≈ 12 √ ζ/24 3/4 , which corresponds to κ u = κ c (Section 5.1).Later, ∆ξ(t) approximately equals the width of the plateau, along which κ drops to κ u (Figure 5).All C + populating the plateau have approximately the same E and so nearly the same ξ i -the plateau forms by the huge stretching (by the factor of ∼ κ −1 ) of an initially small interval in the original wave profile δξ i .This stretching occurs because the characteristics "fall" onto the shock from ξ i where E approaches −B bg /2, i.e. sin(ωξ i ) ≈ −B bg /2E 0 .The plateau begins and ends where sin(ωξ) ≈ −B bg /2E 0 ≈ −R 2 × /r 2 , so its width is ω∆ξ ≈ π − 2 arcsin(R 2 × /r 2 ).Then, from Equation ( 73) we find the plasma Lorentz factor just upstream of the shock, This expression accurately reproduces the simulation results (Figure 6).The maximum The same γ u can be found from energy conservation, as the flow is accelerated at the expense of the electromagnetic wave energy E = c r 2 E 2 dξ (isotropic equivalent).The profile E(ξ) includes the plateau part (E ≈ E p = −B bg /2) and the part where E(ξ) is practically unchanged from the vacuum propagation E(ξ) ≈ E 0 sin(ωξ) (Figure 5).The plateau has the width ∆ξ and ends at the shock ξ sh with a jump , where E u ≈ E p .During time dt ≈ dr/c, the "area" associated with the integral r 2 E 2 dξ changes: (1) the rise of the plateau erases a horizontal stripe |d(r 2 E 2 p )| ∆ξ and (2) the motion of the shock erases a vertical stripe r 2 [E 2 ] dξ sh .So, the wave loses energy The energy lost at ξ < ξ sh (the first term on the r.h.s.) goes into plasma acceleration upstream of the shock, where dN = 4πr 2 cn bg dt = 4πcN d ln r is the particle number (isotropic equivalent) passing through the wave during dt, and we used the fact that E 2 p r 2 ∝ r −4 .Equations ( 77) and ( 75) give the same γ u .
The consideration of the C + flow or energy conservation gives a simple solution for γ(ξ) along the plateau.In particular, using κ ≪ 1 and sin 2 (ωξ i ) ≈ R 4 × /r 4 , we find from Equation (72), The solution reproduces γ(ξ) found in the simulations (Figure 5).It ends at the shock with γ(ξ sh ) = γ u .Note that the plasma magnetization in the wave σ = κσ bg is lowest just upstream of the shock, σ u = κ u σ bg ≈ σ bg /2γ u .Its minimum value is reached soon after the wave crosses R × : σ min ∼ ζσ bg = ωR × /c ≫ 1.This minimum value is independent of σ bg .

Beyond single-fluid MHD
The accelerated flow experiences dramatic expansion in the region ξ sh < ξ < ξ i , by the factor of κ −1 ≫ 1, and the gyro-frequency in the fluid frame ωB = κω B drops.Note that the proper time t (measured in fluid frame) slows down by a similar factor, d t = dt/γ ≈ 2κ dt.
Consider now the profile of γ(ξ) ≫ 1 (Equation 78).A fluid element crosses dξ in time dt = dξ/(1 − β) ≈ dξ/2, and we find that γ changes along the fluid streamline on the proper timescale tev given by t−1 ev ≡ MHD description assumes that the gyration timescale ω−1 B ≈ (ω B /2γ) −1 is much shorter than tev .This condition may be stated as It requires γ ≪ γ MHD , where The MHD solution for the wave propagation problem remains mathematically well behaved even when γ reaches arbitrary high values.However, its applicability to a real plasma is vindicated only if γ ≪ γ MHD , which roughly corresponds to ζ ≫ (2ω/ω × ) 1/2 .Our sample numerical model with N = 10 39 had the peak γ ≈ 10 5 , just below γ MHD ≈ 2 × 10 5 , so the MHD description is marginally valid near the peak and accurate at other radii.The model with N = 10 37 developed γ ≫ γ MHD , breaking the MHD condition.
In waves with ζ < ∼ (2ω/ω × ) 1/2 , where γ reaches γ MHD , the e + and e − motions are no longer the common E ×B drift, and they are no longer coupled to the magnetic field via gyration.The momenta of unmagnetized e + and e − become significantly different and a two-fluid description is required to formulate the electric current j and its effect on the electromagnetic field of the wave.Main features of the two-fluid solution may be anticipated using the following two considerations.
(1) The symmetry of e + and e − motions in an electromagnetic wave implies that they have equal velocities β z along the wavevector k and opposite transverse velocities ±β x with β y = 0 (Beloborodov 2022).We here use local Cartesian coordinates x, y, z with the z axis along e r , y along e θ , and x along −e ϕ .The symmetry implies electric current j = j e x where j = cenβ x .Where vacuum propagation would give E 2 > B 2 , current is excited in the plasma to limit the growth of E 2 to E 2 ≈ B 2 .The two-fluid plasma with γ > γ MHD can enforce E 2 < ∼ B 2 similar to normal MHD.The ceiling of E 2 ≈ B 2 corresponds to the plateau E ≈ −B bg /2, and j needed to sustain E ≈ −B bg /2 is found from Equation ( 23): The plasma accelerated in the wave flows with speed β z ≈ −1 and density n ≈ n bg /(1 − β z ) ≈ n bg /2 (Equation 24), so the current j corresponds to Note that |β x | ≈ γ −1 MHD .If E 2 significantly exceeded B 2 , particles would be accelerated to a much greater β x .This cannot happen, because the small β x ≈ −γ −1 MHD already gives j sufficient to enforce the ceiling E 2 ≈ B 2 .Thus, the wave profile E(ξ) in the two-fluid regime should evolve similarly to MHD -an excess of E 2 > B 2 will be shaved off at E ≈ −B bg /2.
(2) Reducing E 2 to satisfy the ceiling of E 2 ≈ B 2 implies a well-defined loss of electromagnetic energy of the wave and the corresponding energy gain by the plasma (Section 5.2).Energy conservation then determines the profile of γ ahead of the shock just like in the MHD regime, which gives Equation (78) for γ(ξ).
So, we expect that the analytical expressions for γ, β x , and u x = γβ x derived in MHD should carry over to the unmagnetized (two-fluid) regime.The transition from MHD to the two-fluid regime occurs at |u x | ∼ 1.The MHD condition |u x | ≪ 1 may also be stated as which is equivalent to γ ≪ γ MHD (Equation 81).Note that . the fluid energy arises from its bulk motion along z while the internal transverse motions of e ± are negligible.By contrast, in the unmagnetized regime the e ± develop |u x | ≫ 1, i.e. the internal motions become relativistic and significantly contribute to the plasma energy.This change will affect the shock jump conditions; however the shock speed should anyway be regulated to the value given in Equation ( 90) below.
It is also instructive to look at the dynamical equations describing the two-fluid e ± motions.Taking into account the symmetry of e ± velocities, it is sufficient to consider the positrons.The energy equation in the plateau region E ≈ −B bg /2 reads where ω B = eB bg /mc and d/dt is taken along the particle trajectory, so dt = dξ/(1 − β z ) ≈ dξ/2.Substituting β x from Equation ( 83), one finds recovering the solution (78) for γ(ξ).Sustaining β x according to Equation ( 83) requires a small mismatch between β z and β D = E/B, as seen from the x-momentum equation, Let us define ∆ ≡ −(B + E)/B as a measure for the deviation of E from −B; ∆ > 0 corresponds to E 2 > B 2 .It is easy to show that ∆ > 1 + β z ≈ (2γ 2 z ) −1 is required to achieve u x of the unmagnetized regime, so du x /dt ≈ −(ω B /2)∆.The electric current j and the corresponding β x (Equation 83) will be sustained if The tiny positive ∆ shows that E 2 slightly exceeds B 2 in the unmagnetized regime.By contrast, in the MHD regime, ∆ ≈ −(2γ 2 ) −1 < 0, i.e.E 2 < B 2 .In both cases, the ceiling of E 2 ≈ B 2 is strongly enforced.

Shock strength
The upstream flow is decelerated in the shock: we observed in the simulation a huge jump of the plasma proper density q = ρd /ρ u , which implies a huge jump of Lorentz factor γ u /γ d ≫ 1.By contrast, the electric field has a small jump during the peak of the monster shock: . This means that the plasma immediately behind the shock still moves with a large speed: The ratio γ u /γ d is regulated by the shock jump conditions.They require fast motion of the shock relative to the downstream plasma: The value of κ d ≪ 1 may be derived as follows.The condition that This relation controls the shock position ξ sh as a function of time t ≈ r/c and thus determines its speed dξ sh /dt = 1−β sh ≈ (2γ 2 sh ) −1 (when the shock Lorentz factor γ sh ≫ 1).This gives On the other hand, the shock jump conditions require dξ sh /dt given in Equation ( 58).Substituting it into Equation (90), we find where η = (σ T B × /8πe).
The obtained κ u and κ d (Equations ( 75) and ( 91)) This result is in excellent agreement with the numerical simulation (Figure 6).It loses accuracy at large radii where E d significantly deviates from −B bg /2, so that the approximation (89) becomes invalid.This occurs when κ d (Equation 91) increases to ∼ 1, i.e. at radius Note also that κ d ≈ 1 corresponds to q ≈ 2γ u .For the wave simulated in Section 4 with N = 10 39 we find x 1 ≈ 1.85 which corresponds to r 1 = x 1 R × ≈ 2.6 × 10 8 cm.The shock radiative parameter χ is given by Equation (56), as long as χ ≫ 1. Substitution of Equations ( 75) and ( 91) for κ u and κ d yields It reproduces the value of χ observed in the simulation.

Precursor effect
In an accompanying paper we describe an additional effect: the shock emits precursor waves into the upstream, which can significantly decelerate the upstream flow before it reaches the collisionless shock.
The precursor emission begins immediately after shock formation at time t × ≈ R × /c near the coordinate ξ 0 = 3π/2ω.As the wave expands to r > R × , the shock moves to ξ sh > ξ 0 , so the precursor occupies the growing region Note that the plasma flow in the MHD wave crosses half of the plateau E ≈ −B bg /2 ahead of the precursor, at ξ < ξ 0 , and begins to interact with the precursor at ξ 0 with the Lorentz factor determined by Equation ( 78), It is lower by a factor of 2 compared to γ u defined previously at ξ sh (Equation 75).In the accompanying paper we find that the precursor can decelerate the upstream flow as it moves from ξ 0 to ξ sh , so the flow comes to the shock with γ(ξ sh ) ≪ γ(ξ 0 ) (and the shock radiative parameter χ may drop below unity).This pre-deceleration effect depends on somewhat uncertain efficiency of precursor emission, which has not been studied in the extreme regime of the monster shocks.Regardless of the shock structure details, with or without the precursor, the basic picture remains the same: the upstream first develops the huge Lorentz factor γ(ξ 0 ) [or 2γ(ξ 0 ) without the precursor] and then promptly radiates its energy, before completing the shock transition.The emitted radiation is in the X/gamma-ray band, as shown in Section 6 below.

Forward shock
The plateau evolution described above leads to erasing half of the wave oscillation where E < 0. It does not affect the leading half 0 < ωξ < π where E > 0 (Figure 5).This leading part of the wave propagates with negligible distortions at radii r ∼ R × , nearly as in vacuum.However, at larger r ≫ R × the profile of E(ξ) > 0 gradually steepens at the leading edge ξ = 0 and eventually forms a forward shock.Below we find radius R F where this forward shock is launched.
In the zero order of the small parameter dξ + /dt, C + are described by ξ + = ξ i or r + = c(t−ξ i ).The correction ξ + − ξ i = (κ − 1) 2 /8cDK 2 κ (Equation 66) in the leading order may be found iteratively by substituting κ 2 = 1 + 2Kr 2 with the zero-order r = c(t − ξ i ).This gives One can now examine the flow of C + characteristics ξ i → ξ + (ξ i , t) and identify where (∂ξ + /∂ξ i ) t vanishes, launching a shock.For the C + flow with σ bg κ 3 ≫ 1 one can use Equation ( 97) and find We here focus on the interval 0 < ωξ i < π and observe that ∂ξ + /∂ξ i drops fastest at the leading edge ξ i = 0, where κ = 1 is minimum and dK/dξ i = K 0 ω = ω/2R 2 × is maximum.We find at the leading edge (using Hence, the forward shock forms at radius The wave now carries two shocks, at ξ ≈ 0 and ξ > 2π/ω, which slowly shift to larger ξ.The nascent forward shock at r = R F is weak.It becomes ultrarelativistic at r ≫ R F and later turns into a strong blast wave expanding into the wind outside the magnetosphere. 6. GAMMA-RAYS FROM MONSTER SHOCKS First, consider shocks with neglected effects of the precursor on the upstream flow.As the plasma flow crosses the shock, it radiates its energy in curvature photons.The spectrum of curvature radiation from a particle with Lorentz factor γ cuts off exponentially at ω > ω c (Landau & Lifshitz 1975), where The radiated power Ėe = γmc 2 is Lorentz invariant, and may be evaluated in the shock frame, γ = |dγ ′ /dt ′ |, using the solution derived in Appendix B. It shows that γ ′ drops from the upstream value γ ′ u ≈ 2γ u γ sh ≈ γ sh /κ u with rate γ ≈ 4 3 where g = γ ′ /γ ′ u , ω B = eB bg /mc, and Bu = κ u B bg .The emission frequency ω c is maximum at g = (5/8) 2/3 , Substituting the solutions for γ u (r) (Equation 75), γ sh (r) (Equation 90), and σ u = κ u σ bg ≈ σ bg /2γ u , we find where One can see that ω max c is in the far gamma-ray band.Next, consider emission from shocks affected by the precursor (Section 5.5).Main emission now occurs where the upstream flow enters the precursor and begins deceleration with γ ≈ γ(ξ 0 ) given by Equation (96).The emission frequency is still given by Equation ( 101), but now with Ėe being the emitted power that results from the particle interaction with the precursor.Using dξ = (1 − β)dt ≈ 2dt along the plasma streamline, one can see that δξ dec ≡ 2γmc 2 / Ėe defines a characteristic deceleration scale, and we rewrite Equation (101) as During the main phase of shock evolution, the upstream flow reaches the peak Lorentz factor γ(ξ 0 ) ≈ (2ζ) −1 and radiates photons with dimensionless energy where we have normalized δξ dec to the precursor width ξ sh − ξ 0 ∼ ω −1 .Note that the precursor deceleration effect is strong if δξ dec ≪ ω −1 , and so ϵ c ≫ 10 −8 ζ −3/2 .This gives the characteristic ϵ c in the gamma-ray band.
7. DISCUSSION 7.1.Formation of shocks The magnetospheres of neutron stars have a huge magnetization parameter σ bg , and therefore their lowfrequency perturbations are often described as FFE waves, which propagate with the speed of light.This description is excellent near the star, however it fails when perturbations propagate to larger radii and grow in relative amplitude E/B bg .Then, FFE becomes remarkably self-destructive: it pushes itself out of the realm of its applicability E 2 < B 2 , and the wave dynamics in the FFE limit σ bg → ∞ becomes undefined.Therefore, waves should be described using the full MHD framework, with an arbitrarily large but finite σ bg .
In particular, compressive MHD waves (vacuum radio waves in the FFE limit) steepen into shocks at radius R × (Equation 60) for any high σ bg .The higher σ bg the stronger the shock at R × , as we have demonstrated by solving the MHD wave equation.As an example, we presented the evolution of a 10-kHz wave with an initial sine profile E = E 0 sin(ωξ), where ξ = t − r/c.We have also derived the wave evolution analytically to show how it depends on the wave frequency ω and power L. Our main conclusions are as follows.6 1.When the wave reaches R × ∼ 10 8 cm it suddenly begins to pull the background magnetosphere toward the star, creating an ultrarelativistic flow at the wave oscillation phase near 3π/2 (trough), where E 2 approaches B 2 (Figure 5).The flow acceleration at the trough of the wave is accompanied by plasma expansion, reducing the local magnetization parameter σ, so the wave propagation slows down at the trough and the wave "stumbles," steepening into a shock.The accelerated plasma flow forms the upstream of the shock.It develops a huge Lorentz factor γ u ∼ cσ × /ωR × , where σ × ≡ σ bg (R × ).The peak γ u is reached when the wave crosses r ≈ 1.15R × , and here the local magnetization parameter of the flow is reduced to σ u ∼ σ × /γ u ∼ ωR × /c, which independent of σ bg .
2. The monster shock has an unusual micro-structure.The upstream flow experiences strong radiative losses before completing a single Larmor rotation, i.e. before forming the downstream MHD flow.This affects the shock jump conditions.In particular, we found that the shock Lorentz factor relative to the downstream plasma is γ sh|d ≈ (1 + χ) 1/7 √ σ u , where χ is a new dimensionless radiative parameter (Equation 55).The shock structure is further complicated by the precursor emission, which can induce radiative losses of the upstream flow ahead of the shock; this effect will be described elsewhere.
3. As the kHz wave propagates beyond R × , half of its oscillation (π < ωξ < 2π) becomes erased (Figure 5).Thus the wave loses half of its energy after crossing R × .The large wave energy per magnetospheric particle corresponds to the huge Lorentz factor gained by the plasma.The plasma radiates the gained energy in the ultra-relativistic shock, leading to a bright X-ray burst.
4. Waves accelerating the plasma to γ ∼ (2ω B /ω) 1/2 enter the two-fluid regime with unmagnetized particles.We argued in Section 5.3 that the main MHD results will carry over to the two-fluid regime.In particular, the Lorentz factor of upstream particles will be given by the same expression γ ∼ cσ × /ωR × , reaching enormous values in powerful low-frequency waves.
5. The leading half of the wave oscillation (ωξ < π) never approaches the condition E 2 = B 2 .It crosses R × without any significant changes.A wave consisting of only this half-oscillation would not develop the monster radiative shock at R × .However, in any case, at a larger radius R F (a few times 10 9 cm, see Equation ( 100)), the leading edge of the wave steepens into a forward shock.
Both R × and R F are well inside the typical light cylinder of a magnetar, R LC ∼ 10 10 cm.Thus, strong kHz waves from magnetars inevitably launch shocks inside the magnetosphere.The monster shock forms at ωξ ≈ 3π/2; then it weakens and shifts to ωξ ≈ 2π.The forward shock forms at ξ ≈ 0 and gets stronger at large radii; it will evolve into an ultrarelativistic blast wave expanding into the magnetar wind (this evolution will be further described elsewhere).
Our results clarify the relation between two pictures proposed for electromagnetic ejecta from magnetars: a vacuum-like electromagnetic wave (Lyubarsky 2014) vs. a blast wave in the wind (Beloborodov 2017(Beloborodov , 2020)).Blast wave formation is important for the theory of fast radio bursts (FRBs, see Lyubarsky 2021 for a review), because it produces semi-coherent radio emission with sub-ms duration.Kinetic simulations demonstrating the efficiency and polarization of this emission are found in Sironi et al. (2021).In addition, Thompson (2023) recently proposed that the blast wave may produce radio waves by a different mechanism if it expands into a turbulent medium.

Compton drag
Our calculations of magnetospheric MHD waves neglected Compton drag exerted by ambient radiation.Radiation flows from magnetars with luminosities L ⋆ ∼ 10 35 erg/s, peaking at photon energies ε ⋆ mc 2 ∼ 1 keV (Kaspi & Beloborodov 2017).When plasma in the wave develops Lorentz factor γ ≫ 1, it upscatters the keV photons to energies ε sc ∼ γ 2 ε ⋆ as long as γ < ε −1 ⋆ .This gives a maximum power of scattered radiation ⋆ L ⋆ ∼ 3 × 10 40 erg/s (this upper bound assumes that each keV photon is upscattered).Waves with power L ≫ L sc are weakly affected by Compton drag.When γ > ε −1 ⋆ scattering occurs with a smaller cross section σ sc ∼ σ T /γε ⋆ and ε sc ∼ γ.Using the scattering optical depth τ sc ∼ σ sc N /R 2 × , we find This again gives L sc ≪ L, even if the density parameter N = nr 3 is increased by secondary e ± creation (which mainly occurs behind the wave, not inside it, see Section 7.4 below).Thus, we conclude that the magnetar radiation does not prevent the enormous plasma acceleration in magnetosonic waves at r ≈ R × .
7.3.Wave evolution at θ ̸ = π/2 In this paper, we calculated the wave propagation in the equatorial plane of a dipole magnetosphere, θ = π/2.The obtained solution shows a clear physical picture of the wave evolution.It includes two facts that help extend the picture to θ ̸ = π/2: (1) The MHD wave is well described as a vacuum electromagnetic wave superimposed on the background dipole magnetosphere until this superposition hits the condition E 2 = B 2 , launching the monster shock.In the equatorial plane, this condition is reached at R × .In other parts of the magnetosphere, one can find the location of shock formation from the same condition E 2 ≈ B 2 , which corresponds to drift speed |β D | → 1.The drift velocity is given by with E = −Ee ϕ and B = B bg + Ee θ .The wave propagates as in vacuum until β 2 D = E 2 /B 2 approaches unity, which corresponds to B 2 bg +2EB bg,θ = 0.This condition is first approached at the trough of the wave, E = −E 0 .It defines the critical surface r × (θ): A spherical wave with approximately isotropic power reaches this surface first at θ = π/2, r = R × .As the wave continues its propagation to r > R × , the monster shock develops in a range of latitudes, sin θ × < sin θ ≤ 1 (Figure 7), where sin θ × (r) is found by solving Equation (110) for sin θ.
The condition E 2 → B 2 corresponds to the plasma velocity β approaching the unit vector,7 This unit vector represents the direction of the accelerated plasma flow.Note that the radial component of β D changes sign at sin θ = 2/ √ 5. (2) The obtained solution at θ = π/2 demonstrates that the wave evolution outside R × follows a simple principle: the part ∆ξ of the wave oscillation that has hit the ceiling E 2 = B 2 forms the plateau of E ≈ E p while the rest of the oscillation profile E(ξ) follows nearly vacuum propagation (until the wave approaches R F ≫ R × ).A similar description can be used at θ ̸ = π/2.The plateau electric field is set by the condition E 2 ≈ B 2 : It is sustained by electric current j p = −∂ r (rE p )/2πr.The plasma passing through the wave is accelerated on the plateau with rate dγ/dt = E p • j p /ρc 2 , which gives where we used dξ = (1 − β r )dt = (ρ bg /ρ)dt along the plasma streamline.A plateau of width W p = c∆ξ gives The acceleration region W p with E ≈ E p and current j p appears at r = r × and grows at r > r × .For example, waves with initial profile E(ξ) = E 0 sin(ωξ) develop plateaus with W p (r) ≈ (c/ω)[π − 2 arcsin(r 2 × /r 2 )].As a result of the wave evolution at r > r × , the parts of the wave profile with E •(k ×B bg ) < 0 become erased and replaced by the plateau with E p /E 0 → 0 at r ≫ r × .In axisymmetric waves E oscillates along e ϕ , and so the parts with E ϕ > 0 become erased.In particular, oscillating kHz waves with average E ϕ ≈ 0 lose half of their energy.The lost energy is radiated by the shock.

Pair creation and emission of X-ray bursts
The monster shock radiates high-energy photons (Section 6), which can convert to e ± pairs.The upstream flow entering the equatorial shock moves along −e r with a high Lorentz factor γ, so its gamma-ray emission is beamed toward the star.These gamma-rays propagate perpendicular to B bg and will convert to e ± pairs.The mean free path for conversion is given by (Erber 1966), For our sample model in Section 4, l ph (R × ) ≫ R × and l ph (r) will be reduced below r when the gamma-rays propagate to r ≪ R × where B bg ≫ B × .Magnetosonic waves accelerate the plasma radially only in the equatorial plane of the magnetosphere.Waves at different latitudes accelerate the plasma obliquely (Section 7.3), so the gamma-ray emission is oblique, not beamed toward the star, and can avoid absorption by the ultrastrong field.Instead, the oblique gamma-rays will convert to e ± pairs via photon-photon collisions, as they propagate toward the equatorial plane and collide with the symmetric gamma-rays from the lower hemisphere.The collisions will be efficient because of the broad spectrum of curvature radiation (its half-width extends from 0.01ϵ c to 1.5ϵ c , see e.g.Longair (1994)), so the gamma-rays will find counterparts near the threshold for e ± creation with a large cross section ∼ 0.1σ T .This leads to a high pair production rate, and the region behind the shock will become populated with optically thick e ± plasma.
The created e ± plasma will cool by emitting softer synchrotron photons and by scattering the photons as they diffuse out of the e ± cloud behind the shock.Radiation production is controlled by the compactness parameter, The high compactness implies a large optical depth to photon-photon collisions, so most gamma-rays should convert to pairs, which experience fast radiative cooling.Thus, the shock emission will be reprocessed to photons of lower energies.The luminosity of the resulting X-ray burst is comparable to the power of the original kHz wave.The burst resembles radiative processes in compact magnetic flares simulated in Beloborodov (2021b); its spectrum can be found with similar detailed radiative transfer simulations, which we leave for future work.
A minimum duration of the burst is set by the light crossing time r/c, where r is a few R × .A typical value of this minimum duration is ∼ 10 ms.Bursts with high ℓ × may last longer, as it takes time for the reprocessed X-rays to diffuse out of the optically thick e ± cloud.In addition, starquakes may produce multiple shocks, creating multiple bursts.Such composite bursts can last ∼ 100 ms, depending on the crustal quake coupling to the neutron star core (Bransgrove et al. 2020).

X-ray precursor of a binary merger
Consider a neutron-star binary with masses M 1 and M 2 separated by distance 2a.The reduced mass of the binary is M = M 1 M 2 /(M 1 +M 2 ), and its orbital angular velocity is In a tight binary (nearing merger) Ω is in the kHz band.Suppose both stars are strongly magnetized.An interface between their magnetospheres forms where their pressures balance.For instance, consider two stars with equal magnetic dipole moments |µ 1 | = |µ 2 | = µ.Then, the interface is in the middle, at the distance a from each star, and the magnetic field at the interface is B i ≈ µ/a 3 .
If the orbiting stars are not in synchronous rotation, there is differential rotation between the two magnetospheres with an angular frequency Ω diff comparable to the orbital Ω.The two magnetic moments µ 1 and µ 2 in general are not aligned, and then the magnetic pressure at the interface will oscillate with frequency Ω diff .The pressure variation timescale ∼ Ω −1 diff exceeds the Alfvén crossing time (a few a/c), which reduces the efficiency of low-frequency wave excitation.However, the interface is also a source of higher-frequency waves, because it is prone to instabilities of the Kelvin-Helmholtz/diocotron type.The instability will generate vortices of sizes up to a fraction of a (with frequencies ω > c/a), creating traction between the rotating magnetospheres.
The turbulent region around the interface will emit both Alfvén and magnetosonic waves, and further investigation of this process requires numerical simulations.An upper bound on the power deposited into the turbulence may be roughly estimated as L max ∼ Ω diff a 3 B 2 i /8π.The power of magnetosonic wave emission L ∼ ca 2 (δB) 2 /8π < L max implies an upper limit on the emission amplitude, In a more general case the magnetic dipole moments of the two stars are not equal, µ 1 < µ 2 .Then the interface will form at distance R from M 1 , defined by (119) The wave source will have the size of ∼ R < a, and its power may be estimated as The emitted magnetosonic waves will steepen into shocks at radius R × ∼ (2b) −1/2 R (if it is inside the light cylinder of the binary, R × < c/Ω).
The emission of kHz waves accompanied by shocks and turbulent dissipation will continue as long as the two stars orbit each other, until they eventually merge.The binary loses orbital energy to gravitational wave emission and shrinks on the timescale which is also comparable to the time remaining until the merger.The wave energy generated during the remaining time t may be estimated as tL, and we obtain A significant fraction of the kHz wave energy will be dissipated and converted to X-rays.This process can produce a detectable X-ray precursor of the merger if the stars have strong magnetic fields.For instance, a binary with µ 1 = µ 2 = 10 32 G cm 3 may generate Xray luminosity L X ∼ 10 46 erg/s at t ∼ 10 s before the merger.Wave dissipation provides an alternative to the recently proposed precursor emission by magnetic flares (Most & Philippov 2020;Beloborodov 2021b).The latter mechanism was shown to operate in binaries with anti-aligned µ 1 and µ 2 ; it is driven by the over-twisting of magnetic flux tubes connecting the two stars.

Shocks from neutron star collapse
A massive neutron star can be born (e.g. in a merger) with fast rotation, which temporarily supports it against collapse.Such objects are likely strongly magnetized and gradually lose rotation by emitting angular momentum in a magnetized wind.The spindown can eventually lead to the collapse of the massive neutron star into a black hole, producing an electromagnetic transient (Lehner et al. 2012;Falcke & Rezzolla 2014;Most et al. 2018).Numerical simulations of the collapse show that it launches a strong outgoing electromagnetic pulse propagating through the outer dipole magnetosphere.A monster shock will form where the pulse amplitude becomes comparable to the dipole field of the pre-collapse magnetosphere and B 2 − E 2 reaches zero.This condition was indeed seen in vacuum and FFE simulations of the magnetosphere of a collapsing star (see Figure 10 in Lehner et al. (2012)).Demonstrating shock formation and tracking its propagation requires a full MHD calculation, as shown in the present paper.The shock will dissipate a significant fraction of the outgoing electromagnetic pulse and produce a bright X-ray transient.
I am grateful to Alex Chen for sharing their particlein-cell simulations of strong magnetosonic waves prior to publication.This work was prompted in a large part by their results, which demonstrate shock formation through first-principles, fully kinetic calculations.I am also grateful to the anonymous referee for many useful comments, to Yuri Levin and Ashley Bransgrove for discussions of fast magnetosonic waves from quakes, and to Jens Mahlmann, Emanule Sobacchi, and Sasha Philippov for reading the manuscript and providing comments.This work is supported by NSF AST 2009453, NASA 21-ATP21-0056 and Simons Foundation #446228.Each solution ends where the downstream begins, i.e.where the e ± streams develop gyration (then they become unstable, thermalizing their gyration energy in the drift frame).The drift Lorentz factor γD at the end point represents the downstream speed relative to the shock.Right: Relation between g and b.The approximate analytical solution (Equation B39) is shown by the dotted curve and compared with accurate numerical integration (solid curve).Both numerical and analytical solutions are plotted until the e ± streams reach | βx| = 1/2, which corresponds to 60 • gyration in the drift frame.
We set the end point s end of the coherent stream solution where the e ± velocity vector in the drift frame, β± , has rotated by 60 • ( β± x = ∓1/2).We define the approximate downstream quantities at this point.In particular, the drift Lorentz factor γ D = B/ √ B 2 − E 2 at s end approximately represents the Lorentz factor of the downstream relative to the shock, γ d .Figure 4 shows the measured γ d as a function χ.8For shocks with χ ≫ 1, the numerical result approximately follows the relation γ d / √ σ u ≈ χ 1/7 .It can also be derived analytically, as shown below.

B.4. Approximate analytical solution
The radiation-reaction transition admits a simple analytical description.A key fact simplifying the analytical approximation is that the term wχR in dw/ds (Equation B28) remains small throughout the shock transition.When this term is neglected, the equations for dw/ds and db/ds can be integrated for w(b).The resulting relation w(b) is the same as in the adiabatic solution (B32).At b ≫ ϵ it gives Note also that the functions R and F entering the dynamical equations can be simplified.Neglecting the small terms ∝ ϵ in Equation (B31), we obtain where we used w/σ u ≪ g to further simplify F. Like Equation (B35), Equations (B36) approximately hold at b ≫ ϵ, practically throughout the entire shock transition.

Figure 1 .
Figure 1.Excitation of magnetospheric waves by a star quake with ω ≫ c/R⋆ and a vertical (radial) wavevector k.The wave polarization is set by the direction of the electric field E = B × v/c in the crustal shear oscillation with a horizontal velocity v. Left: crustal motion with v ∥ n ≡ k × B bg excites an Alfvén wave E ⊥ n.Right: v ⊥ n excites a magnetosonic wave E ∥ n.The wave amplitude is small, so B ≈ B bg near the star.

Figure 2 .
Figure 2. A magnetosonic wave with wavevector k ⊥ B bg at two oscillation phases: when the wave field Bw = B −B bg is aligned with B bg (left) and anti-aligned with B bg (right).The wave electric field Ew = E is always orthogonal to B bg .The quantity B 2 − E 2 reaches zero when Bw = −B bg /2.At this moment, the plasma drift speed v = c E × B/B 2 approaches c, so it experiences ultra-relativistic acceleration in the direction opposite to the wave propagation direction.

Figure 3 .
Figure 3. Schematic illustration of the flow of characteristics C + in FFE and MHD, shown in the t-ξ plane.Each C +has its initial position ξi and carries a value of K(ξi) = rE/µ determined by the initial wave profile.In the FFE limit, the characteristics are vertical straight lines; they propagate with speed β+ = 1, which corresponds to dξ+/dt = 0.The MHD correction to FFE implies dξ+/dt = 1 − β+ ̸ = 0, so the C + characteristics are no longer static in ξ; they become significantly bent in the region where E 2 approaches B 2 , leading to the formation of a shock (red curve).

Figure 4 .
Figure 4. Lorentz factor of the shock relative to the downstream plasma as a function of the radiative parameter χ (solid curve, calculated as explained in Appendix B).Dashed line shows the analytical result at χ ≫ 1 (Equation54).

Figure 5 .
Figure5.Left: Evolution of the wave profile E(ξ) and κ(ξ), where ξ ≡ t − r/c.The wave has frequency ν = ω/2π = 10 kHz and initial power L = 10 43 erg/s; the magnetosphere has magnetic dipole moment µ = 10 33 G cm 3 and density parameter N = 10 39 .Five snapshots are shown, when the wave packet reaches r/R× = 0.9 (black), 1 (red), 1.1 (green), 1.5 (blue), and 3.7 (magenta).Electric field E is normalized to E0 that would be the wave amplitude if it propagated in vacuum.The plasma Lorentz factor γ is related to the proper compression κ by γ = (1 + κ 2 )/2κ.Black dotted curves show the analytical result for κ(ξ) (Equation78).The simulation neglected the shock precursor effect, which can reduce γ in the interval 3π/2ω < ξ < ξ sh (Section 5.5).Right: Same wave but now launched into the magnetosphere with N = 10 37 .Here, κ becomes extremely small, breaking the MHD description and transitioning to a two-fluid regime.We argue in Section 5.3 that the two-fluid calculation will likely give the same evolution of E(ξ) and κ(ξ) as found in the MHD model.

Figure 7 .
Figure 7. Expansion of a spherical wave front through the magnetosphere.The front thickness (wave duration) is much smaller than r/c.Its shock-free part (sin θ < sin θ×) is plotted in green and the part carrying the shock is in red.The region swept by the shock is shaded in pink; its boundary θ(r) is defined by Equation (110).The wave front is shown at two moments: (1) when it crosses radius R× (the shock has formed in the equatorial plane, θ = π/2); (2) when it expands to a larger radius (the shock now occupies a range of θ).The shock radiates gamma-rays (red arrows) beamed along β D (Equation 111); most of them convert to e ± pairs.