Introduction
Pulsar radio emission originates from the star’s magnetosphere (see [38 ] , for a review). A typical radio pulsar rotates with an angular frequency Ω ⋆ ≡ 2 π / P \Omega_{\star} \equiv2\pi/P Ω ⋆ ≡ 2 π / P P ∼ 1 s P \sim1\ \mathrm{s} P ∼ 1 s R L C ≡ c / Ω ⋆ R_{\mathrm{LC}} \equiv c/\Omega_{\star} R LC ≡ c / Ω ⋆ [20 , 46 ] ). In the closed zone, the magnetic field B \mathbf{B} B B \mathbf{B} B [12 , 21 , 45 , 47 , 26 ] ). We focus on the open field-line region.
The rotating magnetosphere requires a charge density of η c o = − Ω ⋅ B / ( 2 π c ) \eta_{\mathrm{co}} = -\boldsymbol{\Omega} \cdot\mathbf{B}/(2\pi c) η co = − Ω ⋅ B / ( 2 π c ) B \mathbf{B} B [20 ] ). The flow of this plasma can sustain a maximum (field-aligned) co-rotation current density j c o = c η c o j_{\mathrm{co}} = c\eta_{\mathrm{co}} j co = c η co j m = ( c / 4 π ) ∇ × B j_{\mathrm{m}} = (c/4\pi)\nabla\times\mathbf{B} j m = ( c /4 π ) ∇ × B j m j_{\mathrm{m}} j m α ≡ j m / j c o > 1 \alpha\equiv j_{\mathrm{m}}/j_{\mathrm{co}} > 1 α ≡ j m / j co > 1 α < 0 \alpha< 0 α < 0 [43 , 7 ] ).
Here, we focus on the polar gap above the magnetic poles of the star ([40 ] ). Other gaps include, e.g., the slot gap ([2 ] ), the outer gap ([10 ] ), and return-current gap(s) ([9 ] ).
Plasma discharges intermittently screen the gaps ([48 , 49 ] ). Quantum electrodynamic (QED) processes enable electron-positron pairs to interact with high-energy photons, populating the gap with copious plasma clouds of high pair multiplicity n ± / ∣ η c o ∣ ≫ 1 n_{\pm}/|\eta_{\mathrm{co}}| \gg1 n ± /∣ η co ∣ ≫ 1 [16 , 1 , 23 , 50 ] ), where n ± n_{\pm} n ± [14 , 51 , 37 ] ) and simulated QED cascades with radiative particle-in-cell (PIC) ([48 , 49 , 39 , 13 , 11 , 7 ] ), and Velasov methods ([53 ] ).
Here we present a new analytical concurrency model that simultaneously describes the coupled dynamics of plasma multiplicity, photon production, and electric field evolution. The discharge exhibits a stable limit-cycle behavior driven by the interplay between particle acceleration and screening. Comparison with one-dimensional radiative PIC simulations including exact QED processes validates the theoretical model and characterizes the resulting particle spectra and electromagnetic radiation.
Pulsar Magnetospheres
System Parameters
Radio pulsars are neutron stars with a radius R ⋆ ≈ 10 6 c m R_{\star} \approx10^{6}\ \mathrm{cm} R ⋆ ≈ 1 0 6 cm M ⋆ ≈ 1.4 M ⊙ M_{\star} \approx1.4M_{\odot} M ⋆ ≈ 1.4 M ⊙ [33 ] ). They rotate with a period P ∼ 1 s P \sim1\ \mathrm{s} P ∼ 1 s Ω ⋆ ≈ 6 P 0 − 1 s − 1 \Omega_{\star} \approx6P_{0}^{-1}\ \mathrm{s}^{-1} Ω ⋆ ≈ 6 P 0 − 1 s − 1 Q x ≡ Q / 10 x Q_{x} \equiv Q/10^{x} Q x ≡ Q /1 0 x Q Q Q 10 x 10^{x} 1 0 x
A typical radio pulsar has a magnetic field strength of B ⋆ ∼ 10 12 G B_{\star} \sim10^{12}\ \mathrm{G} B ⋆ ∼ 1 0 12 G B Q ≡ m e 2 c 3 / ( ℏ e ) ≈ 4.4 × 10 13 G B_{Q} \equiv m_{\mathrm{e}}^{2}c^{3}/(\hbar e) \approx4.4 \times10^{13}\ \mathrm{G} B Q ≡ m e 2 c 3 / ( ℏ e ) ≈ 4.4 × 1 0 13 G
b = B ⋆ B Q ∼ 0.02 B 12 , (1) b = \frac{B_{\star}}{B_{Q}} \sim0.02 B_{12}, \tag*{(1)} b = B Q B ⋆ ∼ 0.02 B 12 , (1) rendering QED mildly nonlinear. Far from the star, the magnetic field structure resembles a dipole, whereas close to the surface it can have higher multipolar components with deformations and magnetic loops similar to those of the Sun.
The rotating magnetosphere splits into open and closed zones (see Figure 1 ). The light cylinder radius, where azimuthal velocity v ϕ = c v_{\phi}=c v ϕ = c R L C / R ⋆ ∼ 5000 P 0 R_{\mathrm{LC}}/R_{\star} \sim5000P_{0} R LC / R ⋆ ∼ 5000 P 0 μ \boldsymbol{\mu} μ r = R ⋆ r=R_{\star} r = R ⋆ sin θ p c ≈ R ⋆ / R L C \sin\theta_{\mathrm{pc}} \approx\sqrt{R_{\star}/R_{\mathrm{LC}}} sin θ pc ≈ R ⋆ / R LC
Figure 1. Schematic visualization of the pulsar magnetosphere (thin black curves), the studied electric circuit (thick orange curve), and the relevant circuit components (see text). We simulate the one-dimensional circuit from point a to b.
R p c ≈ R ⋆ sin θ p c ≈ R ⋆ ( R ⋆ R L C ) 1 / 2 ∼ 0.01 P 0 − 1 / 2 R ⋆ . (2) R_{\mathrm{pc}} \approx R_{\star}\sin\theta_{\mathrm{pc}} \approx R_{\star}\left(\frac{R_{\star}}{R_{\mathrm{LC}}}\right)^{1/2} \sim0.01P_{0}^{-1/2}R_{\star}. \tag*{(2)} R pc ≈ R ⋆ sin θ pc ≈ R ⋆ ( R LC R ⋆ ) 1/2 ∼ 0.01 P 0 − 1/2 R ⋆ . (2) A curvature radius for a dipole field is R c u r v ≈ 4 r / ( 3 sin θ ) R_{\mathrm{curv}} \approx4r/(3\sin\theta) R curv ≈ 4 r / ( 3 sin θ ) r r r θ \theta θ R c u r v ≈ 4 3 R L C R ⋆ ∼ 9 × 10 7 P 0 1 / 2 c m R_{\mathrm{curv}} \approx\frac{4}{3}\sqrt{R_{\mathrm{LC}}R_{\star}} \sim9 \times10^{7}P_{0}^{1/2}\ \mathrm{cm} R curv ≈ 3 4 R LC R ⋆ ∼ 9 × 1 0 7 P 0 1/2 cm R c u r v ∼ R ⋆ ∼ 10 6 c m R_{\mathrm{curv}} \sim R_{\star} \sim10^{6}\ \mathrm{cm} R curv ∼ R ⋆ ∼ 1 0 6 cm
Assuming that the magnetic axis and spin axis are aligned, the polar cap’s rim rotates with a velocity
β r o t ≈ Ω ⋆ R p c c ≈ 3 × 10 − 6 P 0 − 3 / 2 , (3) \beta_{\mathrm{rot}} \approx\frac{\Omega_{\star}R_{\mathrm{pc}}}{c} \approx3 \times10^{-6}P_{0}^{-3/2}, \tag*{(3)} β rot ≈ c Ω ⋆ R pc ≈ 3 × 1 0 − 6 P 0 − 3/2 , (3) twisting the open field-line tube from the stellar surface. The rotational velocity v \boldsymbol{v} v E = − v × B / c \boldsymbol{E}=-\boldsymbol{v}\times\boldsymbol{B}/c E = − v × B / c
E θ ∼ − Ω ⋆ R p c B r c ≈ − β r o t B ⋆ , (4) E_{\theta} \sim-\frac{\Omega_{\star}R_{\mathrm{pc}}B_{r}}{c} \approx-\beta_{\mathrm{rot}}B_{\star}, \tag*{(4)} E θ ∼ − c Ω ⋆ R pc B r ≈ − β rot B ⋆ , (4) which defines the characteristic rotational field scale E r o t = β r o t B ⋆ E_{\mathrm{rot}}=\beta_{\mathrm{rot}}B_{\star} E rot = β rot B ⋆ ℓ \ell ℓ 1
Δ V = − ∫ E ⋅ d ℓ ∼ β r o t R p c B ⋆ 2 ∼ 7 × 10 12 B 12 P 0 − 2 V . \Delta V=-\int\boldsymbol{E}\cdot\mathrm{d}\boldsymbol{\ell}\sim\frac{\beta_{\mathrm{rot}}R_{\mathrm{pc}}B_{\star}}{2}\sim7\times10^{12}B_{12}P_{0}^{-2}\ \mathrm{V}. Δ V = − ∫ E ⋅ d ℓ ∼ 2 β rot R pc B ⋆ ∼ 7 × 1 0 12 B 12 P 0 − 2 V . When the gap above the polar cap cannot supply the charge density required to maintain E θ E_{\theta} E θ ∇ ⋅ E ≠ 0 \nabla\cdot\boldsymbol{E}\ne0 ∇ ⋅ E = 0 B r ≫ B θ B_{r}\gg B_{\theta} B r ≫ B θ E ∥ E_{\parallel} E ∥
Global Circuit
The rotating magnetosphere forms an electric circuit with two competing elements (see Figure 1 ): rotational electric field and twist current of the field lines. These two sources do not balance in general, resulting in magnetospheric regions with unscreened (parallel) electric fields, known as gaps.
First, rotation of a magnetized conductor induces a surface charge distribution that supports the co-rotation electric field. This distribution determines the sign of charges extracted from the polar caps.2 ∣ n c o ∣ = e n c o \lvert n_{\mathrm{co}}\rvert=e n_{\mathrm{co}} ∣ n co ∣ = e n co 3
n c o = ∣ ∇ ⋅ E ∣ 4 π e ∼ Ω ⋆ B d 2 π c e ≈ β r o t B ⋆ 4 π e R p c ∼ 10 10 B 12 P 0 − 1 c m − 3 , (5) n_{\mathrm{co}}=\frac{\lvert\nabla\cdot\boldsymbol{E}\rvert}{4\pi e}\sim\frac{\Omega_{\star}B_{\mathrm{d}}}{2\pi ce}\approx\frac{\beta_{\mathrm{rot}}B_{\star}}{4\pi eR_{\mathrm{pc}}}\sim10^{10}B_{12}P_{0}^{-1}\ \mathrm{cm}^{-3}, \tag*{(5)} n co = 4 π e ∣ ∇ ⋅ E ∣ ∼ 2 π ce Ω ⋆ B d ≈ 4 π e R pc β rot B ⋆ ∼ 1 0 10 B 12 P 0 − 1 cm − 3 , (5) where B d = B ⋆ / 2 B_{\mathrm{d}}=B_{\star}/2 B d = B ⋆ /2 ω p , c o = 4 π e 2 n c o / m e \omega_{\mathrm{p,co}}=\sqrt{4\pi e^{2}n_{\mathrm{co}}/m_{\mathrm{e}}} ω p , co = 4 π e 2 n co / m e ω p , c o 2 ≈ ω B Ω ⋆ \omega_{\mathrm{p,co}}^{2}\approx\omega_{B}\Omega_{\star} ω p , co 2 ≈ ω B Ω ⋆ ω B ≡ e B ⋆ / ( m e c ) \omega_{B}\equiv eB_{\star}/(m_{\mathrm{e}}c) ω B ≡ e B ⋆ / ( m e c ) Ω ⋆ = 2 π / P ≈ β r o t c / R p c \Omega_{\star}=2\pi/P\approx\beta_{\mathrm{rot}}c/R_{\mathrm{pc}} Ω ⋆ = 2 π / P ≈ β rot c / R pc c / ω p , c o ∼ 3 c m c/\omega_{\mathrm{p,co}}\sim3\ \mathrm{cm} c / ω p , co ∼ 3 cm ω p , c o = σ c o Ω ⋆ \omega_{\mathrm{p,co}}=\sqrt{\sigma_{\mathrm{co}}}\Omega_{\star} ω p , co = σ co Ω ⋆
σ c o = B ⋆ 2 4 π n c o m e c 2 ∼ 10 18 B 12 P 0 , (6) \sigma_{\mathrm{co}}=\frac{B_{\star}^{2}}{4\pi n_{\mathrm{co}}m_{\mathrm{e}}c^{2}}\sim10^{18}B_{12}P_{0}, \tag*{(6)} σ co = 4 π n co m e c 2 B ⋆ 2 ∼ 1 0 18 B 12 P 0 , (6) confirms that the magnetospheric plasma at the co-rotation density n c o n_{\mathrm{co}} n co
Second, the twisted open field-line tube carries an electric current j m = ( c / 4 π ) ∇ × B j_{\mathrm{m}}=(c/4\pi)\nabla\times\boldsymbol{B} j m = ( c /4 π ) ∇ × B [32 ] ; [43 ] ; [5 ] ; [49 ] ) compared to the maximum current sustained by the co-rotation density j c o ≡ ± e n c o c j_{\rm co} \equiv\pm en_{\rm co}c j co ≡ ± e n co c
The parallel electric field along the B B B h h h h = 0 h = 0 h = 0
∂ h E ∥ ≈ 4 π ( η − η c o ) , (7) \partial_h E_{\parallel} \approx4\pi(\eta-\eta_{\rm co}), \tag*{(7)} ∂ h E ∥ ≈ 4 π ( η − η co ) , (7) ∂ t E ∥ ≈ 4 π ( j m − j ± ) . (8) \partial_t E_{\parallel} \approx4\pi(j_{\rm m}-j_{\pm}). \tag*{(8)} ∂ t E ∥ ≈ 4 π ( j m − j ± ) . (8) The field vanishes, E ∥ ( h , t ) = 0 E_{\parallel}(h,t) = 0 E ∥ ( h , t ) = 0 η = η c o \eta= \eta_{\rm co} η = η co j m = j ± = η c β ∥ = η c o c β ∥ = j c o β ∥ j_{\rm m} = j_{\pm} = \eta c\beta_{\parallel} = \eta_{\rm co}c\beta_{\parallel} = j_{\rm co}\beta_{\parallel} j m = j ± = η c β ∥ = η co c β ∥ = j co β ∥ η \eta η j ± j_{\pm} j ± β ∥ \beta_{\parallel} β ∥ B B B E ∥ = 0 E_{\parallel} = 0 E ∥ = 0 ∣ α ∣ < 1 |\alpha| < 1 ∣ α ∣ < 1 ∣ β ∥ ∣ < 1 |\beta_{\parallel}| < 1 ∣ β ∥ ∣ < 1 β ∥ ( h = 0 ) = 0 \beta_{\parallel}(h = 0) = 0 β ∥ ( h = 0 ) = 0 β ∥ \beta_{\parallel} β ∥ α \alpha α
i) Sub-corotating flow. For 0 < α < 1 0 < \alpha< 1 0 < α < 1 p ∥ ∈ [ 0 , p max ] p_{\parallel} \in[0,p_{\max}] p ∥ ∈ [ 0 , p m a x ] p max = 2 α / ( 1 − α 2 ) p_{\max} = 2\alpha/(1-\alpha^2) p m a x = 2 α / ( 1 − α 2 ) [5 ] . While the stationary solution is purely spatially oscillatory, time-dependent kinetic simulations (with free particle supply from the star) show that it is unstable: phase-space mixing develops and the flow evolves into a chaotic trapped-particle state [49 ] . For typical pulsar parameters, p max p_{\max} p m a x
ii) Super-corotating flow. If α > 1 \alpha> 1 α > 1 ∣ η / η c o ∣ > 1 |\eta/\eta_{\rm co}| > 1 ∣ η / η co ∣ > 1
iii) Anti-corotating flow. For α < 0 \alpha< 0 α < 0 α > 1 \alpha> 1 α > 1
Atmospheric Boundary
The star supplies plasma for the circuit. The atmosphere is an ample reservoir of electrons and protons, easily extracted by the unscreened parallel electric field. The atmospheric charges screen the electric field close to the star so that E ∥ ( h < 0 ) ≈ 0 E_{\parallel}(h < 0) \approx0 E ∥ ( h < 0 ) ≈ 0 α > 0 \alpha> 0 α > 0
The atmospheric plasma forms a thin, stratified layer on top of the star’s surface with a scale height H p = p / g ρ H_p = p/g\rho H p = p / g ρ p = p ( ρ ) p = p(\rho) p = p ( ρ ) ρ \rho ρ g g g H p ≈ 6 × 10 − 2 c m H_p \approx6 \times10^{-2}\ {\rm cm} H p ≈ 6 × 1 0 − 2 cm ρ = 10 − 2 g c m − 3 \rho= 10^{-2}\ {\rm g\,cm^{-3}} ρ = 1 0 − 2 g c m − 3 T ≈ 10 5 K T \approx10^5\ {\rm K} T ≈ 1 0 5 K [34 ] . The density drops to n c o n_{\rm co} n co 25 H p ∼ 1.5 c m 25H_p \sim1.5\ {\rm cm} 25 H p ∼ 1.5 cm c / ω p , c o c/\omega_{\rm p,co} c / ω p , co
Gaps with unlimited atmospheric charge have been studied under the Space-Charge-Limited Flow scenario (SCLF; [3 ] ). The opposite case with no charge supply [40 ] is relevant, e.g., to magnetospheres above solid surfaces or around black holes. These historical models have looked for steady-state gap solutions whereas we study time-dependent and intermittent gaps. While there is conceptual difference between these models, we find it irrelevant for the gap dynamics, except that without charge supply pair cascades are ignited even when 0 < α < 1 0 < \alpha< 1 0 < α < 1 [48 ] ; [49 ] ).
QED Processes
QED processes enable photons and pair plasma to interact. For moderate b b b [28 ] ; [29 ] ). The quantum parameter χ \chi χ 4
χ = ∣ p μ G μ ν ∣ B Q , (9) \chi= \frac{|p_{\mu}G^{\mu\nu}|}{B_Q}, \tag*{(9)} χ = B Q ∣ p μ G μν ∣ , (9) where p μ p^{\mu} p μ p μ = ( γ , γ β ) p^{\mu} = (\gamma,\gamma\boldsymbol{\beta}) p μ = ( γ , γ β ) ( x , x ω ^ ) (x,x\hat{\boldsymbol{\omega}}) ( x , x ω ^ ) x x x m e c 2 m_{\rm e}c^2 m e c 2 G μ ν G^{\mu\nu} G μν
χ ± = γ B Q ( B e f f ± ) 2 − ( β ⋅ B ) 2 ≈ γ b sin ψ , (10) \chi_{\pm} = \frac{\gamma}{B_Q}\sqrt{(B_{\rm eff}^{\pm})^2-(\boldsymbol{\beta}\cdot\mathbf{B})^2} \approx\gamma b\sin\psi, \tag*{(10)} χ ± = B Q γ ( B eff ± ) 2 − ( β ⋅ B ) 2 ≈ γ b sin ψ , (10) where B e f f ± ≡ B − β × E \mathbf{B}_{\rm eff}^{\pm} \equiv\mathbf{B}-\boldsymbol{\beta}\times\mathbf{E} B eff ± ≡ B − β × E ψ \psi ψ B \mathbf{B} B (10) corresponds to the ultra-relativistic limit γ ≫ 1 \gamma\gg1 γ ≫ 1
χ γ = x B Q ( B e f f γ ) 2 − ( ω ^ ⋅ B ) 2 ≈ x b sin ψ , (11) \chi_{\gamma} = \frac{x}{B_Q}\sqrt{(B_{\rm eff}^{\gamma})^2-(\hat{\boldsymbol{\omega}}\cdot\mathbf{B})^2} \approx xb\sin\psi, \tag*{(11)} χ γ = B Q x ( B eff γ ) 2 − ( ω ^ ⋅ B ) 2 ≈ x b sin ψ , (11) where B e f f γ ≡ B − ω ^ × E \mathbf{B}_{\rm eff}^{\gamma} \equiv\mathbf{B}-\hat{\boldsymbol{\omega}}\times\mathbf{E} B eff γ ≡ B − ω ^ × E ω ^ \hat{\boldsymbol{\omega}} ω ^
Prominent first-order QED processes include synchrotron emission (also known as nonlinear Compton; e ± → e ± + γ e^{\pm}\to e^{\pm}+\gamma e ± → e ± + γ γ → e + + e − \gamma\to e^+ + e^- γ → e + + e −
n ˙ 1 → 2 ≡ d n d t ∝ α e m b sin ψ t C n 1 T ( χ ) (12) \dot{n}_{1\to2} \equiv\frac{dn}{dt} \propto\frac{\alpha_{\rm em}b\sin\psi}{t_{\rm C}}n_1T(\chi) \tag*{(12)} n ˙ 1 → 2 ≡ d t d n ∝ t C α em b sin ψ n 1 T ( χ ) (12) for an interaction with one incoming and two outgoing particles. Here α e m ≡ e 2 / ( ℏ c ) ≈ 1 / 137 \alpha_{\rm em} \equiv e^2/(\hbar c) \approx1/137 α em ≡ e 2 / ( ℏ c ) ≈ 1/137 t C ≡ λ C / c t_{\rm C} \equiv\lambda_{\rm C}/c t C ≡ λ C / c λ C ≡ ℏ / ( m e c ) \lambda_{\rm C} \equiv\hbar/(m_{\rm e}c) λ C ≡ ℏ/ ( m e c ) T ( χ ) T(\chi) T ( χ ) T s y n T_{\rm syn} T syn T B W T_{\rm BW} T BW n 1 n_1 n 1
Photon Emission
The energy-differential production rate of the synchrotron photons from an e ± e^\pm e ± γ \gamma γ [27 , 8 ] )
d 2 N x d t d x = 2 3 α e m t C χ ± 2 S ( γ , x ) x . (13) \frac{\mathrm{d}^{2}N_x}{\mathrm{d}t\,\mathrm{d}x}=\frac{2}{3}\frac{\alpha_{\mathrm{em}}}{t_{\mathrm{C}}}\chi_{\pm}^{2}\frac{S(\gamma,x)}{x}. \tag*{(13)} d t d x d 2 N x = 3 2 t C α em χ ± 2 x S ( γ , x ) . (13) where S ( γ , x ) S(\gamma,x) S ( γ , x )
P s y n ≡ m e c 2 ∫ 0 ∞ x d 2 N x d t d x d x = 2 3 α e m t C χ ± 2 g Q ( χ ± ) m e c 2 . (14) P_{\mathrm{syn}}\equiv m_{\mathrm{e}}c^{2}\int_{0}^{\infty}x\frac{\mathrm{d}^{2}N_x}{\mathrm{d}t\,\mathrm{d}x}\,\mathrm{d}x=\frac{2}{3}\frac{\alpha_{\mathrm{em}}}{t_{\mathrm{C}}}\chi_{\pm}^{2}g_{\mathrm{Q}}(\chi_{\pm})m_{\mathrm{e}}c^{2}. \tag*{(14)} P syn ≡ m e c 2 ∫ 0 ∞ x d t d x d 2 N x d x = 3 2 t C α em χ ± 2 g Q ( χ ± ) m e c 2 . (14) A Gaunt factor g Q ( χ ± ) ≡ ∫ S d x ≤ 1 g_{\mathrm{Q}}(\chi_{\pm})\equiv\int S\,\mathrm{d}x\leq1 g Q ( χ ± ) ≡ ∫ S d x ≤ 1 S S S g Q → 1 g_{\mathrm{Q}}\to1 g Q → 1 χ ± ≪ 1 \chi_{\pm}\ll1 χ ± ≪ 1 g Q = 1 g_{\mathrm{Q}}=1 g Q = 1
We obtain the differential optical depth (emission rate) for a single particle to emit a synchrotron photon, using the relation χ ± / χ x = γ / x \chi_{\pm}/\chi_x=\gamma/x χ ± / χ x = γ / x
d τ s y n d t = ∫ 0 χ ± d 2 N x d t d x d χ x = α e m b sin ψ t C T s y n ( χ ± ) . (15) \frac{\mathrm{d}\tau_{\mathrm{syn}}}{\mathrm{d}t}=\int_{0}^{\chi_{\pm}}\frac{\mathrm{d}^{2}N_x}{\mathrm{d}t\,\mathrm{d}x}\,\mathrm{d}\chi_x=\frac{\alpha_{\mathrm{em}}b\sin\psi}{t_{\mathrm{C}}}T_{\mathrm{syn}}(\chi_{\pm}). \tag*{(15)} d t d τ syn = ∫ 0 χ ± d t d x d 2 N x d χ x = t C α em b sin ψ T syn ( χ ± ) . (15) For χ ± ≪ 1 \chi_{\pm}\ll1 χ ± ≪ 1 T s y n ( χ ± ) → 5 / ( 2 3 ) ≈ 1.44 T_{\mathrm{syn}}(\chi_{\pm})\to5/(2\sqrt{3})\approx1.44 T syn ( χ ± ) → 5/ ( 2 3 ) ≈ 1.44 T s y n T_{\mathrm{syn}} T syn ∼ 1 % \sim1\% ∼ 1% χ ± ≲ 30 \chi_{\pm}\lesssim30 χ ± ≲ 30
T s y n ( χ ± ) ≈ 5 2 3 ( 1 + 16 3 5 χ ± + 1.1 χ ± 2 ) − 1 / 6 . (16) T_{\mathrm{syn}}(\chi_{\pm})\approx\frac{5}{2\sqrt{3}}\left(1+\frac{16\sqrt{3}}{5}\chi_{\pm}+1.1\chi_{\pm}^{2}\right)^{-1/6}. \tag*{(16)} T syn ( χ ± ) ≈ 2 3 5 ( 1 + 5 16 3 χ ± + 1.1 χ ± 2 ) − 1/6 . (16) The linear coefficient 16 3 / 5 ≈ 5.5 16\sqrt{3}/5\approx5.5 16 3 /5 ≈ 5.5 T s y n = 5 2 3 ( 1 − 8 3 15 χ ± ) + O ( χ ± 2 ) T_{\mathrm{syn}}=\frac{5}{2\sqrt{3}}(1-\frac{8\sqrt{3}}{15}\chi_{\pm})+O(\chi_{\pm}^{2}) T syn = 2 3 5 ( 1 − 15 8 3 χ ± ) + O ( χ ± 2 )
The synchrotron emission process reduces to curvature radiation by substituting χ ± = γ b sin ψ \chi_{\pm}=\gamma b\sin\psi χ ± = γ b sin ψ P s y n P_{\mathrm{syn}} P syn sin ψ → γ r g / R c u r v \sin\psi\to\gamma r_g/R_{\mathrm{curv}} sin ψ → γ r g / R curv r g = c / ω B r_g=c/\omega_B r g = c / ω B
P c u r v m e c 2 = 2 3 α e m λ C b 2 r g 2 R c u r v 2 γ 4 . (17) \frac{P_{\mathrm{curv}}}{m_{\mathrm{e}}c^{2}}=\frac{2}{3}\frac{\alpha_{\mathrm{em}}}{\lambda_{\mathrm{C}}}b^{2}\frac{r_g^{2}}{R_{\mathrm{curv}}^{2}}\gamma^{4}. \tag*{(17)} m e c 2 P curv = 3 2 λ C α em b 2 R curv 2 r g 2 γ 4 . (17) Pair Production
The differential optical depth for a gamma ray to pair-create via the one-photon Breit–Wheeler process is
d τ B W d t = α e m b sin ψ t C T B W ( χ x ) , (18) \frac{\mathrm{d}\tau_{\mathrm{BW}}}{\mathrm{d}t}=\frac{\alpha_{\mathrm{em}}b\sin\psi}{t_{\mathrm{C}}}T_{\mathrm{BW}}(\chi_x), \tag*{(18)} d t d τ BW = t C α em b sin ψ T BW ( χ x ) , (18) where T B W ( χ x ) → 0.23 exp ( − 8 / ( 3 χ x ) ) T_{\mathrm{BW}}(\chi_x)\to0.23\exp(-8/(3\chi_x)) T BW ( χ x ) → 0.23 exp ( − 8/ ( 3 χ x )) χ x ≪ 1 \chi_x\ll1 χ x ≪ 1 [17 ] ) and T B W ( χ x ) → 0.38 χ x − 1 / 3 T_{\mathrm{BW}}(\chi_x)\to0.38\chi_x^{-1/3} T BW ( χ x ) → 0.38 χ x − 1/3 χ x ≫ 1 \chi_x\gg1 χ x ≫ 1 [27 ] ). We combine these limits and approximate T B W T_{\mathrm{BW}} T BW ∼ 10 % \sim10\% ∼ 10% χ x ≲ 100 \chi_x\lesssim100 χ x ≲ 100
T B W ( χ x ) ≈ 0.23 exp ( − 8 3 χ x ) ( 1 + 0.22 χ x ) − 1 / 3 . (19) T_{\mathrm{BW}}(\chi_x)\approx0.23\exp\left(-\frac{8}{3\chi_x}\right)(1+0.22\chi_x)^{-1/3}. \tag*{(19)} T BW ( χ x ) ≈ 0.23 exp ( − 3 χ x 8 ) ( 1 + 0.22 χ x ) − 1/3 . (19) where the coefficient 0.22 ≈ ( 0.23 / 0.38 ) 3 0.22\approx(0.23/0.38)^3 0.22 ≈ ( 0.23/0.38 ) 3
The corresponding mean free path l m f p l_{\mathrm{mfp}} l mfp l m f p ≈ R c u r v sin ψ a l_{\mathrm{mfp}}\approx R_{\mathrm{curv}}\sin\psi_a l mfp ≈ R curv sin ψ a ψ a \psi_a ψ a l m f p l_{\mathrm{mfp}} l mfp χ x ≈ x b sin ψ \chi_x\approx xb\sin\psi χ x ≈ x b sin ψ l m f p ≈ R c u r v χ x a / x b l_{\mathrm{mfp}}\approx R_{\mathrm{curv}}\chi_x^a/xb l mfp ≈ R curv χ x a / x b χ x a \chi_x^a χ x a τ B W ( χ x a ) = 1 \tau_{\mathrm{BW}}(\chi_x^a)=1 τ BW ( χ x a ) = 1 χ x ≪ 1 \chi_x\ll1 χ x ≪ 1 5
l m f p ≈ 8 R c u r v 9 x b W [ 2 ( λ C x a b 2 α e m R c u r v ) − 1 / 3 − 1 ] − 1 ≈ 8 R c u r v 9 x b { ln [ ( α e m R c u r v 16 λ C x a b 2 ) 1 / 3 ] − ln ln [ ( α e m R c u r v 16 λ C x a b 2 ) 1 / 3 ] } − 1 , (20) \begin{aligned}
l_{\mathrm{mfp}}\approx\frac{8R_{\mathrm{curv}}}{9xb}W\left[\frac{2}{\left(\frac{\lambda_{\mathrm{C}}x^{a}b^{2}}{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}\right)^{-1/3}-1}\right]^{-1} \\
&\approx\frac{8R_{\mathrm{curv}}}{9xb}\left\{\ln\left[\left(\frac{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}{16\lambda_{\mathrm{C}}x^{a}b^{2}}\right)^{1/3}\right]-\ln\ln\left[\left(\frac{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}{16\lambda_{\mathrm{C}}x^{a}b^{2}}\right)^{1/3}\right]\right\}^{-1}, \tag*{(20)}
\end{aligned} l mfp ≈ 9 x b 8 R curv W ( α em R curv λ C x a b 2 ) − 1/3 − 1 2 − 1 ≈ 9 x b 8 R curv { ln [ ( 16 λ C x a b 2 α em R curv ) 1/3 ] − ln ln [ ( 16 λ C x a b 2 α em R curv ) 1/3 ] } − 1 , (20) where W W W W W W z z z W ( z ) ≈ ln z − ln ln z W(z)\approx\ln z-\ln\ln z W ( z ) ≈ ln z − ln ln z l m f p ∼ 10 − 2 R p c l_{\mathrm{mfp}}\sim10^{-2}R_{\mathrm{pc}} l mfp ∼ 1 0 − 2 R pc x ∼ x c u r v x\sim x_{\mathrm{curv}} x ∼ x curv (25) ).
Neglected Processes
We ignore other first-order processes besides the synchrotron and one-photon pair creation. These include, e.g., synchrotron absorption (e ± + γ → e ± e^\pm+\gamma\to e^\pm e ± + γ → e ± e − + e + → γ e^-+e^+\to\gamma e − + e + → γ σ a n n ≈ 0 \sigma_{\mathrm{ann}}\approx0 σ ann ≈ 0 b ≲ 1 b\lesssim1 b ≲ 1
We also ignore all the higher-order processes. These include, e.g., QED processes with two vertices and 2 + 2 2+2 2 + 2 e ± + γ → e ± + γ e^\pm+\gamma\to e^\pm+\gamma e ± + γ → e ± + γ γ + γ → e − + e + \gamma+\gamma\to e^-+e^+ γ + γ → e − + e + e − + e + → γ + γ ′ e^-+e^+\to\gamma+\gamma' e − + e + → γ + γ ′ b ≪ 1 b\ll1 b ≪ 1 [41 ] ). Other effects, such as photon polarization, if included, enhance the pair yield in the cascade by a few percent ([44 ] ).
Particle Dynamics in the Gap
We consider individual particle dynamics in the gap. The time evolution of the proper velocity (the spatial part of the four-velocity) u = γ β c \boldsymbol{u}=\gamma\boldsymbol{\beta}c u = γ β c
d u d t = q m e ( E + v c × B ) − F r a d m e , (21) \frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t}=\frac{q}{m_{\mathrm{e}}}\left(\boldsymbol{E}+\frac{\boldsymbol{v}}{c}\times\boldsymbol{B}\right)-\frac{\boldsymbol{F}_{\mathrm{rad}}}{m_{\mathrm{e}}}, \tag*{(21)} d t d u = m e q ( E + c v × B ) − m e F rad , (21) where q = − e q=-e q = − e q = + e q=+e q = + e e > 0 e>0 e > 0 F r a d \boldsymbol{F}_{\mathrm{rad}} F rad
τ B W ( χ x a ) = ∫ 0 l m f p α e m t C χ x x T B W ( χ x ) d l ≈ 0.23 α e m λ C ∫ 0 l m f p χ x x e − 8 / ( 3 χ x ) d l \tau_{\mathrm{BW}}(\chi_x^a)=\int_{0}^{l_{\mathrm{mfp}}}\frac{\alpha_{\mathrm{em}}}{t_{\mathrm{C}}}\frac{\chi_x}{x}T_{\mathrm{BW}}(\chi_x)\,\mathrm{d}l\approx0.23\frac{\alpha_{\mathrm{em}}}{\lambda_{\mathrm{C}}}\int_{0}^{l_{\mathrm{mfp}}}\frac{\chi_x}{x}e^{-8/(3\chi_x)}\,\mathrm{d}l τ BW ( χ x a ) = ∫ 0 l mfp t C α em x χ x T BW ( χ x ) d l ≈ 0.23 λ C α em ∫ 0 l mfp x χ x e − 8/ ( 3 χ x ) d l ≈ 0.23 α e m R c u r v λ C x b 2 ∫ 0 χ x a χ x e − 8 / ( 3 χ x ) d χ x = 0.23 α e m R c u r v λ C x b 2 ( χ x a ) 2 E 3 ( 8 3 χ x a ) \approx0.23\frac{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}{\lambda_{\mathrm{C}}xb^{2}}\int_{0}^{\chi_x^a}\chi_xe^{-8/(3\chi_x)}\,\mathrm{d}\chi_x=0.23\frac{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}{\lambda_{\mathrm{C}}xb^{2}}(\chi_x^a)^2E_3\left(\frac{8}{3\chi_x^a}\right) ≈ 0.23 λ C x b 2 α em R curv ∫ 0 χ x a χ x e − 8/ ( 3 χ x ) d χ x = 0.23 λ C x b 2 α em R curv ( χ x a ) 2 E 3 ( 3 χ x a 8 ) ≈ 0.086 α e m R c u r v λ C x b 2 ( χ x a ) 3 e − 8 / ( 3 χ x a ) , \approx0.086\frac{\alpha_{\mathrm{em}}R_{\mathrm{curv}}}{\lambda_{\mathrm{C}}xb^{2}}(\chi_x^a)^3e^{-8/(3\chi_x^a)}, ≈ 0.086 λ C x b 2 α em R curv ( χ x a ) 3 e − 8/ ( 3 χ x a ) , where E 3 E_3 E 3 z z z z 2 E 3 ( 8 / 3 z ) ≈ 8 3 z 3 e − 8 / ( 3 z ) + O ( z 4 e − 8 / ( 3 z ) ) z^2E_3(8/3z)\approx\frac{8}{3}z^3e^{-8/(3z)}+O(z^4e^{-8/(3z)}) z 2 E 3 ( 8/3 z ) ≈ 3 8 z 3 e − 8/ ( 3 z ) + O ( z 4 e − 8/ ( 3 z ) ) h h h B \mathbf{B} B u ≡ β γ c → γ c u \equiv\beta\gamma c \to\gamma c u ≡ β γ c → γ c P r a d = F r a d ⋅ β c → F r a d c P_{\rm rad} = \mathbf{F}_{\rm rad} \cdot\boldsymbol{\beta}c \to F_{\rm rad}c P rad = F rad ⋅ β c → F rad c
d γ d t = q E ∥ m e c − P r a d m e c 2 . (22) \frac{\mathrm{d}\gamma}{\mathrm{d}t} = \frac{qE_{\parallel}}{m_{\mathrm{e}}c} - \frac{P_{\rm rad}}{m_{\mathrm{e}}c^{2}} . \tag*{(22)} d t d γ = m e c q E ∥ − m e c 2 P rad . (22) The parallel electric field at the polar cap can be approximated with a linear-voltage gap (e.g., [16 ] ; [48 ] )
E ∥ = { − Δ V H g a p , h ≤ H g a p 0 , h > H g a p , (23) E_{\parallel} =
\begin{cases}
-\dfrac{\Delta V}{H_{\mathrm{gap}}}, & h \le H_{\mathrm{gap}} \\
0, & h > H_{\mathrm{gap}},
\end{cases} \tag*{(23)} E ∥ = ⎩ ⎨ ⎧ − H gap Δ V , 0 , h ≤ H gap h > H gap , (23) where H g a p H_{\mathrm{gap}} H gap h = 0 h = 0 h = 0 H g a p = R p c H_{\mathrm{gap}} = R_{\mathrm{pc}} H gap = R pc
Curvature radiation provides the radiative force with an emission power
d γ d t ∣ c u r v = − P c u r v m e c 2 . (24) \left.\frac{\mathrm{d}\gamma}{\mathrm{d}t}\right|_{\mathrm{curv}} = -\frac{P_{\mathrm{curv}}}{m_{\mathrm{e}}c^{2}} . \tag*{(24)} d t d γ curv = − m e c 2 P curv . (24) The characteristic photon energy6
x c u r v ≡ ℏ ω c u r v m e c 2 = 3 2 b r g R c u r v γ 3 , (25) x_{\mathrm{curv}} \equiv\frac{\hbar\omega_{\mathrm{curv}}}{m_{\mathrm{e}}c^{2}} = \frac{3}{2}b\frac{r_{\mathrm{g}}}{R_{\mathrm{curv}}}\gamma^{3}, \tag*{(25)} x curv ≡ m e c 2 ℏ ω curv = 2 3 b R curv r g γ 3 , (25) where r g ∼ 10 − 9 c m r_{\mathrm{g}} \sim10^{-9}\,\mathrm{cm} r g ∼ 1 0 − 9 cm r g / R c u r v ∼ 10 − 15 r_{\mathrm{g}}/R_{\mathrm{curv}} \sim10^{-15} r g / R curv ∼ 1 0 − 15
The particle dynamics in terms of distance from the surface h h h d h ≈ c d t \mathrm{d}h \approx c\mathrm{d}t d h ≈ c d t
d γ d h = q E ∥ m e c 2 − F r a d m e c 2 ≈ − q Δ V m e c 2 H g a p − 2 3 α e m λ C r g 2 R c u r v 2 b 2 γ 4 . (26) \frac{\mathrm{d}\gamma}{\mathrm{d}h} = \frac{qE_{\parallel}}{m_{\mathrm{e}}c^{2}} - \frac{F_{\mathrm{rad}}}{m_{\mathrm{e}}c^{2}} \approx-\frac{q\Delta V}{m_{\mathrm{e}}c^{2}H_{\mathrm{gap}}} - \frac{2}{3}\frac{\alpha_{\mathrm{em}}}{\lambda_{\mathrm{C}}}\frac{r_{\mathrm{g}}^{2}}{R_{\mathrm{curv}}^{2}}b^{2}\gamma^{4}. \tag*{(26)} d h d γ = m e c 2 q E ∥ − m e c 2 F rad ≈ − m e c 2 H gap q Δ V − 3 2 λ C α em R curv 2 r g 2 b 2 γ 4 . (26) The maximum Lorentz factor of an electron experiencing the full potential drop across the polar cap (i.e., h h h H g a p H_{\mathrm{gap}} H gap F r a d → 0 F_{\mathrm{rad}} \to0 F rad → 0
γ g a p = ∫ 0 H g a p d γ d h d h = e Δ V m e c 2 ≈ e B ⋆ β r o t R p c 2 m e c 2 . (27) \gamma_{\mathrm{gap}} = \int_{0}^{H_{\mathrm{gap}}}\frac{\mathrm{d}\gamma}{\mathrm{d}h}\,\mathrm{d}h = \frac{e\Delta V}{m_{\mathrm{e}}c^{2}} \approx\frac{eB_{\star}\beta_{\mathrm{rot}}R_{\mathrm{pc}}}{2m_{\mathrm{e}}c^{2}} . \tag*{(27)} γ gap = ∫ 0 H gap d h d γ d h = m e c 2 e Δ V ≈ 2 m e c 2 e B ⋆ β rot R pc . (27) We have γ g a p ∼ 10 7 \gamma_{\mathrm{gap}} \sim10^{7} γ gap ∼ 1 0 7
The steady-state Lorentz factor (d γ / d t = 0 \mathrm{d}\gamma/\mathrm{d}t = 0 d γ / d t = 0 [16 ] ),
γ r a d = ( 3 e Δ V B Q 2 R c u r v 2 2 m e c 2 B ⋆ 2 r g 2 λ C α e m H g a p ) 1 / 4 = γ g a p 1 / 4 b − 1 / 2 ( 3 λ C 2 α e m H g a p ) 1 / 4 ( r g R c u r v ) − 1 / 2 . (28) \begin{aligned}
\gamma_{\mathrm{rad}} = \left(\frac{3e\Delta V B_{Q}^{2}R_{\mathrm{curv}}^{2}}{2m_{\mathrm{e}}c^{2}B_{\star}^{2}r_{\mathrm{g}}^{2}}\frac{\lambda_{\mathrm{C}}}{\alpha_{\mathrm{em}}H_{\mathrm{gap}}}\right)^{1/4} \\
&= \gamma_{\mathrm{gap}}^{1/4}b^{-1/2}\left(\frac{3\lambda_{\mathrm{C}}}{2\alpha_{\mathrm{em}}H_{\mathrm{gap}}}\right)^{1/4}\left(\frac{r_{\mathrm{g}}}{R_{\mathrm{curv}}}\right)^{-1/2}. \tag*{(28)}
\end{aligned} γ rad = ( 2 m e c 2 B ⋆ 2 r g 2 3 e Δ V B Q 2 R curv 2 α em H gap λ C ) 1/4 = γ gap 1/4 b − 1/2 ( 2 α em H gap 3 λ C ) 1/4 ( R curv r g ) − 1/2 . (28) and it is obtained over a distance7 ( γ r a d / γ g a p ) H g a p (\gamma_{\mathrm{rad}}/\gamma_{\mathrm{gap}})H_{\mathrm{gap}} ( γ rad / γ gap ) H gap γ r a d ∼ 10 7 \gamma_{\mathrm{rad}} \sim10^{7} γ rad ∼ 1 0 7 γ r a d / γ g a p ≈ 0.6 \gamma_{\mathrm{rad}}/\gamma_{\mathrm{gap}} \approx0.6 γ rad / γ gap ≈ 0.6
The accelerated pairs emit curvature photons with a characteristic energy x c u r v x_{\mathrm{curv}} x curv x ≳ x t h r ≃ 2 x \gtrsim x_{\mathrm{thr}} \simeq2 x ≳ x thr ≃ 2
γ t h r ∼ ( 4 R c u r v 3 b r g ) 1 / 3 , (29) \gamma_{\mathrm{thr}} \sim\left(\frac{4R_{\mathrm{curv}}}{3br_{\mathrm{g}}}\right)^{1/3}, \tag*{(29)} γ thr ∼ ( 3 b r g 4 R curv ) 1/3 , (29) which is γ t h r ∼ 3 × 10 5 \gamma_{\mathrm{thr}} \sim3 \times10^{5} γ thr ∼ 3 × 1 0 5 γ r a d ≫ γ t h r \gamma_{\mathrm{rad}} \gg\gamma_{\mathrm{thr}} γ rad ≫ γ thr
Since χ x \chi_x χ x sin ψ ≈ l / R c u r v \sin\psi\approx l/R_{\mathrm{curv}} sin ψ ≈ l / R curv l m f p ( x ) ≲ H g a p l_{\mathrm{mfp}}(x) \lesssim H_{\mathrm{gap}} l mfp ( x ) ≲ H gap l m f p ≈ R c u r v χ x a / ( x b ) l_{\mathrm{mfp}} \approx R_{\mathrm{curv}}\chi_x^{a}/(xb) l mfp ≈ R curv χ x a / ( x b ) x min ≈ χ x a R c u r v / ( H g a p b ) ≈ 500 x_{\min} \approx\chi_x^{a}R_{\mathrm{curv}}/(H_{\mathrm{gap}}b) \approx500 x m i n ≈ χ x a R curv / ( H gap b ) ≈ 500 χ x a \chi_x^{a} χ x a x x x
γ min = x min 2 ≈ 250. (30) \gamma_{\min} = \frac{x_{\min}}{2} \approx250. \tag*{(30)} γ m i n = 2 x m i n ≈ 250. (30) We conclude that there is a hierarchy of energy scales
1 ≪ γ min ≪ γ t h r < γ r a d ≲ γ g a p , (31) 1 \ll\gamma_{\min} \ll\gamma_{\mathrm{thr}} < \gamma_{\mathrm{rad}} \lesssim\gamma_{\mathrm{gap}}, \tag*{(31)} 1 ≪ γ m i n ≪ γ thr < γ rad ≲ γ gap , (31) which also needs to be fulfilled in numerical simulations.
Discharge Dynamics
Particles energized in the gap regions can trigger explosive plasma production and intermittent limit-cycle behavior. We construct a concurrent model for the discharge dynamics and plasma generation in the gap regions. See also Appendix B for the full derivation of the analytic expressions.
Characteristic Timescales
The gap height H g a p ∼ R p c H_{\mathrm{gap}} \sim R_{\mathrm{pc}} H gap ∼ R pc
t e s c = H g a p c ∼ 5 × 10 − 7 s . (32) t_{\mathrm{esc}} = \frac{H_{\mathrm{gap}}}{c} \sim5 \times10^{-7}\,\mathrm{s}. \tag*{(32)} t esc = c H gap ∼ 5 × 1 0 − 7 s . (32) Particle acceleration sets the fastest discharge timescale. The energization is rapid,
d γ d t ∼ Δ γ t a c c ∼ q E ∥ m e c , (33) \frac{\mathrm{d}\gamma}{\mathrm{d}t} \sim\frac{\Delta\gamma}{t_{\mathrm{acc}}} \sim\frac{qE_{\parallel}}{m_{\mathrm{e}}c}, \tag*{(33)} d t d γ ∼ t acc Δ γ ∼ m e c q E ∥ , (33) with a timescale
t a c c ∼ ∣ Δ γ m e c e E ∥ ∣ ≈ m e c e β r o t B ⋆ = β r o t − 1 ω B − 1 ∼ 10 − 14 s , (34) t_{\mathrm{acc}} \sim\left|\frac{\Delta\gamma m_{\mathrm{e}}c}{eE_{\parallel}}\right| \approx\frac{m_{\mathrm{e}}c}{e\beta_{\mathrm{rot}}B_{\star}} = \beta_{\mathrm{rot}}^{-1}\omega_{B}^{-1} \sim10^{-14}\,\mathrm{s}, \tag*{(34)} t acc ∼ e E ∥ Δ γ m e c ≈ e β rot B ⋆ m e c = β rot − 1 ω B − 1 ∼ 1 0 − 14 s , (34) where E ∥ ≈ − β r o t B ⋆ E_{\parallel} \approx-\beta_{\mathrm{rot}}B_{\star} E ∥ ≈ − β rot B ⋆ Δ γ = 1 \Delta\gamma= 1 Δ γ = 1 m e c 2 m_{\mathrm{e}}c^{2} m e c 2
Curvature-photon emission and one-photon pair creation proceed on intermediate timescales. High-energy
d l = d γ A − B γ 4 , \mathrm{d}l = \frac{\mathrm{d}\gamma}{A-B\gamma^{4}}, d l = A − B γ 4 d γ , and so
l = ∫ d l = ∫ d γ A − B γ 4 ≈ γ A + B γ 5 5 A 2 + O ( γ 9 ) , (35) l = \int\mathrm{d}l = \int\frac{\mathrm{d}\gamma}{A-B\gamma^{4}} \approx\frac{\gamma}{A} + \frac{B\gamma^{5}}{5A^{2}} + \mathcal{O}(\gamma^{9}), \tag*{(35)} l = ∫ d l = ∫ A − B γ 4 d γ ≈ A γ + 5 A 2 B γ 5 + O ( γ 9 ) , (35) solved via Taylor expansion of the integral (convergent for γ < γ r a d \gamma< \gamma_{\mathrm{rad}} γ < γ rad A ≫ B γ 4 A \gg B\gamma^{4} A ≫ B γ 4 x ≳ 2 x \gtrsim2 x ≳ 2 γ > γ t h r \gamma> \gamma_{\mathrm{thr}} γ > γ thr t ± = l m f p / c ∼ 10 − 8 s t_{\pm}=l_{\mathrm{mfp}}/c\sim10^{-8}\,\mathrm{s} t ± = l mfp / c ∼ 1 0 − 8 s l m f p l_{\mathrm{mfp}} l mfp (20) , evaluated at x = x c u r v ( γ r a d ) x=x_{\mathrm{curv}}(\gamma_{\mathrm{rad}}) x = x curv ( γ rad ) e ± e^{\pm} e ± (15) as8
d τ c u r v d t ∼ t x − 1 ≈ α e m b c λ ˉ C r g R c u r v γ r a d (36) \frac{\mathrm{d}\tau_{\mathrm{curv}}}{\mathrm{d}t}\sim t_x^{-1}\approx\alpha_{\mathrm{em}}b\frac{c}{\bar{\lambda}_{\mathrm{C}}}\frac{r_g}{R_{\mathrm{curv}}}\gamma_{\mathrm{rad}} \tag*{(36)} d t d τ curv ∼ t x − 1 ≈ α em b λ ˉ C c R curv r g γ rad (36) to give
t x ∼ λ ˉ C α e m b c R c u r v r g 1 γ r a d ≈ 5 × 10 − 10 s . (37) t_x\sim\frac{\bar{\lambda}_{\mathrm{C}}}{\alpha_{\mathrm{em}}bc}\frac{R_{\mathrm{curv}}}{r_g}\frac{1}{\gamma_{\mathrm{rad}}}\approx5\times10^{-10}\,\mathrm{s}. \tag*{(37)} t x ∼ α em b c λ ˉ C r g R curv γ rad 1 ≈ 5 × 1 0 − 10 s . (37) Hence, for radio pulsars, t a c c ≪ t x ≲ t ± ≪ t e s c t_{\mathrm{acc}}\ll t_x\lesssim t_{\pm}\ll t_{\mathrm{esc}} t acc ≪ t x ≲ t ± ≪ t esc t x ∼ t ± t_x\sim t_{\pm} t x ∼ t ±
The gap with ∣ E ∥ / E r o t ∣ > 0 \lvert E_{\parallel}/E_{\mathrm{rot}}\rvert>0 ∣ E ∥ / E rot ∣ > 0 j ± = e n ± c j_{\pm}=en_{\pm}c j ± = e n ± c
∂ t E ∥ ∼ E ∥ t s c r e e n ≈ − 4 π j ± , (38) \partial_t E_{\parallel}\sim\frac{E_{\parallel}}{t_{\mathrm{screen}}}\approx-4\pi j_{\pm}, \tag*{(38)} ∂ t E ∥ ∼ t screen E ∥ ≈ − 4 π j ± , (38) where we have ignored j m j_{\mathrm{m}} j m j m ≪ j ± j_{\mathrm{m}}\ll j_{\pm} j m ≪ j ±
t s c r e e n ∼ ∣ E ∥ 4 π j ± ∣ ≈ t e s c M ± ′ , (39) t_{\mathrm{screen}}\sim\left\lvert\frac{E_{\parallel}}{4\pi j_{\pm}}\right\rvert\approx\frac{t_{\mathrm{esc}}}{M'_{\pm}}, \tag*{(39)} t screen ∼ 4 π j ± E ∥ ≈ M ± ′ t esc , (39) where M ± ′ = n ± / n c o ≫ 1 M'_{\pm}=n_{\pm}/n_{\mathrm{co}}\gg1 M ± ′ = n ± / n co ≫ 1
In addition to the discharge timescales, the plasma has two kinetic timescales. The gyro timescale satisfies ω B − 1 ≪ t a c c \omega_B^{-1}\ll t_{\mathrm{acc}} ω B − 1 ≪ t acc B B B ω p , c o − 1 ≲ t s c r e e n \omega_{\mathrm{p,co}}^{-1}\lesssim t_{\mathrm{screen}} ω p , co − 1 ≲ t screen
The discharge dynamics then has a hierarchy of timescales set as
ω B − 1 ≪ t a c c ≪ t s c r e e n ∼ ω p , c o − 1 ≲ t x ≲ t ± ≪ t e s c . (40) \omega_B^{-1}\ll t_{\mathrm{acc}}\ll t_{\mathrm{screen}}\sim\omega_{\mathrm{p,co}}^{-1}\lesssim t_x\lesssim t_{\pm}\ll t_{\mathrm{esc}}. \tag*{(40)} ω B − 1 ≪ t acc ≪ t screen ∼ ω p , co − 1 ≲ t x ≲ t ± ≪ t esc . (40) Discharge Dynamics Model
We express the pair number density n ± n_{\pm} n ± n x n_x n x
d n ± d t = 2 n x t ± − n ± t e s c (41) \frac{\mathrm{d}n_{\pm}}{\mathrm{d}t}=2\frac{n_x}{t_{\pm}}-\frac{n_{\pm}}{t_{\mathrm{esc}}} \tag*{(41)} d t d n ± = 2 t ± n x − t esc n ± (41) d n x d t = − n x t ± + n ± t x H ( ∣ E ∥ E r o t ∣ − 1 ) − n x t e s c , (42) \frac{\mathrm{d}n_x}{\mathrm{d}t}=-\frac{n_x}{t_{\pm}}+\frac{n_{\pm}}{t_x}\mathcal{H}\left(\left\lvert\frac{E_{\parallel}}{E_{\mathrm{rot}}}\right\rvert-1\right)-\frac{n_x}{t_{\mathrm{esc}}}, \tag*{(42)} d t d n x = − t ± n x + t x n ± H ( E rot E ∥ − 1 ) − t esc n x , (42) where H \mathcal{H} H ∣ E ∥ / E r o t ∣ = 1 \lvert E_{\parallel}/E_{\mathrm{rot}}\rvert=1 ∣ E ∥ / E rot ∣ = 1
The electric field in the gap evolves according to Maxwell’s law
d E ∥ d t ≈ 4 π ( j m − j ± ) . (43) \frac{\mathrm{d}E_{\parallel}}{\mathrm{d}t}\approx4\pi(j_{\mathrm{m}}-j_{\pm}). \tag*{(43)} d t d E ∥ ≈ 4 π ( j m − j ± ) . (43) The current from the pair plasma, j ± = ∑ s q s n s β ∥ , s c j_{\pm}=\sum_s q_s n_s\beta_{\parallel,s}c j ± = ∑ s q s n s β ∥ , s c ∣ β ∥ ∣ \lvert\beta_{\parallel}\rvert ∣ β ∥ ∣ ∣ j ± ∣ = e n ± ∣ β ∥ ∣ c \lvert j_{\pm}\rvert=en_{\pm}\lvert\beta_{\parallel}\rvert c ∣ j ± ∣ = e n ± ∣ β ∥ ∣ c
The current dynamics follow from the 1D Vlasov equation for the particle distribution f s ( h , u , t ) f_s(h,u,t) f s ( h , u , t )
∂ f s ∂ t = − c β ∂ f s ∂ h − q s E ∥ m e ∂ f s ∂ u . (44) \frac{\partial f_s}{\partial t}=-c\beta\frac{\partial f_s}{\partial h}-\frac{q_sE_{\parallel}}{m_e}\frac{\partial f_s}{\partial u}. \tag*{(44)} ∂ t ∂ f s = − c β ∂ h ∂ f s − m e q s E ∥ ∂ u ∂ f s . (44) Here we have ignored the source term from the particle generation. The current equation follows from differentiating its definition and substituting the Vlasov equation (see also [51 ] ),
∂ j ± ∂ t = ∑ s q s ∫ c β ∂ f s ∂ t d u = ∑ s ( − q s ∫ c 2 β 2 ∂ f s ∂ h d u + q s 2 m e ∫ c E ∥ f s ∂ β ∂ u d u ) ≈ − c L j ± + e 2 E ∥ n ± m e ⟨ 1 γ 3 ⟩ , (45) \begin{aligned}
\frac{\partial j_{\pm}}{\partial t}=\sum_s q_s\int c\beta\frac{\partial f_s}{\partial t}\,\mathrm{d}u \\
&=\sum_s\left(-q_s\int c^2\beta^2\frac{\partial f_s}{\partial h}\,\mathrm{d}u+\frac{q_s^2}{m_e}\int cE_{\parallel}f_s\frac{\partial\beta}{\partial u}\,\mathrm{d}u\right) \\
&\approx-\frac{c}{L}j_{\pm}+\frac{e^2E_{\parallel}n_{\pm}}{m_e}\left\langle\frac{1}{\gamma^3}\right\rangle, \tag*{(45)}
\end{aligned} ∂ t ∂ j ± = s ∑ q s ∫ c β ∂ t ∂ f s d u = s ∑ ( − q s ∫ c 2 β 2 ∂ h ∂ f s d u + m e q s 2 ∫ c E ∥ f s ∂ u ∂ β d u ) ≈ − L c j ± + m e e 2 E ∥ n ± ⟨ γ 3 1 ⟩ , (45) where n ± = ∫ f d u n_{\pm}=\int f\,\mathrm{d}u n ± = ∫ f d u ⟨ γ − 3 ⟩ ≡ n ± − 1 ∫ f γ − 3 d u \langle\gamma^{-3}\rangle\equiv n_{\pm}^{-1}\int f\gamma^{-3}\,\mathrm{d}u ⟨ γ − 3 ⟩ ≡ n ± − 1 ∫ f γ − 3 d u f = ∑ s f s f=\sum_s f_s f = ∑ s f s ∂ h → c / L \partial_h\to c/L ∂ h → c / L L ∼ H g a p L\sim H_{\mathrm{gap}} L ∼ H gap j ± j_{\pm} j ± β → 1 \beta\to1 β → 1
Next, we express the number densities as n ± → m ± n c o n_{\pm}\to m_{\pm}n_{\mathrm{co}} n ± → m ± n co n x → m x n c o n_x\to m_xn_{\mathrm{co}} n x → m x n co m ± m_{\pm} m ± m x m_x m x j ± → j e n c o c j_{\pm}\to j_{\mathrm{en}_{\mathrm{co}}c} j ± → j en co c ∣ j ∣ = m ± ∣ β ∥ ∣ \lvert j\rvert=m_{\pm}\lvert\beta_{\parallel}\rvert ∣ j ∣ = m ± ∣ β ∥ ∣ j m → α j c o j_{\mathrm{m}}\to\alpha j_{\mathrm{co}} j m → α j co E ∥ → ε ∣ E r o t ∣ E_{\parallel}\to\varepsilon\lvert E_{\mathrm{rot}}\rvert E ∥ → ε ∣ E rot ∣ ε \varepsilon ε
m ˙ ± = 2 m x t ± − m ± t e s c (46) \dot{m}_{\pm}=2\frac{m_x}{t_{\pm}}-\frac{m_{\pm}}{t_{\mathrm{esc}}} \tag*{(46)} m ˙ ± = 2 t ± m x − t esc m ± (46) m ˙ x = − m x t ± + m ± t x H ( ∣ ε ∣ − 1 ) − m x t e s c (47) \dot{m}_x=-\frac{m_x}{t_{\pm}}+\frac{m_{\pm}}{t_x}\mathcal{H}(\lvert\varepsilon\rvert-1)-\frac{m_x}{t_{\mathrm{esc}}} \tag*{(47)} m ˙ x = − t ± m x + t x m ± H (∣ ε ∣ − 1 ) − t esc m x (47) ε ˙ = α j t e s c (48) \dot{\varepsilon}=\frac{\alpha j}{t_{\mathrm{esc}}} \tag*{(48)} ε ˙ = t esc α j (48) j ˙ = − j t e s c + m ± ε t a c c ⟨ 1 γ 3 ⟩ , (49) \dot{j}=-\frac{j}{t_{\mathrm{esc}}}+\frac{m_{\pm}\varepsilon}{t_{\mathrm{acc}}}\left\langle\frac{1}{\gamma^3}\right\rangle, \tag*{(49)} j ˙ = − t esc j + t acc m ± ε ⟨ γ 3 1 ⟩ , (49) where overdots denote time derivatives (Q ˙ ≡ d Q / d t \dot{Q}\equiv\mathrm{d}Q/\mathrm{d}t Q ˙ ≡ d Q / d t
The system also requires an expression for the evolution of the plasma-inertia measure ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ∂ t ⟨ γ − 3 ⟩ = 0 \partial_t\langle\gamma^{-3}\rangle=0 ∂ t ⟨ γ − 3 ⟩ = 0 ⟨ γ − 3 ⟩ i g n ∼ γ m i n − 3 \langle\gamma^{-3}\rangle_{\mathrm{ign}}\sim\gamma_{\mathrm{min}}^{-3} ⟨ γ − 3 ⟩ ign ∼ γ min − 3
The coupled system has five qualitatively different phases that we discuss next; the last three overlap in time.
Limit Cycle Phases
Loading. We analyze the early loading phase of the discharge cycle with t ≲ t ± t x t \lesssim\sqrt{t_{\pm}t_x} t ≲ t ± t x ε ≈ 1 \varepsilon\approx1 ε ≈ 1 m ± ( t = 0 ) ≈ 1 m_{\pm}(t=0) \approx1 m ± ( t = 0 ) ≈ 1 m x ( t = 0 ) = 0 m_x(t=0)=0 m x ( t = 0 ) = 0 m x m_x m x
m ˙ x = m ± t x (50) \dot{m}_x = \frac{m_{\pm}}{t_x} \tag*{(50)} m ˙ x = t x m ± (50) with a linear-growth solution of
m x ( t ) = t t x . (51) m_x(t) = \frac{t}{t_x}. \tag*{(51)} m x ( t ) = t x t . (51) During this time, the pair density increases only slightly according to
m ˙ ± = 2 m x ( t ) t ± = 2 t t x t ± , (52) \dot{m}_{\pm} = 2\frac{m_x(t)}{t_{\pm}} = \frac{2t}{t_xt_{\pm}}, \tag*{(52)} m ˙ ± = 2 t ± m x ( t ) = t x t ± 2 t , (52) with a solution
m ± ( t ) = 1 + t 2 t ± t x . (53) m_{\pm}(t) = 1 + \frac{t^2}{t_{\pm}t_x}. \tag*{(53)} m ± ( t ) = 1 + t ± t x t 2 . (53) The system remains in the loading phase until m x ( t ) ∼ 1 m_x(t) \sim1 m x ( t ) ∼ 1 t ∼ t x t \sim t_x t ∼ t x
Ignition. The ignition phase starts when the photon multiplicity has become comparable to the plasma multiplicity, m x ( t ) ∼ 1 m_x(t) \sim1 m x ( t ) ∼ 1 ∣ ε ∣ ≥ 1 |\varepsilon| \ge1 ∣ ε ∣ ≥ 1 t ≪ t e s c t \ll t_{\mathrm{esc}} t ≪ t esc
m ˙ ± ≈ 2 m x t ± , (54) \dot{m}_{\pm} \approx2\frac{m_x}{t_{\pm}}, \tag*{(54)} m ˙ ± ≈ 2 t ± m x , (54) m ˙ x ≈ − m x t ± + m ± t x . (55) \dot{m}_x \approx-\frac{m_x}{t_{\pm}}+\frac{m_{\pm}}{t_x}. \tag*{(55)} m ˙ x ≈ − t ± m x + t x m ± . (55) Importantly, for t ≳ t ± t \gtrsim t_{\pm} t ≳ t ±
m ± / x ( t ) ∝ exp ( t t ± ) , (56) m_{\pm/x}(t) \propto\exp\left(\frac{t}{t_{\pm}}\right), \tag*{(56)} m ± / x ( t ) ∝ exp ( t ± t ) , (56) when t ± ∼ t x t_{\pm} \sim t_x t ± ∼ t x m ± ≫ 1 m_{\pm} \gg1 m ± ≫ 1 j ± ∝ m ± j_{\pm} \propto m_{\pm} j ± ∝ m ± ε = 1 → 0 \varepsilon=1 \to0 ε = 1 → 0
Oscillations. The sudden screening of the gap leads to electric field oscillations because the plasma-induced current has inertia and does not instantly respond to the electric field changes. For m ± m_{\pm} m ± m x m_x m x [31 ] , [51 ] , [15 ] , [37 ] .
Combining the equations for ε ˙ \dot{\varepsilon} ε ˙ j ˙ \dot{j} j ˙ m ± ( t ) ∼ m ± m_{\pm}(t) \sim m_{\pm} m ± ( t ) ∼ m ±
ε ¨ + ε ˙ t e s c + ω o s c 2 ε = α t e s c 2 , (57) \ddot{\varepsilon}+\frac{\dot{\varepsilon}}{t_{\mathrm{esc}}}+\omega_{\mathrm{osc}}^2\varepsilon=\frac{\alpha}{t_{\mathrm{esc}}^2}, \tag*{(57)} ε ¨ + t esc ε ˙ + ω osc 2 ε = t esc 2 α , (57) where the oscillation frequency
ω o s c 2 = m ± t a c c t e s c ⟨ 1 γ 3 ⟩ = ω p , c o 2 m ± ⟨ γ − 3 ⟩ . (58) \omega_{\mathrm{osc}}^2=\frac{m_{\pm}}{t_{\mathrm{acc}}t_{\mathrm{esc}}}\left\langle\frac{1}{\gamma^3}\right\rangle=\omega_{\mathrm{p,co}}^2m_{\pm}\left\langle\gamma^{-3}\right\rangle. \tag*{(58)} ω osc 2 = t acc t esc m ± ⟨ γ 3 1 ⟩ = ω p , co 2 m ± ⟨ γ − 3 ⟩ . (58) The discharge rings at the plasma frequency of the gap, suppressed by the effective inertia γ 3 m e \gamma^3m_e γ 3 m e ζ = ( 2 ω o s c t e s c ) − 1 ≪ 1 \zeta=(2\omega_{\mathrm{osc}}t_{\mathrm{esc}})^{-1}\ll1 ζ = ( 2 ω osc t esc ) − 1 ≪ 1 m ± m_{\pm} m ± t d a m p = 2 t e s c t_{\mathrm{damp}}=2t_{\mathrm{esc}} t damp = 2 t esc ∼ ω o s c \sim\omega_{\mathrm{osc}} ∼ ω osc ∼ t e s c \sim t_{\mathrm{esc}} ∼ t esc
During the oscillation phase, ε \varepsilon ε ω o s c − 1 ≪ t e s c \omega_{\mathrm{osc}}^{-1}\ll t_{\mathrm{esc}} ω osc − 1 ≪ t esc
j ( t ) ≈ m ± ε t e s c t a c c ⟨ 1 γ 3 ⟩ . (59) j(t)\approx m_{\pm}\varepsilon\frac{t_{\mathrm{esc}}}{t_{\mathrm{acc}}}\left\langle\frac{1}{\gamma^3}\right\rangle. \tag*{(59)} j ( t ) ≈ m ± ε t acc t esc ⟨ γ 3 1 ⟩ . (59) Therefore, we can compare this expression to the earlier definition of ∣ j ∣ = m ± ∣ β ∥ ∣ |j|=m_{\pm}|\beta_{\parallel}| ∣ j ∣ = m ± ∣ β ∥ ∣
∣ β ∥ ∣ ≈ ∣ ε ∣ t e s c t a c c ⟨ 1 γ 3 ⟩ . (60) |\beta_{\parallel}|\approx|\varepsilon|\frac{t_{\mathrm{esc}}}{t_{\mathrm{acc}}}\left\langle\frac{1}{\gamma^3}\right\rangle. \tag*{(60)} ∣ β ∥ ∣ ≈ ∣ ε ∣ t acc t esc ⟨ γ 3 1 ⟩ . (60) A good proxy for the current is then j ≈ C 1 ε m ± j\approx C_1\varepsilon m_{\pm} j ≈ C 1 ε m ± C 1 C_1 C 1 ∣ β ∥ ∣ < 1 |\beta_{\parallel}|<1 ∣ β ∥ ∣ < 1 ∣ j ∣ → m ± |j|\to m_{\pm} ∣ j ∣ → m ±
Saturation. The sudden screening of the gap leads to saturation of the plasma and photon generation. The screened electric field halts curvature-photon production since H ( ∣ ε ∣ − 1 ) = 0 \mathcal{H}(|\varepsilon|-1)=0 H ( ∣ ε ∣ − 1 ) = 0
m ˙ ± = 2 m x t ± − m ± t e s c , (61) \dot{m}_{\pm}=2\frac{m_x}{t_{\pm}}-\frac{m_{\pm}}{t_{\mathrm{esc}}}, \tag*{(61)} m ˙ ± = 2 t ± m x − t esc m ± , (61) m ˙ x = − m x t ± − m x t e s c , (62) \dot{m}_x=-\frac{m_x}{t_{\pm}}-\frac{m_x}{t_{\mathrm{esc}}}, \tag*{(62)} m ˙ x = − t ± m x − t esc m x , (62) with an exponential-decay solution for the photons,
m x ( t ) = M x ′ exp ( − t t e s c − t t ± ) , (63) m_x(t)=M'_x\exp\left(-\frac{t}{t_{\mathrm{esc}}}-\frac{t}{t_{\pm}}\right), \tag*{(63)} m x ( t ) = M x ′ exp ( − t esc t − t ± t ) , (63) and a round-topped solution for the pairs,
m ± ( t ) = m x ( t ) { exp ( t t ± ) [ M ± ′ M x ′ + 2 ] − 2 } , (64) m_{\pm}(t)=m_x(t)\left\{\exp\left(\frac{t}{t_{\pm}}\right)\left[\frac{M'_{\pm}}{M'_x}+2\right]-2\right\}, \tag*{(64)} m ± ( t ) = m x ( t ) { exp ( t ± t ) [ M x ′ M ± ′ + 2 ] − 2 } , (64) where M x ′ M'_x M x ′ M ± ′ M'_{\pm} M ± ′ m ± ( t ) m_{\pm}(t) m ± ( t ) t d = t ± ln { 2 M x ′ ( t e s c + t ± ) / [ t ± ( M ± ′ + 2 M x ′ ) ] } t_d=t_{\pm}\ln\{2M'_x(t_{\mathrm{esc}}+t_{\pm})/[t_{\pm}(M'_{\pm}+2M'_x)]\} t d = t ± ln { 2 M x ′ ( t esc + t ± ) / [ t ± ( M ± ′ + 2 M x ′ )]} 2 M x ′ / M ± ′ > t ± / t e s c 2M'_x/M'_{\pm}>t_{\pm}/t_{\mathrm{esc}} 2 M x ′ / M ± ′ > t ± / t esc
As a caveat, this phase models the pair cascade only up to the second generation. In contrast, our PIC simulations find that the cascade generates up to six generations of pairs (cf. the five to eight generations reported by [50 ] ). Therefore, a realistic saturation phase can last longer and produce more pairs.
Decay. In the last decay phase, pair production fades as most high-energy photons are already converted to pairs. At late times (t ≫ t ± t \gg t_{\pm} t ≫ t ±
m ± ( t ) ≈ ( M ± ′ + 2 M x ′ ) exp ( − t t e s c ) ≈ M ± exp ( − t t e s c ) , (65) m_{\pm}(t) \approx\left(M'_{\pm} + 2M'_{x}\right) \exp\left(-\frac{t}{t_{\mathrm{esc}}}\right) \approx M_{\pm} \exp\left(-\frac{t}{t_{\mathrm{esc}}}\right), \tag*{(65)} m ± ( t ) ≈ ( M ± ′ + 2 M x ′ ) exp ( − t esc t ) ≈ M ± exp ( − t esc t ) , (65) where we have introduced the peak multiplicity M ± ≈ M ± ′ + 2 M x ′ M_{\pm} \approx M'_{\pm} + 2M'_{x} M ± ≈ M ± ′ + 2 M x ′
As the particle density fades away, the gap starts to reform and the electric field grows as ε ˙ ≈ α / t e s c \dot{\varepsilon} \approx\alpha/t_{\mathrm{esc}} ε ˙ ≈ α / t esc ε ( t ) = α ( t / t e s c ) \varepsilon(t) = \alpha(t/t_{\mathrm{esc}}) ε ( t ) = α ( t / t esc ) ∣ j ∣ ≈ m ± ∣ β ∥ ∣ < ∣ α ∣ |j| \approx m_{\pm}|\beta_{\parallel}| < |\alpha| ∣ j ∣ ≈ m ± ∣ β ∥ ∣ < ∣ α ∣ ∣ β ∥ ∣ → 1 |\beta_{\parallel}| \to1 ∣ β ∥ ∣ → 1 m ± ≲ ∣ α ∣ ∼ m_{\pm} \lesssim|\alpha| \sim m ± ≲ ∣ α ∣ ∼
Macroscopic Properties
The peak multiplicity M ± M_{\pm} M ±
M ± ∼ 10 t a c c t x γ min − 3 . (66) M_{\pm} \sim10 \frac{t_{\mathrm{acc}}}{t_x \gamma_{\min}^{-3}}. \tag*{(66)} M ± ∼ 10 t x γ m i n − 3 t acc . (66) For the typical pulsar parameters, we obtain M ± ∼ 5 × 10 3 M_{\pm} \sim5 \times10^{3} M ± ∼ 5 × 1 0 3
The total cycle duration t c y c l e t_{\mathrm{cycle}} t cycle m ± ( t c y c l e ) = ∣ α ∣ m_{\pm}(t_{\mathrm{cycle}}) = |\alpha| m ± ( t cycle ) = ∣ α ∣ (65) —is dominated by the decay phase and is given as
t c y c l e ≈ ln ( M ± ∣ α ∣ ) t e s c . (67) t_{\mathrm{cycle}} \approx\ln\left(\frac{M_{\pm}}{|\alpha|}\right)t_{\mathrm{esc}}. \tag*{(67)} t cycle ≈ ln ( ∣ α ∣ M ± ) t esc . (67) The time-averaged multiplicity is
⟨ m ± ⟩ = 1 t c y c l e ∫ 0 t c y c l e m ± ( t ) d t ≈ M ± ln M ± , (68) \langle m_{\pm}\rangle= \frac{1}{t_{\mathrm{cycle}}}\int_{0}^{t_{\mathrm{cycle}}}m_{\pm}(t)\,\mathrm{d}t \approx\frac{M_{\pm}}{\ln M_{\pm}}, \tag*{(68)} ⟨ m ± ⟩ = t cycle 1 ∫ 0 t cycle m ± ( t ) d t ≈ ln M ± M ± , (68) where we have used (65) for m ± ( t ) m_{\pm}(t) m ± ( t )
Numerical Solution
The full numerical solution of the ODE system is shown in Figure 2 . We use the Radau method (implicit Runge–Kutta method of the Radau IIA family of order 5) to solve the stiff ODE system. The solution is numerically straightforward to obtain and is not sensitive to the initial conditions.
Figure 2. Numerical solution of the analytic discharge model for the fiducial parameters α = 2 \alpha=2 α = 2 t x = t ± = 10 − 3 t e s c t_x=t_{\pm}=10^{-3}t_{\mathrm{esc}} t x = t ± = 1 0 − 3 t esc t a c c = 4 × 10 − 8 t e s c t_{\mathrm{acc}}=4\times10^{-8}t_{\mathrm{esc}} t acc = 4 × 1 0 − 8 t esc ⟨ γ − 3 ⟩ = 8 × 10 − 8 \langle\gamma^{-3}\rangle=8\times10^{-8} ⟨ γ − 3 ⟩ = 8 × 1 0 − 8 α \alpha α ⟨ β ⟩ = ± j / m ± \langle\beta\rangle=\pm j/m_{\pm} ⟨ β ⟩ = ± j / m ± ε ( t e s c / t a c c ) ⟨ γ − 3 ⟩ \varepsilon(t_{\mathrm{esc}}/t_{\mathrm{acc}})\langle\gamma^{-3}\rangle ε ( t esc / t acc ) ⟨ γ − 3 ⟩ ∣ β ∣ = 1 |\beta|=1 ∣ β ∣ = 1
The numerical result exhibits all the phases described in the Sect. 5.3 and has limit-cycle behavior. The measured peak multiplicity, M ± ≈ 5 × 10 3 M_{\pm} \approx5 \times10^{3} M ± ≈ 5 × 1 0 3 ∼ 20 % \sim20\% ∼ 20% (66) . The measured cycle length ( 9.1 ± 0.1 ) t e s c (9.1 \pm0.1)t_{\mathrm{esc}} ( 9.1 ± 0.1 ) t esc t c y c l e ≈ ln ( M ± / ∣ α ∣ ) t e s c ≈ 7.8 t e s c t_{\mathrm{cycle}} \approx\ln(M_{\pm}/|\alpha|)t_{\mathrm{esc}} \approx7.8\,t_{\mathrm{esc}} t cycle ≈ ln ( M ± /∣ α ∣ ) t esc ≈ 7.8 t esc t ± = t x t_{\pm}=t_x t ± = t x t ± ≈ 2 × 10 − 2 t e s c t_{\pm} \approx2 \times10^{-2}t_{\mathrm{esc}} t ± ≈ 2 × 1 0 − 2 t esc t x = 10 − 3 t e s c t_x=10^{-3}t_{\mathrm{esc}} t x = 1 0 − 3 t esc
Kinetic Plasma Simulations
We compare the theoretical analysis to numerical simulations. We study the radiative plasma dynamics of the magnetospheric circuit with fully kinetic PIC simulations using the RUNKO framework ([33 ] ). The simulations use realistic system parameters together with exact QED cross sections.
Numerical Setup
We study the circuit with one-dimensional electrodynamic simulations that evolve the magnetic and electric fields, B ( h ) \mathbf{B}(h) B ( h ) E ( h ) \mathbf{E}(h) E ( h ) h h h L b o x ∼ 3 × 10 5 Δ x L_{\rm box} \sim3 \times10^{5}\Delta x L box ∼ 3 × 1 0 5 Δ x Δ x \Delta x Δ x H g a p = 0.8 L b o x H_{\rm gap} = 0.8L_{\rm box} H gap = 0.8 L box c / ω p , c o ≈ 33 Δ x c/\omega_{\rm p,co} \approx33\Delta x c / ω p , co ≈ 33Δ x H g a p / ( c / ω p , c o ) ≈ 7400 H_{\rm gap}/(c/\omega_{\rm p,co}) \approx7400 H gap / ( c / ω p , co ) ≈ 7400 ≈ 5000 \approx5000 ≈ 5000 c ^ = c Δ t / Δ x = 0.45 \hat{c} = c\Delta t/\Delta x = 0.45 c ^ = c Δ t /Δ x = 0.45 ∼ 20 t e s c \sim20t_{\rm esc} ∼ 20 t esc 10 7 10^{7} 1 0 7 m p / m e = 1836 m_{\rm p}/m_{\rm e} = 1836 m p / m e = 1836
We take the initial charge density η = 0 \eta= 0 η = 0 h < H g a p h < H_{\rm gap} h < H gap η = η c o \eta= \eta_{\rm co} η = η co h ≥ H g a p h \ge H_{\rm gap} h ≥ H gap h < H g a p h < H_{\rm gap} h < H gap ∇ ⋅ E = 4 π ( η − η c o ) \nabla\cdot\mathbf{E} = 4\pi(\eta- \eta_{\rm co}) ∇ ⋅ E = 4 π ( η − η co ) E ( h ) E(h) E ( h ) E h ( 0 ) = 0 E_{h}(0) = 0 E h ( 0 ) = 0 E h ( H g a p ) = E r o t E_{h}(H_{\rm gap}) = E_{\rm rot} E h ( H gap ) = E rot E h ≡ E ∥ E_{h} \equiv E_{\parallel} E h ≡ E ∥ h ^ \hat{\mathbf{h}} h ^ B = B ⋆ h ^ \mathbf{B} = B_{\star}\hat{\mathbf{h}} B = B ⋆ h ^
We advance the electromagnetic fields with a second-order field solver. A first-order method interpolates particle positions and a ZigZag algorithm deposits the current. We push the particles with a Boris particle pusher. The current is smoothed with a digital filter (with a binomial 3-point shape) with 8 passes per step. In addition to the standard PIC algorithm, we add a virtual current to model the system in a co-rotating frame with β r o t = 3 × 10 − 6 \beta_{\rm rot} = 3 \times10^{-6} β rot = 3 × 1 0 − 6 j m = 2 j c o j_{\rm m} = 2j_{\rm co} j m = 2 j co α = 2 \alpha= 2 α = 2 α = − 2 \alpha= -2 α = − 2
The left simulation boundary includes the star and the atmosphere. The atmosphere supplies particles to the magnetosphere such that the resulting current will always screen the electric field for h < 0 h < 0 h < 0 k B T / m e c 2 = 0.3 k_{\rm B}T/m_{\rm e}c^{2} = 0.3 k B T / m e c 2 = 0.3
The QED module and the adaptive Monte Carlo routines are discussed in [34 ] (Supplementary Material). Here, we take into account synchrotron radiation and one-photon (nonlinear Breit–Wheeler) pair creation. We re-scale the QED processes by a factor N m p = 2 × 10 11 \mathcal{N}_{\rm mp} = 2 \times10^{11} N mp = 2 × 1 0 11 N m p \mathcal{N}_{\rm mp} N mp r e / Δ x = ( 4 π N m p c ^ 2 ) − 1 r_{\rm e}/\Delta x = (4\pi\mathcal{N}_{\rm mp}\hat{c}^{2})^{-1} r e /Δ x = ( 4 π N mp c ^ 2 ) − 1 r e / H g a p ∼ 10 − 17 r_{\rm e}/H_{\rm gap} \sim10^{-17} r e / H gap ∼ 1 0 − 17 Δ t > t r \Delta t > t_{r} Δ t > t r B y B_{y} B y B y ∝ sin [ ( h − h 0 ) / R c u r v ] B_{y} \propto\sin[(h-h_{0})/R_{\rm curv}] B y ∝ sin [( h − h 0 ) / R curv ] h 0 h_{0} h 0 h > H g a p h > H_{\rm gap} h > H gap b = 0.02 b = 0.02 b = 0.02 R c u r v = 30 H g a p R_{\rm curv} = 30H_{\rm gap} R curv = 30 H gap γ r a d / γ g a p ≈ 0.3 \gamma_{\rm rad}/\gamma_{\rm gap} \approx0.3 γ rad / γ gap ≈ 0.3 λ m f p / H g a p ≈ 0.01 \lambda_{\rm mfp}/H_{\rm gap} \approx0.01 λ mfp / H gap ≈ 0.01 γ g a p ≈ 3 × 10 7 \gamma_{\rm gap} \approx3 \times10^{7} γ gap ≈ 3 × 1 0 7 t a c c / t e s c ≈ 2 × 10 − 8 t_{\rm acc}/t_{\rm esc} \approx2 \times10^{-8} t acc / t esc ≈ 2 × 1 0 − 8 t ± / t e s c = λ m f p / H g a p ≈ 10 − 2 t_{\pm}/t_{\rm esc} = \lambda_{\rm mfp}/H_{\rm gap} \approx10^{-2} t ± / t esc = λ mfp / H gap ≈ 1 0 − 2
The simulation setup differs only slightly from previous simulations (notably [48 ] ; [49 ] ) by having realistic system parameters, employing exact differential cross sections for the QED processes, including a thin, gravitationally bound atmosphere, and solving the EM fields in the co-rotating frame. These improvements result in a more realistic particle distribution during the discharge, which, in turn, give a more realistic inverse-plasma-inertia measure ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩
Observed Dynamics
The particle dynamics in the gap follow the picture described in Sect. 4. The evolution is visualized in Figure 3 . In the beginning, we observe a rapid acceleration of the particles to γ r a d \gamma_{\rm rad} γ rad ∼ λ m f p \sim\lambda_{\rm mfp} ∼ λ mfp
Figure 3. Simulation snapshots of the discharge cycle at t / t e s c = 0 t/t_{\mathrm{esc}} = 0 t / t esc = 0 1st column ), 0.1 0.1 0.1 2nd column ), 0.3 0.3 0.3 3rd column ), and 0.4 0.4 0.4 4th column ) for α = 2 \alpha= 2 α = 2 Rows 1–4: Particle phase space, d n ( p ) / d log p \mathrm{d}n(p)/\mathrm{d}\log p d n ( p ) / d log p p = γ β p = \gamma\beta p = γ β 1st row ), positrons (2nd ), and protons (3rd ), and the photon energy spectrum, x d n x ( x ) / d log x x\,\mathrm{d}n_x(x)/\mathrm{d}\log x x d n x ( x ) / d log x 4th ). The values range from low to high as blue to green to red. Vertical thin lines show the γ g a p \gamma_{\mathrm{gap}} γ gap γ r a d \gamma_{\mathrm{rad}} γ rad x c u r v x_{\mathrm{curv}} x curv p = 0 p = 0 p = 0 5th row: Multiplicity of the electrons (blue curve), positrons (orange), and high-energy photons with x > 2 x > 2 x > 2 6th row: Mean Lorentz factor of electrons (blue) and positrons (orange). 7th row: Inverse plasma-inertia measure, ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ 8th row: Gap electric field E ∥ E_{\parallel} E ∥ E r o t E_{\mathrm{rot}} E rot
The subsequent gap dynamics follow the phases outlined in Sect. 5. The empty gap is loaded with photons. The gap stays in the loading and ignition phases until the plasma cloud develops a moderate bulk velocity and an electric current forms. The current then screens the electric field. The abrupt screening leads to overshooting the equilibrium value and the electric field in the gap is left ringing with the frequency ω o s c \omega_{\rm osc} ω osc Figure 4 visualizes the E ∥ 2 ( h , t ) E_{\parallel}^{2}(h,t) E ∥ 2 ( h , t ) t ≈ 0.2 t e s c t \approx0.2t_{\rm esc} t ≈ 0.2 t esc h ∈ [ 0.5 , 1 ] H g a p h \in[0.5,1]H_{\rm gap} h ∈ [ 0.5 , 1 ] H gap ∼ t e s c \sim t_{\mathrm{esc}} ∼ t esc (67) because the effective gap length is closer to ∼ 0.2 H g a p \sim0.2H_{\mathrm{gap}} ∼ 0.2 H gap
Figure 4. Space-time visualization of the squared electric field E ∥ 2 / E r o t 2 E_{\parallel}^{2}/E_{\rm rot}^{2} E ∥ 2 / E rot 2 t > 0.2 t e s c t > 0.2t_{\rm esc} t > 0.2 t esc
The global gap dynamics are visualized in Figure 5 for the first discharge-episode and the results can be compared to the numerical solutions of the concurrency model (see Figure 2 ). The results are in reasonable agreement. The peak multiplicity evolves broadly as presented in Sect. 5: it grows exponentially, saturates, reaches a maximum, and, finally, decays exponentially. The main difference is the prolonged plateau-like saturation phase, caused by the inclusion of particle populations beyond the second generation. Another visual difference is in the decay phase: the pair cloud is moving and slides out of the finite simulation domain at around t ≈ 1.1 t e s c t \approx1.1t_{\mathrm{esc}} t ≈ 1.1 t esc
Figure 5. Measurements of the discharge region quantities from a full kinetic simulation during the first pair cascade event. Top panel shows the local peak multiplicity M ± M_{\pm} M ± ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ⟨ γ − 3 ⟩ ∼ 10 − 7 \langle\gamma^{-3}\rangle\sim10^{-7} ⟨ γ − 3 ⟩ ∼ 1 0 − 7 ∼ 10 − 3 \sim10^{-3} ∼ 1 0 − 3 ∼ 0.1 \sim0.1 ∼ 0.1
The inverse plasma inertia ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ⟨ γ − 3 ⟩ i n i t ≲ 10 − 10 \langle\gamma^{-3}\rangle_{\mathrm{init}} \lesssim10^{-10} ⟨ γ − 3 ⟩ init ≲ 1 0 − 10 m ± ∼ 1 m_{\pm} \sim1 m ± ∼ 1 ⟨ γ − 3 ⟩ i g n ∼ 10 − 7 \langle\gamma^{-3}\rangle_{\mathrm{ign}} \sim10^{-7} ⟨ γ − 3 ⟩ ign ∼ 1 0 − 7 t ∼ t e s c t \sim t_{\mathrm{esc}} t ∼ t esc ⟨ γ ⟩ ≈ γ min ∼ 400 \langle\gamma\rangle\approx\gamma_{\min} \sim400 ⟨ γ ⟩ ≈ γ m i n ∼ 400 ⟨ γ − 3 ⟩ s a t ∼ 1 2 ⟨ γ ⟩ − 1 ∼ 10 − 3 \langle\gamma^{-3}\rangle_{\mathrm{sat}} \sim\frac{1}{2}\langle\gamma\rangle^{-1} \sim10^{-3} ⟨ γ − 3 ⟩ sat ∼ 2 1 ⟨ γ ⟩ − 1 ∼ 1 0 − 3 t ∼ 2 t e s c t \sim2t_{\mathrm{esc}} t ∼ 2 t esc ⟨ γ ⟩ ∼ 10 \langle\gamma\rangle\sim10 ⟨ γ ⟩ ∼ 10 ⟨ γ − 3 ⟩ l a t e ∼ 0.1 \langle\gamma^{-3}\rangle_{\mathrm{late}} \sim0.1 ⟨ γ − 3 ⟩ late ∼ 0.1 ⟨ γ − 3 ⟩ l a t e \langle\gamma^{-3}\rangle_{\mathrm{late}} ⟨ γ − 3 ⟩ late
The general features of the discharges are very similar for the α < 0 \alpha< 0 α < 0 [48 ] ; [49 ] ). In particular, we have tested α = − 2 \alpha=-2 α = − 2
Peak Multiplicity
The measured peak multiplicity is M ± ∼ 10 3 − 10 4 M_{\pm} \sim10^{3}-10^{4} M ± ∼ 1 0 3 − 1 0 4 ⟨ M ± ⟩ ∼ \langle M_{\pm}\rangle\sim ⟨ M ± ⟩ ∼ × 10 2 \times10^{2} × 1 0 2 ⟨ γ − 3 ⟩ i g n ∼ 10 − 7 \langle\gamma^{-3} \rangle_{\mathrm{ign}} \sim10^{-7} ⟨ γ − 3 ⟩ ign ∼ 1 0 − 7 (66) gives an estimate for the peak multiplicity of M ± ∼ 5 × 10 3 M_{\pm} \sim5 \times10^{3} M ± ∼ 5 × 1 0 3
The measured pair yield varies depending on the simulation input parameters and modeling choices, such as γ g a p \gamma_{\mathrm{gap}} γ gap γ r a d \gamma_{\mathrm{rad}} γ rad E E E
Particle Distributions
Particle spectra during the discharge are visualized in Figure 6 . The particles are spatially averaged in a region h ∈ [ 0 , 1 ] H g a p h \in[0,1]H_{\mathrm{gap}} h ∈ [ 0 , 1 ] H gap t ≲ 0.2 t e s c t \lesssim0.2t_{\mathrm{esc}} t ≲ 0.2 t esc γ 1 ∼ 10 7 \gamma_1 \sim10^{7} γ 1 ∼ 1 0 7 ≈ γ r a d \approx\gamma_{\mathrm{rad}} ≈ γ rad γ r a d \gamma_{\mathrm{rad}} γ rad E h E_h E h E r o t E_{\mathrm{rot}} E rot γ 2 ∼ 10 5 \gamma_2 \sim10^{5} γ 2 ∼ 1 0 5 ∼ 10 − 2 γ r a d \sim10^{-2}\gamma_{\mathrm{rad}} ∼ 1 0 − 2 γ rad t ∼ [ 0.5 − 1 ] t e s c t \sim[0.5-1]t_{\mathrm{esc}} t ∼ [ 0.5 − 1 ] t esc
Figure 6. Average particle spectra in the gap as a function of time. Top panel shows the average momentum distribution for the pair plasma across the full gap (i.e., h ∈ [ 0 , H g a p ] h \in[0,H_{\mathrm{gap}}] h ∈ [ 0 , H gap ]
The photon spectra have similar peaks at x 1 ∼ 10 − 2 γ 1 x_1 \sim10^{-2}\gamma_1 x 1 ∼ 1 0 − 2 γ 1 x 2 ∼ 10 − 1 γ 2 x_2 \sim10^{-1}\gamma_2 x 2 ∼ 1 0 − 1 γ 2 t ≲ 0.5 t e s c t \lesssim0.5t_{\mathrm{esc}} t ≲ 0.5 t esc x ∼ 10 2 x \sim10^{2} x ∼ 1 0 2 10 5 10^{5} 1 0 5 x ≲ 10 x \lesssim10 x ≲ 10 t ≳ t e s c t \gtrsim t_{\mathrm{esc}} t ≳ t esc
The particle distributions are not exactly the same for electrons vs. positrons and for inward- vs. outward-moving populations; instead they evolve as a function of time. At early times (t ≲ 0.2 t e s c t \lesssim0.2t_{\mathrm{esc}} t ≲ 0.2 t esc E E E t ∼ [ 0.3 − 1 ] t e s c t \sim[0.3-1]t_{\mathrm{esc}} t ∼ [ 0.3 − 1 ] t esc Figure 3 ). For late times (t ≳ t e s c t \gtrsim t_{\mathrm{esc}} t ≳ t esc
Return-Particle Flux
The characteristics of the return-particle flux are important because the in-falling particles cause heating of the atmosphere [41 ] . The particle distribution for the in-falling pairs is visualized in Figure 7 . The distribution closely resembles the overall particle spectra in the gap (see Figure 6 ) because the majority of the particles are produced at the late saturation stage of the burst, during which the particle populations have had time to mix.
Figure 7. Return-particle spectra as a function of time (colored curves). The time-averaged spectrum of particles bombarding the atmosphere is shown with a thick black curve. The vertical dashed line corresponds to γ = 400 \gamma= 400 γ = 400 γ − 0.3 \gamma^{-0.3} γ − 0.3
A thermal distribution with ⟨ γ ⟩ ∼ 400 \langle\gamma\rangle\sim400 ⟨ γ ⟩ ∼ 400 ∼ 0.1 γ r a d \sim0.1\gamma_{\mathrm{rad}} ∼ 0.1 γ rad d n ± / d ln γ ∝ γ − s d n_{\pm}/d\ln\gamma\propto\gamma^{-s} d n ± / d ln γ ∝ γ − s s ≈ 0.3 s \approx0.3 s ≈ 0.3
Electromagnetic Oscillations
The discharge excites strong electromagnetic oscillations in the gap. These oscillations excite (super-luminal-branch) O-mode waves in 2D/3D when the discharge cloud front is tilted w.r.t. the local magnetic field, i.e., ∂ y n ± ≠ 0 \partial_y n_{\pm} \ne0 ∂ y n ± = 0 B = B y ^ \mathbf{B} = B\hat{\mathbf{y}} B = B y ^ [39 ] . Therefore, it is interesting to analyze the properties of these oscillations because they can be connected to the pulsar radio emission.
The oscillation phase of the discharge cycle excites primarily E h E_h E h δ E h / E r o t ≡ η E ∼ 10 − 2 \delta E_h/E_{\mathrm{rot}} \equiv\eta_E \sim10^{-2} δ E h / E rot ≡ η E ∼ 1 0 − 2 Figure 8 ). The wavenumber (k k k ∣ E ~ h ∣ 2 ∝ k − 2 |\tilde{E}_h|^2 \propto k^{-2} ∣ E ~ h ∣ 2 ∝ k − 2 S ν ∝ ν − 1.4 S_{\nu} \propto\nu^{-1.4} S ν ∝ ν − 1.4 ν \nu ν [4 ] {cite:p}`timokhin_2013,kumar_2025``. The measured power spectrum has additional time-variable features (dips and spikes) at values of k k k k o s c ≈ ω o s c / c k_{\mathrm{osc}} \approx\omega_{\mathrm{osc}}/c k osc ≈ ω osc / c
Figure 8. Top panel shows the electric field oscillations in the gap for E h E_h E h E y E_y E y E z E_z E z ∣ E ~ h ∣ 2 |\tilde{E}_h|^2 ∣ E ~ h ∣ 2 ∣ E ~ y ∣ 2 |\tilde{E}_y|^2 ∣ E ~ y ∣ 2 ∣ E ~ z ∣ 2 |\tilde{E}_z|^2 ∣ E ~ z ∣ 2 y y y z z z × 10 6 \times10^6 × 1 0 6 ∣ E ~ h ∣ 2 ∝ k − 2 |\tilde{E}_h|^2 \propto k^{-2} ∣ E ~ h ∣ 2 ∝ k − 2 ∣ E ~ y ∣ 2 ∝ k 0 |\tilde{E}_y|^2 \propto k^0 ∣ E ~ y ∣ 2 ∝ k 0 ∣ E ~ z ∣ 2 ∝ k 0 |\tilde{E}_z|^2 \propto k^0 ∣ E ~ z ∣ 2 ∝ k 0 k o s c = ω o s c / c k_{\mathrm{osc}} = \omega_{\mathrm{osc}}/c k osc = ω osc / c M ± = 5 × 10 2 M_{\pm}=5 \times10^2 M ± = 5 × 1 0 2 ⟨ γ − 3 ⟩ = 10 − 3 \langle\gamma^{-3} \rangle= 10^{-3} ⟨ γ − 3 ⟩ = 1 0 − 3 t = 0.5 t e s c t=0.5t_{\mathrm{esc}} t = 0.5 t esc
As the discharge evolves, the frequency changes. At ignition, the plasma density has only started to grow from the co-rotation value and the plasma is dominated by the freshly accelerated beam, ⟨ γ − 3 ⟩ i g n ∼ 10 − 7 \langle\gamma^{-3} \rangle_{\mathrm{ign}} \sim10^{-7} ⟨ γ − 3 ⟩ ign ∼ 1 0 − 7 m ± ∼ 10 m_{\pm} \sim10 m ± ∼ 10 ω o s c t e s c ∼ ( m ± t e s c t a c c ⟨ γ − 3 ⟩ i g n ) 1 / 2 ∼ 5 \omega_{\mathrm{osc}}t_{\mathrm{esc}} \sim\left(m_{\pm}\frac{t_{\mathrm{esc}}}{t_{\mathrm{acc}}}\langle\gamma^{-3} \rangle_{\mathrm{ign}}\right)^{1/2} \sim5 ω osc t esc ∼ ( m ± t acc t esc ⟨ γ − 3 ⟩ ign ) 1/2 ∼ 5 ω o s c / 2 π ∼ 1 M H z \omega_{\mathrm{osc}}/2\pi\sim1\ \mathrm{MHz} ω osc /2 π ∼ 1 MHz M ± ∼ 10 4 M_{\pm} \sim10^{4} M ± ∼ 1 0 4 ⟨ γ − 3 ⟩ s a t ∼ 10 − 3 \langle\gamma^{-3} \rangle_{\mathrm{sat}} \sim10^{-3} ⟨ γ − 3 ⟩ sat ∼ 1 0 − 3 ω o s c t e s c ∼ 10 4 \omega_{\mathrm{osc}}t_{\mathrm{esc}} \sim10^{4} ω osc t esc ∼ 1 0 4 ω o s c / 2 π ∼ 3 G H z \omega_{\mathrm{osc}}/2\pi\sim3\ \mathrm{GHz} ω osc /2 π ∼ 3 GHz m ± → 1 m_{\pm} \to1 m ± → 1 ⟨ γ − 3 ⟩ l a t e ∼ 0.1 \langle\gamma^{-3} \rangle_{\mathrm{late}} \sim0.1 ⟨ γ − 3 ⟩ late ∼ 0.1 ω o s c t e s c ∼ 10 3 \omega_{\mathrm{osc}}t_{\mathrm{esc}} \sim10^{3} ω osc t esc ∼ 1 0 3 ω o s c / 2 π ∼ 400 M H z \omega_{\mathrm{osc}}/2\pi\sim400\ \mathrm{MHz} ω osc /2 π ∼ 400 MHz
The excited oscillations are strong. We can quantify their strength with the nonlinearity parameter
a ≡ e δ E h m e c ω o s c = η E ( 2 γ g a p m ± ⟨ γ − 3 ⟩ ) 1 / 2 (69) a \equiv\frac{e\delta E_h}{m_e c\omega_{\mathrm{osc}}} = \eta_E\left(\frac{2\gamma_{\mathrm{gap}}}{m_{\pm}\langle\gamma^{-3}\rangle}\right)^{1/2} \tag*{(69)} a ≡ m e c ω osc eδ E h = η E ( m ± ⟨ γ − 3 ⟩ 2 γ gap ) 1/2 (69) ∼ 10 η E , − 2 B 12 1 / 2 P 0 1 / 2 ( m ± 10 4 ) − 1 / 2 ( ⟨ γ − 3 ⟩ 10 − 3 ) − 1 / 2 , \sim10\eta_{E,-2}\frac{B_{12}^{1/2}}{P_0^{1/2}}\left(\frac{m_{\pm}}{10^4}\right)^{-1/2}\left(\frac{\langle\gamma^{-3}\rangle}{10^{-3}}\right)^{-1/2}, ∼ 10 η E , − 2 P 0 1/2 B 12 1/2 ( 1 0 4 m ± ) − 1/2 ( 1 0 − 3 ⟨ γ − 3 ⟩ ) − 1/2 , expressed using the saturation-phase values. The parameter is a ∼ 10 4 η E , − 2 \sim10^{4}\eta_{E,-2} ∼ 1 0 4 η E , − 2 ∼ 100 η E , − 2 \sim100\eta_{E,-2} ∼ 100 η E , − 2 L r a d i o ∼ 10 28 e r g s − 1 L_{\mathrm{radio}} \sim10^{28}\ \mathrm{erg\ s^{-1}} L radio ∼ 1 0 28 erg s − 1
a ∼ 4 π e L r a d i o / c m e c ν r a d i o R p c ∼ 300 L 28 1 / 2 P 0 1 / 2 ( ν r a d i o 1 G H z ) − 1 . (70) a \sim\frac{4\pi e\sqrt{L_{\mathrm{radio}}/c}}{m_e c\nu_{\mathrm{radio}}R_{\mathrm{pc}}} \sim300L_{28}^{1/2}P_0^{1/2}\left(\frac{\nu_{\mathrm{radio}}}{1\ \mathrm{GHz}}\right)^{-1}. \tag*{(70)} a ∼ m e c ν radio R pc 4 π e L radio / c ∼ 300 L 28 1/2 P 0 1/2 ( 1 GHz ν radio ) − 1 . (70) The luminosity-implied relative amplitude, δ E / E r o t ∼ 4 × 10 − 2 L 28 1 / 2 P 0 1 / 2 B 12 − 1 \delta E/E_{\mathrm{rot}} \sim4 \times10^{-2}L_{28}^{1/2}P_0^{1/2}B_{12}^{-1} δ E / E rot ∼ 4 × 1 0 − 2 L 28 1/2 P 0 1/2 B 12 − 1 η E ∼ 10 − 2 \eta_E \sim10^{-2} η E ∼ 1 0 − 2 a ≫ 1 a \gg1 a ≫ 1 p ∼ a m e c p \sim a m_e c p ∼ a m e c
The rotating frame (see Appendix A) also enables excitation of the E y E_y E y δ E y ∼ 10 − 6 δ E h \delta E_y \sim10^{-6}\delta E_h δ E y ∼ 1 0 − 6 δ E h ∣ E ^ y ∣ 2 ∝ k 0 |\hat{E}_y|^2 \propto k_0 ∣ E ^ y ∣ 2 ∝ k 0
Summary
Radio-pulsar polar cap discharges generate the pair plasma that powers coherent radio emission. We studied these discharges with analytic estimates and numerical simulations to understand the dynamics of the magnetospheric circuit.
Sects. 2–4 reviewed the plasma physics of polar cap dynamics and re-expressed the historical formulas in unified modern notation. A new concurrency model (Sect. 5) simultaneously describes the time-evolution of the plasma multiplicity, photon density, and electric field. The analytical model encapsulates all key features of pair discharge dynamics: limit-cycle behavior, intermittent pair-plasma generation, and excitation of electric field oscillations. These cycles connect to the microstructures seen in pulsar radio data.
One-dimensional radiative PIC simulations (Sect. 6) validated the analytic theory. The simulations employ exact QED cross sections and realistic plasma parameters, enabling direct analysis of the generated particle distributions. The numerical setup builds on that of [48 ] ; [49 ] —adding realistic system parameters, exact QED cross sections, the atmosphere, and the co-rotating frame—and our results agree well with their work. We characterized the particle spectra of the discharge clouds and the return-particle showers falling back onto the polar caps. The electric field oscillations produced by the discharge have more realistic properties than in previous simulations because they are induced by realistic particle distributions (i.e., realistic ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩
This work focused on synchrotron-based cascades, although other processes can also contribute to the discharge physics. Compton-based cascades, relevant for millisecond pulsars, will be reported in [41 ] (in prep.).
ACKNOWLEDGMENTS
J.N. and T.S. acknowledge useful discussions with A. Timokhin. The authors used Anthropic Claude (Opus-5 and Fable-5) models for assistance with manuscript editing and equation checking.
This work is supported by an ERC grant (ILLUMINATOR, 101114623) and by the Research Council of Finland Centre of Excellence in Neutron-Star Physics (project 374063). TS acknowledges funding by the Research Council of Finland grant No. 368807. The views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. The authors thank the Finnish Computing Competence Infrastructure (FCCI) for supporting this project with computational and data storage resources and for maintaining the HILE cluster.
Software. We use RUNKO ([33 ] ) for the kinetic simulations. Data analysis and visualization software used in this study: NUMPY ([22 ] ), MATPLOTLIB ([24 ] ), PYBIND11 ([25 ] ), and SCIPY ([52 ] ).
Appendix
Electrodynamics in a Moving Frame
The electrodynamics in a moving frame was first discussed by [42 ] [18 ] and [48 ] applied these equations to pulsar magnetospheres. Here, we present their main equations and perform a scaling analysis of the terms.
Consider a rotating frame with an angular velocity Ω \boldsymbol{\Omega} Ω r = R r = R r = R β = Ω × R / c \boldsymbol{\beta} = \boldsymbol{\Omega} \times\boldsymbol{R}/c β = Ω × R / c t t t ∇ \nabla ∇ E E E B B B j j j η \eta η
Maxwell’s equations in the rotating frame at r = R r = R r = R
simplification:
∇ ⋅ E = 4 π ( η − η v i r ) (71) \nabla\cdot E = 4\pi(\eta- \eta_{\mathrm{vir}}) \tag*{(71)} ∇ ⋅ E = 4 π ( η − η vir ) (71) ∇ ⋅ B = 0 (72) \nabla\cdot B = 0 \tag*{(72)} ∇ ⋅ B = 0 (72) ∂ B ∂ t = − c ∇ × E (73) \frac{\partial B}{\partial t} = -c\nabla\times E \tag*{(73)} ∂ t ∂ B = − c ∇ × E (73) ∂ E ∂ t = c ∇ × B − 4 π ( j − j v i r ) . (74) \frac{\partial E}{\partial t} = c\nabla\times B - 4\pi(j - j_{\mathrm{vir}}). \tag*{(74)} ∂ t ∂ E = c ∇ × B − 4 π ( j − j vir ) . (74) The equations appear similar to the standard Maxwell equations except for a virtual charge density
η v i r = − Ω ⋅ B 2 π c + r 4 π R 3 β ⋅ ∇ × B , (75) \eta_{\mathrm{vir}} = -\frac{\boldsymbol{\Omega}\cdot B}{2\pi c} + \frac{r}{4\pi R^3}\boldsymbol{\beta}\cdot\nabla\times B, \tag*{(75)} η vir = − 2 π c Ω ⋅ B + 4 π R 3 r β ⋅ ∇ × B , (75) and current
j v i r = − c r 4 π R { ∇ × ( β × E ) + β × ( ∇ × E ) − ∇ × [ β × ( β × B ) ] } . (76) j_{\mathrm{vir}} = -\frac{cr}{4\pi R}\left\{\nabla\times(\boldsymbol{\beta}\times E) + \boldsymbol{\beta}\times(\nabla\times E) - \nabla\times\left[\boldsymbol{\beta}\times(\boldsymbol{\beta}\times B)\right]\right\}. \tag*{(76)} j vir = − 4 π R cr { ∇ × ( β × E ) + β × ( ∇ × E ) − ∇ × [ β × ( β × B ) ] } . (76) Here the first term in Eq. (75) is the co-rotation charge density η c o = − Ω ⋅ B / ( 2 π c ) \eta_{\rm co}=-\boldsymbol{\Omega}\cdot\mathbf{B}/(2\pi c) η co = − Ω ⋅ B / ( 2 π c ) j c o ≡ η c o c j_{\rm co}\equiv\eta_{\rm co}c j co ≡ η co c
Next, we consider a one-dimensional system along h ^ \hat{h} h ^ B = B ( h ) \mathbf{B}=\mathbf{B}(h) B = B ( h ) E = E ( h ) \mathbf{E}=\mathbf{E}(h) E = E ( h ) h ^ \hat{h} h ^ B = B ∥ h ^ + B ⊥ \mathbf{B}=B_{\parallel}\hat{h}+\mathbf{B}_{\perp} B = B ∥ h ^ + B ⊥ E = E ∥ h ^ + E ⊥ \mathbf{E}=E_{\parallel}\hat{h}+\mathbf{E}_{\perp} E = E ∥ h ^ + E ⊥ ∣ ∇ ∣ ∼ l ∥ − 1 |\nabla|\sim l_{\parallel}^{-1} ∣∇∣ ∼ l ∥ − 1 l ∥ l_{\parallel} l ∥ h ^ \hat{h} h ^ r ∼ R r\sim R r ∼ R
∣ η v i r ∣ ∣ η c o ∣ ∼ O ( 1 ) + O ( B ⊥ B ∥ R l ∥ ) (77) \frac{|\eta_{\rm vir}|}{|\eta_{\rm co}|}\sim\mathcal{O}(1)+\mathcal{O}\left(\frac{B_{\perp}}{B_{\parallel}}\frac{R}{l_{\parallel}}\right) \tag*{(77)} ∣ η co ∣ ∣ η vir ∣ ∼ O ( 1 ) + O ( B ∥ B ⊥ l ∥ R ) (77) and
∣ j v i r ∣ ∣ j c o ∣ ∼ O ( E ∥ B ∥ R l ∥ ) + O ( E ⊥ B ∥ R l ∥ ) + O ( β B ⊥ B ∥ R l ∥ ) . (78) \frac{|j_{\rm vir}|}{|j_{\rm co}|}\sim\mathcal{O}\left(\frac{E_{\parallel}}{B_{\parallel}}\frac{R}{l_{\parallel}}\right)+\mathcal{O}\left(\frac{E_{\perp}}{B_{\parallel}}\frac{R}{l_{\parallel}}\right)+\mathcal{O}\left(\beta\frac{B_{\perp}}{B_{\parallel}}\frac{R}{l_{\parallel}}\right). \tag*{(78)} ∣ j co ∣ ∣ j vir ∣ ∼ O ( B ∥ E ∥ l ∥ R ) + O ( B ∥ E ⊥ l ∥ R ) + O ( β B ∥ B ⊥ l ∥ R ) . (78) Assuming B ⊥ / B ∥ ≪ l ∥ / R B_{\perp}/B_{\parallel}\ll l_{\parallel}/R B ⊥ / B ∥ ≪ l ∥ / R E ∥ / B ∥ ∼ β E_{\parallel}/B_{\parallel}\sim\beta E ∥ / B ∥ ∼ β β ≪ l ∥ / R \beta\ll l_{\parallel}/R β ≪ l ∥ / R E ⊥ ∼ B ⊥ E_{\perp}\sim B_{\perp} E ⊥ ∼ B ⊥ η v i r ≈ η c o \eta_{\rm vir}\approx\eta_{\rm co} η vir ≈ η co ∣ j v i r ∣ ≪ ∣ j c o ∣ |j_{\rm vir}|\ll|j_{\rm co}| ∣ j vir ∣ ≪ ∣ j co ∣ l ∥ , min ∼ c / ω p ≪ R l_{\parallel,\min}\sim c/\omega_{\rm p}\ll R l ∥ , m i n ∼ c / ω p ≪ R (76) scale as R / l ∥ R/l_{\parallel} R / l ∥ E ∥ ≫ E ⊥ E_{\parallel}\gg E_{\perp} E ∥ ≫ E ⊥ B ⊥ ≈ 0 B_{\perp}\approx0 B ⊥ ≈ 0 δ j 1 ∝ ∇ × ( β × E ) \delta\mathbf{j}_1\propto\nabla\times(\boldsymbol{\beta}\times\mathbf{E}) δ j 1 ∝ ∇ × ( β × E ) δ j 1 ⊥ h ^ \delta\mathbf{j}_1\perp\hat{h} δ j 1 ⊥ h ^
Analysis of the Discharge Model
Here we present a self-contained dynamical-systems analysis of the coupled ODE system (Eqs. (46) –(49) ). We define the shorthand g 3 ≡ ⟨ γ − 3 ⟩ g_3\equiv\langle\gamma^{-3}\rangle g 3 ≡ ⟨ γ − 3 ⟩
The ordering of timescales given in Eq. (40) reveals three well-separated regimes among the ODE timescales: t a c c t_{\rm acc} t acc t ± ∼ t x t_{\pm}\sim t_x t ± ∼ t x t e s c t_{\rm esc} t esc
Block Structure and Linear Stability
The system (Eqs. (46) –(49) ) can be expressed in quasi-linear form as x ˙ = J ( x ) x + s \dot{\mathbf{x}}=\mathcal{J}(\mathbf{x})\mathbf{x}+\mathbf{s} x ˙ = J ( x ) x + s x = ( m ± , m x , ϵ , j ) T \mathbf{x}=(m_{\pm},m_x,\epsilon,j)^{\mathsf{T}} x = ( m ± , m x , ϵ , j ) T s = ( 0 , 0 , α / t e s c , 0 ) T \mathbf{s}=(0,0,\alpha/t_{\rm esc},0)^{\mathsf{T}} s = ( 0 , 0 , α / t esc , 0 ) T
J = ( − 1 t e s c 2 t ± 0 0 H t x − 1 t ± − 1 t e s c 0 0 0 0 0 − 1 t e s c 0 0 g 3 m ± t a c c − 1 t e s c ) . (79) \mathcal{J}=
\begin{pmatrix}
-\frac{1}{t_{\rm esc}} & \frac{2}{t_{\pm}} & 0 & 0 \\
\frac{\mathcal{H}}{t_x} & -\frac{1}{t_{\pm}}-\frac{1}{t_{\rm esc}} & 0 & 0 \\
0 & 0 & 0 & -\frac{1}{t_{\rm esc}} \\
0 & 0 & \frac{g_3m_{\pm}}{t_{\rm acc}} & -\frac{1}{t_{\rm esc}}
\end{pmatrix}. \tag*{(79)} J = − t esc 1 t x H 0 0 t ± 2 − t ± 1 − t esc 1 0 0 0 0 0 t acc g 3 m ± 0 0 − t esc 1 − t esc 1 . (79) The Heaviside gate H ≡ H ( ∣ ϵ ∣ − 1 ) \mathcal{H}\equiv\mathcal{H}(|\epsilon|-1) H ≡ H ( ∣ ϵ ∣ − 1 ) m ± m_{\pm} m ± H \mathcal{H} H 2 × 2 2\times2 2 × 2
Screened regime (H = 0 \mathcal{H}=0 H = 0 Setting H = 0 \mathcal{H}=0 H = 0 (79) gives the pair-photon block of the coefficient matrix
J s c r = ( − 1 / t e s c 2 / t ± 0 − 1 / t ± − 1 / t e s c ) . (80) \mathcal{J}_{\rm scr}=
\begin{pmatrix}
-1/t_{\rm esc} & 2/t_{\pm} \\
0 & -1/t_{\pm}-1/t_{\rm esc}
\end{pmatrix}. \tag*{(80)} J scr = ( − 1/ t esc 0 2/ t ± − 1/ t ± − 1/ t esc ) . (80) This upper-triangular matrix has eigenvalues − 1 / t e s c -1/t_{\rm esc} − 1/ t esc − ( t e s c + t ± ) / ( t e s c t ± ) -(t_{\rm esc}+t_{\pm})/(t_{\rm esc}t_{\pm}) − ( t esc + t ± ) / ( t esc t ± )
J f l d = ( 0 − 1 / t e s c g 3 m ± / t a c c − 1 / t e s c ) (81) \mathcal{J}_{\rm fld}=
\begin{pmatrix}
0 & -1/t_{\rm esc} \\
g_3m_{\pm}/t_{\rm acc} & -1/t_{\rm esc}
\end{pmatrix} \tag*{(81)} J fld = ( 0 g 3 m ± / t acc − 1/ t esc − 1/ t esc ) (81) has the eigenvalues
− 1 2 t e s c ± 1 4 t e s c 2 − g 3 m ± t e s c t a c c . (82) -\frac{1}{2t_{\rm esc}}\pm\sqrt{\frac{1}{4t_{\rm esc}^{2}}-\frac{g_3m_{\pm}}{t_{\rm esc}t_{\rm acc}}}. \tag*{(82)} − 2 t esc 1 ± 4 t esc 2 1 − t esc t acc g 3 m ± . (82) whose real part is always negative. These eigenvalues become complex conjugates, corresponding to damped oscillations, if
4 g 3 m ± t a c c > 1 t e s c . (83) \frac{4g_3m_{\pm}}{t_{\rm acc}}>\frac{1}{t_{\rm esc}}. \tag*{(83)} t acc 4 g 3 m ± > t esc 1 . (83) Thus, all modes decay exponentially and the screened state acts as a local attractor. Pair and photon perturbations decay exponentially, while the current relaxes towards α \alpha α
Unscreened regime (H = 1 \mathcal{H}=1 H = 1 During ignition, when t ≪ t e s c t\ll t_{\rm esc} t ≪ t esc
J u n s c r = ( 0 2 / t ± 1 / t x − 1 / t ± ) . (84) \mathcal{J}_{\rm unscr}=
\begin{pmatrix}
0 & 2/t_{\pm} \\
1/t_x & -1/t_{\pm}
\end{pmatrix}. \tag*{(84)} J unscr = ( 0 1/ t x 2/ t ± − 1/ t ± ) . (84) This matrix has a positive eigenvalue Γ 1 > 0 \Gamma_1>0 Γ 1 > 0
Together with particle escape and subsequent field recovery, the alternation between cascaded growth and screening produces the limit-cycle behavior discussed in the main text.
Relaxation Oscillation Structure
The discharge cycle has fast phases (ignition and saturation, timescale ∼ t ± \sim t_{\pm} ∼ t ± ∼ t e s c \sim t_{\rm esc} ∼ t esc t ± / t e s c = 10 − 3 t_{\pm}/t_{\rm esc}=10^{-3} t ± / t esc = 1 0 − 3 m ± m_{\pm} m ± ∼ 1 \sim1 ∼ 1 ∼ M ± \sim M_{\pm} ∼ M ±
The pair-photon sector ( m ± , m x ) (m_{\pm},m_x) ( m ± , m x ) ( ϵ , j ) (\epsilon,j) ( ϵ , j )
Current generation: the term m ± g 3 / t a c c m_{\pm}g_3/t_{\rm acc} m ± g 3 / t acc j j j
Heaviside gate: the factor H ( ∣ ϵ ∣ − 1 ) \mathcal{H}(|\epsilon|-1) H ( ∣ ϵ ∣ − 1 ) m ˙ x \dot{m}_x m ˙ x
This coupling is asymmetric in timescale: pairs evolve on t ± t_{\pm} t ± ω o s c − 1 \omega_{\rm osc}^{-1} ω osc − 1 t e s c t_{\rm esc} t esc
During decay, ε ˙ ≈ α / t e s c \dot{\varepsilon} \approx\alpha/t_{\mathrm{esc}} ε ˙ ≈ α / t esc ε ( t ) = α t / t e s c \varepsilon(t) = \alpha t/t_{\mathrm{esc}} ε ( t ) = α t / t esc ∣ ε ∣ = 1 |\varepsilon| = 1 ∣ ε ∣ = 1 t e s c / α t_{\mathrm{esc}}/\alpha t esc / α t c y c l e ≈ t e s c ln ( M ± / ∣ α ∣ ) t_{\mathrm{cycle}} \approx t_{\mathrm{esc}}\ln(M_{\pm}/|\alpha|) t cycle ≈ t esc ln ( M ± /∣ α ∣ ) (67) ).
Ignition Phase: Eigenvalue Analysis and Analytic Solutions
For t ≪ t e s c t \ll t_{\mathrm{esc}} t ≪ t esc H = 1 \mathcal{H} = 1 H = 1 2 × 2 2 \times2 2 × 2 J u n s c r \mathcal{J}_{\mathrm{unscr}} J unscr (84) ):
m ˙ ± = 2 m x t ± (85) \dot{m}_{\pm} = \frac{2m_x}{t_{\pm}} \tag*{(85)} m ˙ ± = t ± 2 m x (85) m ˙ x = − m x t ± + m + t x . (86) \dot{m}_x = -\frac{m_x}{t_{\pm}} + \frac{m_+}{t_x}. \tag*{(86)} m ˙ x = − t ± m x + t x m + . (86) The characteristic equation is
Γ 2 + Γ t ± − 2 t ± t x = 0 , (87) \Gamma^2 + \frac{\Gamma}{t_{\pm}} - \frac{2}{t_{\pm}t_x} = 0, \tag*{(87)} Γ 2 + t ± Γ − t ± t x 2 = 0 , (87) with eigenvalues
Γ 1 , 2 = − 1 ± 1 + 8 t ± / t x 2 t ± . (88) \Gamma_{1,2} = \frac{-1 \pm\sqrt{1+8t_{\pm}/t_x}}{2t_{\pm}}. \tag*{(88)} Γ 1 , 2 = 2 t ± − 1 ± 1 + 8 t ± / t x . (88) The positive root gives the growth rate Γ 1 = 1 / ( 2 t i g n ) − 1 / ( 2 t ± ) \Gamma_1 = 1/(2t_{\mathrm{ign}}) - 1/(2t_{\pm}) Γ 1 = 1/ ( 2 t ign ) − 1/ ( 2 t ± )
t i g n ≡ t ± t x 8 t ± + t x . (89) t_{\mathrm{ign}} \equiv\frac{t_{\pm}\sqrt{t_x}}{\sqrt{8t_{\pm}+t_x}}. \tag*{(89)} t ign ≡ 8 t ± + t x t ± t x . (89) For t ± ∼ t x t_{\pm} \sim t_x t ± ∼ t x t i g n ≈ t ± / 3 t_{\mathrm{ign}} \approx t_{\pm}/3 t ign ≈ t ± /3 Γ 1 = 1 / t ± \Gamma_1 = 1/t_{\pm} Γ 1 = 1/ t ± Γ 2 = − 2 / t ± \Gamma_2 = -2/t_{\pm} Γ 2 = − 2/ t ± Γ 1 \Gamma_1 Γ 1 m x / m ± = Γ 1 t ± / 2 m_x/m_{\pm} = \Gamma_1t_{\pm}/2 m x / m ± = Γ 1 t ± /2 t ± = t x t_{\pm} = t_x t ± = t x 1 / 2 1/2 1/2
With initial conditions m ± ( 0 ) = m ± , 0 m_{\pm}(0) = m_{\pm,0} m ± ( 0 ) = m ± , 0 m x ( 0 ) = 0 m_x(0) = 0 m x ( 0 ) = 0
m ± ( t ) = m ± , 0 exp ( − t 2 t ± − t 2 t i g n ) { t i g n 2 t ± [ exp ( t t i g n ) − 1 ] + 1 2 [ exp ( t t i g n ) + 1 ] } (90) m_{\pm}(t) = m_{\pm,0}\exp\left(-\frac{t}{2t_{\pm}}-\frac{t}{2t_{\mathrm{ign}}}\right)
\left\{\frac{t_{\mathrm{ign}}}{2t_{\pm}}\left[\exp\left(\frac{t}{t_{\mathrm{ign}}}\right)-1\right]+\frac{1}{2}\left[\exp\left(\frac{t}{t_{\mathrm{ign}}}\right)+1\right]\right\} \tag*{(90)} m ± ( t ) = m ± , 0 exp ( − 2 t ± t − 2 t ign t ) { 2 t ± t ign [ exp ( t ign t ) − 1 ] + 2 1 [ exp ( t ign t ) + 1 ] } (90) and
m x ( t ) = m ± , 0 t i g n t x exp ( − t 2 t ± − t 2 t i g n ) [ exp ( t t i g n ) − 1 ] . (91) m_x(t) = m_{\pm,0}\frac{t_{\mathrm{ign}}}{t_x}\exp\left(-\frac{t}{2t_{\pm}}-\frac{t}{2t_{\mathrm{ign}}}\right)\left[\exp\left(\frac{t}{t_{\mathrm{ign}}}\right)-1\right]. \tag*{(91)} m x ( t ) = m ± , 0 t x t ign exp ( − 2 t ± t − 2 t ign t ) [ exp ( t ign t ) − 1 ] . (91) For t ≳ t i g n t \gtrsim t_{\mathrm{ign}} t ≳ t ign t ≳ t ± t \gtrsim t_{\pm} t ≳ t ± t i g n ∼ t ± t_{\mathrm{ign}} \sim t_{\pm} t ign ∼ t ±
m ± ( t ) ∝ m x ( t ) ∝ exp ( t 2 t i g n − t 2 t ± ) . (92) m_{\pm}(t) \propto m_x(t) \propto\exp\left(\frac{t}{2t_{\mathrm{ign}}}-\frac{t}{2t_{\pm}}\right). \tag*{(92)} m ± ( t ) ∝ m x ( t ) ∝ exp ( 2 t ign t − 2 t ± t ) . (92) For t ± = t x t_{\pm} = t_x t ± = t x exp ( t / t ± ) \exp(t/t_{\pm}) exp ( t / t ± )
J u n s c r f u l l = ( − 1 / t e s c 2 / t ± 1 / t x − 1 / t ± − 1 / t e s c ) . (93) \mathcal{J}_{\mathrm{unscr}}^{\mathrm{full}} =
\begin{pmatrix}
-1/t_{\mathrm{esc}} & 2/t_{\pm}\\
1/t_x & -1/t_{\pm}-1/t_{\mathrm{esc}}
\end{pmatrix}. \tag*{(93)} J unscr full = ( − 1/ t esc 1/ t x 2/ t ± − 1/ t ± − 1/ t esc ) . (93) The eigenvalues shift by − 1 / t e s c -1/t_{\mathrm{esc}} − 1/ t esc J u n s c r \mathcal{J}_{\mathrm{unscr}} J unscr Γ 1 \Gamma_1 Γ 1 ≈ 0.1 % \approx0.1\% ≈ 0.1%
Growth requires the shifted growth eigenvalue to remain positive, equivalently det ( J u n s c r f u l l ) < 0 \det(\mathcal{J}_{\mathrm{unscr}}^{\mathrm{full}}) < 0 det ( J unscr full ) < 0
t x < t x , c r i t = 2 t e s c 2 t e s c + t ± ≈ 2 t e s c , (94) t_x < t_{x,\mathrm{crit}} = \frac{2t_{\mathrm{esc}}^2}{t_{\mathrm{esc}}+t_{\pm}} \approx2t_{\mathrm{esc}}, \tag*{(94)} t x < t x , crit = t esc + t ± 2 t esc 2 ≈ 2 t esc , (94) a bifurcation on the Heaviside switching boundary. The fiducial value t x / t e s c ≪ 2 t_x/t_{\mathrm{esc}} \ll2 t x / t esc ≪ 2
Saturation and Decay Phases
Saturation. When H → 0 \mathcal{H} \to0 H → 0
m ˙ ± = 2 m x t ± − m ± t e s c , (95) \dot{m}_{\pm} = 2\frac{m_x}{t_{\pm}}-\frac{m_{\pm}}{t_{\mathrm{esc}}}, \tag*{(95)} m ˙ ± = 2 t ± m x − t esc m ± , (95) m ˙ x = − m x t ± − m x t e s c . (96) \dot{m}_x = -\frac{m_x}{t_{\pm}}-\frac{m_x}{t_{\mathrm{esc}}}. \tag*{(96)} m ˙ x = − t ± m x − t esc m x . (96) The photon solution is
m x ( t ) = M x ′ exp ( − t t e s c − t t ± ) , (97) m_x(t) = M'_x\exp\left(-\frac{t}{t_{\mathrm{esc}}}-\frac{t}{t_{\pm}}\right), \tag*{(97)} m x ( t ) = M x ′ exp ( − t esc t − t ± t ) , (97) and the pair solution is
m ± ( t ) = m x ( t ) { exp ( t t ± ) [ M ± ′ M x ′ + 2 ] − 2 } , (98) m_{\pm}(t) = m_x(t)\left\{\exp\left(\frac{t}{t_{\pm}}\right)\left[\frac{M'_{\pm}}{M'_x}+2\right]-2\right\}, \tag*{(98)} m ± ( t ) = m x ( t ) { exp ( t ± t ) [ M x ′ M ± ′ + 2 ] − 2 } , (98) where M x ′ M'_x M x ′ M ± ′ M'_{\pm} M ± ′ m ± ( t ) m_{\pm}(t) m ± ( t ) t d = t ± ln { 2 M x ′ ( t e s c + t ± ) / [ t ± ( M ± ′ + 2 M x ′ ) ] } t_d = t_{\pm}\ln\{2M'_x(t_{\mathrm{esc}}+t_{\pm})/[t_{\pm}(M'_{\pm}+2M'_x)]\} t d = t ± ln { 2 M x ′ ( t esc + t ± ) / [ t ± ( M ± ′ + 2 M x ′ )]} 2 M x ′ / M ± ′ > t ± / t e s c 2M'_x/M'_{\pm} > t_{\pm}/t_{\mathrm{esc}} 2 M x ′ / M ± ′ > t ± / t esc
Decay. At late times, t ≫ t ± t \gg t_{\pm} t ≫ t ± (98) reduces to
m ± ( t ) ≈ M ± exp ( − t t e s c ) , (99) m_{\pm}(t) \approx M_{\pm}\exp\left(-\frac{t}{t_{\mathrm{esc}}}\right), \tag*{(99)} m ± ( t ) ≈ M ± exp ( − t esc t ) , (99) where the peak multiplicity is M ± ≈ M ± ′ + 2 M x ′ ≈ 2 M ± ′ M_{\pm} \approx M'_{\pm}+2M'_x \approx2M'_{\pm} M ± ≈ M ± ′ + 2 M x ′ ≈ 2 M ± ′ t x ∼ t ± t_x \sim t_{\pm} t x ∼ t ±
Current solution. The electric field oscillates on timescale ω o s c − 1 ≪ t e s c \omega_{\mathrm{osc}}^{-1} \ll t_{\mathrm{esc}} ω osc − 1 ≪ t esc ε \varepsilon ε
j ( t ) = e − t / t e s c ( j d , 0 + M ± ϵ g 3 t t a c c ) , (100) j(t) = e^{-t/t_{\mathrm{esc}}}\left(j_{\mathrm{d},0}+M_{\pm}\epsilon g_3\frac{t}{t_{\mathrm{acc}}}\right), \tag*{(100)} j ( t ) = e − t / t esc ( j d , 0 + M ± ϵ g 3 t acc t ) , (100) where j d , 0 ≡ j ( t = 0 ) j_{\mathrm{d},0} \equiv j(t=0) j d , 0 ≡ j ( t = 0 ) t = t e s c t=t_{\mathrm{esc}} t = t esc j d , 0 j_{\mathrm{d},0} j d , 0
j ( t e s c ) ≈ m ± ( t e s c ) ε g 3 t e s c t a c c . (101) j(t_{\mathrm{esc}}) \approx m_{\pm}(t_{\mathrm{esc}})\,\varepsilon\frac{g_3t_{\mathrm{esc}}}{t_{\mathrm{acc}}}. \tag*{(101)} j ( t esc ) ≈ m ± ( t esc ) ε t acc g 3 t esc . (101) Field-Current Oscillator
During the oscillation and decay phases, m ± m_{\pm} m ± j j j ( ε , j ) (\varepsilon,j) ( ε , j )
ε ¨ + ε ˙ t e s c + ω o s c 2 ε = α t e s c 2 . (102) \ddot{\varepsilon}+\frac{\dot{\varepsilon}}{t_{\mathrm{esc}}}+\omega_{\mathrm{osc}}^2\varepsilon= \frac{\alpha}{t_{\mathrm{esc}}^2}. \tag*{(102)} ε ¨ + t esc ε ˙ + ω osc 2 ε = t esc 2 α . (102) with natural frequency, damping ratio, and quality factor
ω o s c ( t ) = g 3 m ± ( t ) t a c c t e s c , ζ = 1 2 ω o s c t e s c , Q = ω o s c t e s c . (103) \omega_{\rm osc}(t)=\sqrt{\frac{g_3m_\pm(t)}{t_{\rm acc}t_{\rm esc}}},\qquad\zeta=\frac{1}{2\omega_{\rm osc}t_{\rm esc}},\qquad Q=\omega_{\rm osc}t_{\rm esc}. \tag*{(103)} ω osc ( t ) = t acc t esc g 3 m ± ( t ) , ζ = 2 ω osc t esc 1 , Q = ω osc t esc . (103) The amplitude damping time is t d a m p = 2 t e s c t_{\rm damp}=2t_{\rm esc} t damp = 2 t esc m ± m_\pm m ± g 3 g_3 g 3 α t a c c / ( m ± ( t ) g 3 t e s c ) ≪ 1 \alpha t_{\rm acc}/(m_\pm(t)g_3t_{\rm esc})\ll1 α t acc / ( m ± ( t ) g 3 t esc ) ≪ 1
Underdamped condition. The system oscillates (ζ < 1 \zeta<1 ζ < 1 m ± > t a c c / ( 4 g 3 t e s c ) m_\pm>t_{\rm acc}/(4g_3t_{\rm esc}) m ± > t acc / ( 4 g 3 t esc ) g 3 = γ min − 3 ∼ 10 − 7 g_3=\gamma_{\min}^{-3}\sim10^{-7} g 3 = γ m i n − 3 ∼ 1 0 − 7 m ± ∼ 0.1 m_\pm\sim0.1 m ± ∼ 0.1 Q ∼ 100 Q\sim100 Q ∼ 100
Frequency chirp during decay. For a constant g 3 g_3 g 3 m ± ( t ) ∝ e − t / t e s c m_\pm(t)\propto e^{-t/t_{\rm esc}} m ± ( t ) ∝ e − t / t esc (65) ), the oscillation frequency chirps downward:
ω o s c ( t ) ∝ m ± ( t ) ∝ e − t / ( 2 t e s c ) . (104) \omega_{\rm osc}(t)\propto\sqrt{m_\pm(t)}\propto e^{-t/(2t_{\rm esc})}. \tag*{(104)} ω osc ( t ) ∝ m ± ( t ) ∝ e − t / ( 2 t esc ) . (104) When the concurrent growth of g 3 g_3 g 3 g 3 g_3 g 3 m ± m_\pm m ± ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩
Peak Multiplicity
We derive an estimate for M ± M_\pm M ± ε = 1 \varepsilon=1 ε = 1
Quasi-static current. Setting j ˙ = 0 \dot{j}=0 j ˙ = 0 ε = 1 \varepsilon=1 ε = 1 (49) defines the quasi-static current
j q s ≡ m ± g 3 t e s c t a c c . (105) j_{\rm qs}\equiv\frac{m_\pm g_3t_{\rm esc}}{t_{\rm acc}}. \tag*{(105)} j qs ≡ t acc m ± g 3 t esc . (105) Current lag. During ignition, m ± ( t ) = m ± , 0 e Γ 1 t m_\pm(t)=m_{\pm,0}e^{\Gamma_1t} m ± ( t ) = m ± , 0 e Γ 1 t ε ≈ 1 \varepsilon\approx1 ε ≈ 1 (105) , the current equation becomes j ˙ = − ( j − j q s ) / t e s c \dot{j}=-(j-j_{\rm qs})/t_{\rm esc} j ˙ = − ( j − j qs ) / t esc t e s c t_{\rm esc} t esc j p ( t ) = C e Γ 1 t j_p(t)=Ce^{\Gamma_1t} j p ( t ) = C e Γ 1 t
C = j q s , 0 Γ 1 t e s c + 1 , (106) C=\frac{j_{{\rm qs},0}}{\Gamma_1t_{\rm esc}+1}, \tag*{(106)} C = Γ 1 t esc + 1 j qs , 0 , (106) where j q s , 0 = m ± , 0 g 3 t e s c / t a c c j_{{\rm qs},0}=m_{\pm,0}g_3t_{\rm esc}/t_{\rm acc} j qs , 0 = m ± , 0 g 3 t esc / t acc
j p j q s = C j q s , 0 = 1 Γ 1 t e s c + 1 ≈ t ± t e s c ≪ 1. (107) \frac{j_p}{j_{\rm qs}}=\frac{C}{j_{{\rm qs},0}}=\frac{1}{\Gamma_1t_{\rm esc}+1}\approx\frac{t_\pm}{t_{\rm esc}}\ll1. \tag*{(107)} j qs j p = j qs , 0 C = Γ 1 t esc + 1 1 ≈ t esc t ± ≪ 1. (107) Thus, the current remains suppressed relative to j q s j_{\rm qs} j qs
Screening from the ε \varepsilon ε At the onset of ignition, the current has relaxed to j 0 ≪ α j_0\ll\alpha j 0 ≪ α j < α j<\alpha j < α ε \varepsilon ε ε ˙ = ( α − j ) / t e s c \dot{\varepsilon}=(\alpha-j)/t_{\rm esc} ε ˙ = ( α − j ) / t esc
ε ( t ) ≈ 1 + ( α − j 0 ) t t e s c − C ( e Γ 1 t − 1 ) Γ 1 t e s c , (108) \varepsilon(t)\approx1+\frac{(\alpha-j_0)t}{t_{\rm esc}}-\frac{C\left(e^{\Gamma_1t}-1\right)}{\Gamma_1t_{\rm esc}}, \tag*{(108)} ε ( t ) ≈ 1 + t esc ( α − j 0 ) t − Γ 1 t esc C ( e Γ 1 t − 1 ) , (108) where the first correction is a linear rise (current deficit) and the second is the exponentially growing screening term. Pair creation ceases when ε \varepsilon ε H ( ∣ ε ∣ − 1 ) → 0 \mathcal{H}(|\varepsilon|-1)\to0 H ( ∣ ε ∣ − 1 ) → 0 ε ( t ∗ ) = 1 \varepsilon(t^*)=1 ε ( t ∗ ) = 1 t ∗ = ln M ± ′ / Γ 1 t^*=\ln M'_\pm/\Gamma_1 t ∗ = ln M ± ′ / Γ 1 M ± ′ = m ± , 0 e Γ 1 t ∗ M'_\pm=m_{\pm,0}e^{\Gamma_1t^*} M ± ′ = m ± , 0 e Γ 1 t ∗ m ± , 0 = 1 m_{\pm,0}=1 m ± , 0 = 1
M ± ′ = M 0 ln M ± ′ + 1 , M 0 ≡ ( α − j 0 ) Γ 1 t a c c g 3 . (109) M'_\pm=M_0\ln M'_\pm+1,\qquad M_0\equiv\frac{(\alpha-j_0)\Gamma_1t_{\rm acc}}{g_3}. \tag*{(109)} M ± ′ = M 0 ln M ± ′ + 1 , M 0 ≡ g 3 ( α − j 0 ) Γ 1 t acc . (109) The current deficit is α − j 0 = O ( 1 ) \alpha-j_0=O(1) α − j 0 = O ( 1 ) M 0 ≈ Γ 1 t a c c / g 3 M_0\approx\Gamma_1t_{\rm acc}/g_3 M 0 ≈ Γ 1 t acc / g 3 α \alpha α M ± ∼ 10 3 M_\pm\sim10^3 M ± ∼ 1 0 3 10 4 10^4 1 0 4 ln M ± ′ ∼ 7 \ln M'_\pm\sim7 ln M ± ′ ∼ 7 9 9 9 M ± ≈ 2 M ± ′ ∼ 10 M 0 M_\pm\approx2M'_\pm\sim10M_0 M ± ≈ 2 M ± ′ ∼ 10 M 0 t ± = t x t_\pm=t_x t ± = t x Γ 1 = 1 / t ± \Gamma_1=1/t_\pm Γ 1 = 1/ t ±
M ± ∼ 10 t a c c t x g 3 . (110) M_\pm\sim10\frac{t_{\rm acc}}{t_xg_3}. \tag*{(110)} M ± ∼ 10 t x g 3 t acc . (110) Note that the g 3 g_3 g 3 ⟨ γ − 3 ⟩ i g n \langle\gamma^{-3}\rangle_{\rm ign} ⟨ γ − 3 ⟩ ign ⟨ γ − 3 ⟩ i n i t ∼ γ r a d − 3 \langle\gamma^{-3}\rangle_{\rm init}\sim\gamma_{\rm rad}^{-3} ⟨ γ − 3 ⟩ init ∼ γ rad − 3 ⟨ γ − 3 ⟩ s c r \langle\gamma^{-3}\rangle_{\rm scr} ⟨ γ − 3 ⟩ scr
Empirical fitting formula. A systematic scan of 36 ODE solutions across five parameter sweeps (g 3 g_3 g 3 t ± t_\pm t ± t x t_x t x t a c c t_{\rm acc} t acc α \alpha α t a c c / ( g 3 t x ) t_{\rm acc}/(g_3t_x) t acc / ( g 3 t x ) t ± t_\pm t ± Γ 1 \Gamma_1 Γ 1 t ± t_\pm t ± ∝ Γ 1 \propto\Gamma_1 ∝ Γ 1 1 + 2 m x / m ± = 1 + Γ 1 t ± 1+2m_x/m_\pm=1+\Gamma_1t_\pm 1 + 2 m x / m ± = 1 + Γ 1 t ± ∝ Γ 1 ( 1 + Γ 1 t ± ) ≈ 2 / t x \propto\Gamma_1(1+\Gamma_1t_\pm)\approx2/t_x ∝ Γ 1 ( 1 + Γ 1 t ± ) ≈ 2/ t x
M ± ≈ 11 t a c c t x g 3 α 1 / 4 , (111) M_\pm\approx11\frac{t_{\rm acc}}{t_xg_3}\alpha^{1/4}, \tag*{(111)} M ± ≈ 11 t x g 3 t acc α 1/4 , (111) which tracks the ODE solutions to ≲ 30 % \lesssim30\% ≲ 30%
On the Estimation of γ − 3 \gamma^{-3} γ − 3
Approximating the moment ⟨ γ − 3 ⟩ \langle\gamma^{-3}\rangle ⟨ γ − 3 ⟩ ⟨ γ ⟩ − 3 \langle\gamma\rangle^{-3} ⟨ γ ⟩ − 3 γ ˉ = ⟨ γ ⟩ \bar{\gamma}=\langle\gamma\rangle γ ˉ = ⟨ γ ⟩ δ γ = γ − γ ˉ \delta\gamma=\gamma-\bar{\gamma} δ γ = γ − γ ˉ [37 ] ),
⟨ γ − 3 ⟩ = ⟨ ( γ ˉ + δ γ ) − 3 ⟩ = γ ˉ − 3 ( 1 + 6 ⟨ δ γ 2 ⟩ γ ˉ 2 − 10 ⟨ δ γ 3 ⟩ γ ˉ 3 + ⋯ ) (112) \langle\gamma^{-3}\rangle=\left\langle(\bar{\gamma}+\delta\gamma)^{-3}\right\rangle
=\bar{\gamma}^{-3}\left(1+6\frac{\langle\delta\gamma^2\rangle}{\bar{\gamma}^2}-10\frac{\langle\delta\gamma^3\rangle}{\bar{\gamma}^3}+\cdots\right) \tag*{(112)} ⟨ γ − 3 ⟩ = ⟨ ( γ ˉ + δ γ ) − 3 ⟩ = γ ˉ − 3 ( 1 + 6 γ ˉ 2 ⟨ δ γ 2 ⟩ − 10 γ ˉ 3 ⟨ δ γ 3 ⟩ + ⋯ ) (112) shows that large relative variance in relativistic plasmas renders the series slow to converge. Jensen’s inequality further guarantees that ⟨ γ − 3 ⟩ ≥ ⟨ γ ⟩ − 3 \langle\gamma^{-3}\rangle\geq\langle\gamma\rangle^{-3} ⟨ γ − 3 ⟩ ≥ ⟨ γ ⟩ − 3 f ( γ ) ≈ Θ − 1 e − γ / Θ f(\gamma)\approx\Theta^{-1}e^{-\gamma/\Theta} f ( γ ) ≈ Θ − 1 e − γ /Θ ⟨ γ − 3 ⟩ ≈ ( 2 Θ ) − 1 \langle\gamma^{-3}\rangle\approx(2\Theta)^{-1} ⟨ γ − 3 ⟩ ≈ ( 2Θ ) − 1 ⟨ γ ⟩ − 3 ≈ Θ − 3 \langle\gamma\rangle^{-3}\approx\Theta^{-3} ⟨ γ ⟩ − 3 ≈ Θ − 3 ⟨ γ − 3 ⟩ / ⟨ γ ⟩ − 3 ≈ Θ 2 / 2 \langle\gamma^{-3}\rangle/\langle\gamma\rangle^{-3}\approx\Theta^2/2 ⟨ γ − 3 ⟩ / ⟨ γ ⟩ − 3 ≈ Θ 2 /2 ∼ 8 × 10 4 \sim8\times10^4 ∼ 8 × 1 0 4 ( w h e r e (where ( w h er e a n d and an d a r e t h e p o w e r l a w i n d i c e s , a n d are the powerlaw indices, and a r e t h e p o w er l a w in d i ces , an d i s t h e s m o o t h i n g f a c t o r ) , g i v e s is the smoothing factor), gives i s t h es m oo t hin g f a c t or ) , g i v es
This paper was built using the Open Journal of Astrophysics \LaTeX\ template. The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike. Learn more at http://astro.theoj.org .
