sweet_parker.py

Download
python 332 lines 8.9 KB
  1#!/usr/bin/env python3
  2"""
  3Sweet-Parker Reconnection Model
  4
  5Analyzes steady-state Sweet-Parker current sheet reconnection.
  6
  7Key results:
  8- Reconnection rate ~ S^(-1/2) where S is Lundquist number
  9- Too slow for solar flares (S ~ 10^12)
 10- Comparison with Petschek model
 11
 12Author: Claude
 13Date: 2026-02-14
 14"""
 15
 16import numpy as np
 17import matplotlib.pyplot as plt
 18
 19
 20def sweet_parker_rate(S):
 21    """
 22    Sweet-Parker reconnection rate.
 23
 24    Rate = 1 / sqrt(S)
 25
 26    where S = L V_A / η is the Lundquist number.
 27
 28    Parameters
 29    ----------
 30    S : float or array
 31        Lundquist number
 32
 33    Returns
 34    -------
 35    rate : float or array
 36        Reconnection rate (normalized)
 37    """
 38    return 1.0 / np.sqrt(S)
 39
 40
 41def petschek_rate(S):
 42    """
 43    Petschek reconnection rate.
 44
 45    Rate = π / (8 ln(S))
 46
 47    Parameters
 48    ----------
 49    S : float or array
 50        Lundquist number
 51
 52    Returns
 53    -------
 54    rate : float or array
 55        Reconnection rate
 56    """
 57    return np.pi / (8 * np.log(S + 1))  # +1 to avoid log(0)
 58
 59
 60def sheet_parameters(S, L, B0, n, T):
 61    """
 62    Compute Sweet-Parker current sheet parameters.
 63
 64    Parameters
 65    ----------
 66    S : float
 67        Lundquist number
 68    L : float
 69        System size (m)
 70    B0 : float
 71        Magnetic field (T)
 72    n : float
 73        Density (m⁻³)
 74    T : float
 75        Temperature (eV)
 76
 77    Returns
 78    -------
 79    dict
 80        Sheet parameters
 81    """
 82    mu_0 = 4 * np.pi * 1e-7
 83    m_i = 1.67e-27  # Proton mass
 84
 85    # Alfven speed
 86    V_A = B0 / np.sqrt(mu_0 * n * m_i)
 87
 88    # Reconnection rate
 89    rate = sweet_parker_rate(S)
 90
 91    # Inflow velocity
 92    V_in = rate * V_A
 93
 94    # Sheet length
 95    sheet_length = L
 96
 97    # Sheet width
 98    delta = L / np.sqrt(S)
 99
100    # Outflow velocity
101    V_out = V_A
102
103    # Current density (approximate)
104    J = B0 / (mu_0 * delta)
105
106    # Resistive diffusion time
107    tau_resist = L**2 / (V_A * L / S)  # = S * tau_A
108
109    # Reconnection time
110    tau_reconnect = L / V_in
111
112    return {
113        'V_A': V_A,
114        'rate': rate,
115        'V_in': V_in,
116        'V_out': V_out,
117        'delta': delta,
118        'L': sheet_length,
119        'J': J,
120        'tau_resist': tau_resist,
121        'tau_reconnect': tau_reconnect
122    }
123
124
125def plot_rate_vs_lundquist():
126    """
127    Plot reconnection rate vs Lundquist number.
128    """
129    S = np.logspace(2, 14, 100)  # 10² to 10¹⁴
130
131    rate_SP = sweet_parker_rate(S)
132    rate_Petschek = petschek_rate(S)
133
134    # Hall MHD rate (approximately constant)
135    rate_Hall = 0.1 * np.ones_like(S)
136
137    # Plasmoid-mediated rate
138    rate_plasmoid = 0.01 * np.ones_like(S)
139
140    fig, ax = plt.subplots(figsize=(11, 7))
141
142    ax.loglog(S, rate_SP, 'b-', linewidth=3, label='Sweet-Parker ~ S^(-1/2)')
143    ax.loglog(S, rate_Petschek, 'r-', linewidth=3,
144             label='Petschek ~ π/(8 ln S)')
145    ax.axhline(rate_Hall[0], color='g', linestyle='--', linewidth=2.5,
146              label='Hall MHD ~ 0.1')
147    ax.axhline(rate_plasmoid[0], color='m', linestyle='--', linewidth=2.5,
148              label='Plasmoid-mediated ~ 0.01')
149
150    # Mark solar flare regime
151    ax.axvline(1e12, color='orange', linestyle=':', linewidth=2,
152              alpha=0.7, label='Solar flares (S~10¹²)')
153
154    # Mark observed rates
155    ax.fill_between([1e2, 1e14], 0.01, 0.1, alpha=0.2, color='yellow',
156                    label='Observed rates')
157
158    ax.set_xlabel('Lundquist number S', fontsize=13)
159    ax.set_ylabel('Reconnection rate (V_in / V_A)', fontsize=13)
160    ax.set_title('Magnetic Reconnection Rate vs Lundquist Number',
161                 fontsize=14, fontweight='bold')
162    ax.legend(fontsize=10, loc='best')
163    ax.grid(True, alpha=0.3, which='both')
164    ax.set_ylim([1e-7, 1])
165
166    plt.tight_layout()
167    plt.savefig('reconnection_rate_vs_S.png', dpi=150, bbox_inches='tight')
168    plt.show()
169
170
171def plot_sheet_geometry():
172    """
173    Visualize Sweet-Parker current sheet geometry.
174    """
175    # Sheet dimensions (normalized)
176    L = 1.0
177    S = 1e4
178    delta = L / np.sqrt(S)
179
180    fig, ax = plt.subplots(figsize=(12, 6))
181
182    # Draw sheet
183    sheet_x = [-L/2, L/2, L/2, -L/2, -L/2]
184    sheet_y = [-delta/2, -delta/2, delta/2, delta/2, -delta/2]
185    ax.fill(sheet_x, sheet_y, color='orange', alpha=0.5, label='Current sheet')
186    ax.plot(sheet_x, sheet_y, 'k-', linewidth=2)
187
188    # Draw magnetic field lines
189    y = np.linspace(-3*delta, 3*delta, 20)
190    for yi in y:
191        if np.abs(yi) > delta/2:
192            # Field lines outside sheet
193            x_left = np.linspace(-L, -L/2, 50)
194            x_right = np.linspace(L/2, L, 50)
195
196            # Approximate field lines
197            ax.plot(x_left, yi * np.ones_like(x_left), 'b-', linewidth=1, alpha=0.7)
198            ax.plot(x_right, yi * np.ones_like(x_right), 'b-', linewidth=1, alpha=0.7)
199
200    # Arrows for flow
201    # Inflow
202    ax.arrow(0, 2*delta, 0, -0.5*delta, head_width=0.02, head_length=0.1*delta,
203            fc='red', ec='red', linewidth=2)
204    ax.text(0.05, 2*delta, 'Inflow\n$V_{in}$', fontsize=11, color='red',
205           fontweight='bold')
206
207    # Outflow
208    ax.arrow(L/2, 0, 0.2*L, 0, head_width=0.1*delta, head_length=0.05*L,
209            fc='green', ec='green', linewidth=2)
210    ax.text(0.7*L, 0.3*delta, 'Outflow\n$V_{out} \\sim V_A$',
211           fontsize=11, color='green', fontweight='bold')
212
213    # Annotations
214    ax.annotate('', xy=(L/2, -delta), xytext=(-L/2, -delta),
215               arrowprops=dict(arrowstyle='<->', lw=2))
216    ax.text(0, -1.5*delta, f'L = {L:.1f}', fontsize=12, ha='center')
217
218    ax.annotate('', xy=(L/2+0.1, delta/2), xytext=(L/2+0.1, -delta/2),
219               arrowprops=dict(arrowstyle='<->', lw=2))
220    ax.text(L/2+0.2, 0, f'δ = L/√S\\n={delta:.3f}',
221           fontsize=11, ha='left')
222
223    ax.set_xlabel('x / L', fontsize=13)
224    ax.set_ylabel('y / L', fontsize=13)
225    ax.set_title(f'Sweet-Parker Current Sheet (S={S:.0e})',
226                 fontsize=14, fontweight='bold')
227    ax.set_aspect('equal')
228    ax.grid(True, alpha=0.3)
229    ax.legend(fontsize=11)
230
231    plt.tight_layout()
232    plt.savefig('sweet_parker_geometry.png', dpi=150, bbox_inches='tight')
233    plt.show()
234
235
236def plot_scaling_comparison():
237    """
238    Compare different reconnection timescales.
239    """
240    S_arr = np.logspace(2, 14, 100)
241
242    # Alfven time (reference)
243    tau_A = 1.0  # Normalized
244
245    # Resistive time
246    tau_resist = S_arr * tau_A
247
248    # Sweet-Parker reconnection time
249    tau_SP = np.sqrt(S_arr) * tau_A
250
251    # Petschek reconnection time
252    tau_Petschek = (8 * np.log(S_arr) / np.pi) * tau_A
253
254    # Fast (observed) reconnection
255    tau_fast = 10 * tau_A * np.ones_like(S_arr)
256
257    fig, ax = plt.subplots(figsize=(11, 7))
258
259    ax.loglog(S_arr, tau_resist / tau_A, 'k--', linewidth=2.5,
260             label='Resistive ~ S τ_A')
261    ax.loglog(S_arr, tau_SP / tau_A, 'b-', linewidth=3,
262             label='Sweet-Parker ~ √S τ_A')
263    ax.loglog(S_arr, tau_Petschek / tau_A, 'r-', linewidth=3,
264             label='Petschek ~ ln(S) τ_A')
265    ax.axhline(10, color='g', linestyle='--', linewidth=2.5,
266              label='Fast reconnection ~ 10 τ_A')
267
268    # Solar flare constraint
269    ax.axvline(1e12, color='orange', linestyle=':', linewidth=2,
270              alpha=0.7)
271    ax.fill_between([1e12, 1e14], 1, 100, alpha=0.2, color='yellow',
272                    label='Solar flare requirement')
273
274    ax.set_xlabel('Lundquist number S', fontsize=13)
275    ax.set_ylabel('Reconnection time / τ_A', fontsize=13)
276    ax.set_title('Reconnection Timescale Comparison',
277                 fontsize=14, fontweight='bold')
278    ax.legend(fontsize=10, loc='upper left')
279    ax.grid(True, alpha=0.3, which='both')
280
281    plt.tight_layout()
282    plt.savefig('reconnection_timescales.png', dpi=150, bbox_inches='tight')
283    plt.show()
284
285
286def main():
287    """Main execution function."""
288    print("Sweet-Parker Reconnection Model")
289    print("=" * 60)
290
291    # Example parameters (solar corona)
292    S = 1e12  # Lundquist number
293    L = 1e7  # 10,000 km
294    B0 = 0.01  # 100 G
295    n = 1e15  # m⁻³
296    T = 100  # eV
297
298    params = sheet_parameters(S, L, B0, n, T)
299
300    print(f"\nSolar corona parameters:")
301    print(f"  Lundquist number S: {S:.1e}")
302    print(f"  System size L: {L/1e3:.0f} km")
303    print(f"  Magnetic field B₀: {B0*1e4:.0f} G")
304    print(f"  Density n: {n:.1e} m⁻³")
305    print()
306
307    print(f"Sweet-Parker sheet:")
308    print(f"  Alfven speed: {params['V_A']/1e3:.0f} km/s")
309    print(f"  Reconnection rate: {params['rate']:.2e}")
310    print(f"  Inflow velocity: {params['V_in']/1e3:.1f} km/s")
311    print(f"  Sheet width δ: {params['delta']/1e3:.1f} km")
312    print(f"  Reconnection time: {params['tau_reconnect']:.1e} s ({params['tau_reconnect']/60:.1f} min)")
313    print()
314
315    print("Problem: Solar flares occur in ~100 s, but Sweet-Parker predicts ~weeks!")
316    print()
317
318    # Plots
319    print("Generating plots...")
320    plot_rate_vs_lundquist()
321    print("  Saved 'reconnection_rate_vs_S.png'")
322
323    plot_sheet_geometry()
324    print("  Saved 'sweet_parker_geometry.png'")
325
326    plot_scaling_comparison()
327    print("  Saved 'reconnection_timescales.png'")
328
329
330if __name__ == '__main__':
331    main()