1#!/usr/bin/env python3
2"""
3Magnetorotational Instability (MRI) - Shearing Box Analysis
4
5This module analyzes the magnetorotational instability (MRI), which is crucial
6for understanding angular momentum transport in accretion disks around compact
7objects (black holes, neutron stars, young stars).
8
9Key physics:
10- MRI destabilizes otherwise stable Keplerian disks when threaded by
11 weak magnetic fields
12- Growth rate: γ_max ≈ (3/4)Ω for ideal MHD
13- Requires weak vertical field: B_z > B_crit
14- Resistivity provides a cutoff at small scales
15
16The analysis uses local shearing box approximation with dispersion relation.
17
18Author: MHD Course Examples
19License: MIT
20"""
21
22import numpy as np
23import matplotlib.pyplot as plt
24from scipy.optimize import fsolve
25
26
27class MRIShearingBox:
28 """
29 MRI dispersion relation solver in shearing box approximation.
30
31 Attributes:
32 Omega (float): Angular velocity at fiducial radius
33 q (float): Shear parameter (q=3/2 for Keplerian)
34 B_z (float): Vertical magnetic field
35 rho (float): Density
36 eta (float): Magnetic diffusivity
37 nu (float): Kinematic viscosity
38 """
39
40 def __init__(self, Omega=1.0, q=1.5, B_z=0.1, rho=1.0,
41 eta=0.0, nu=0.0):
42 """
43 Initialize MRI shearing box model.
44
45 Parameters:
46 Omega (float): Angular velocity
47 q (float): Shear parameter (q=3/2 for Keplerian)
48 B_z (float): Vertical field strength
49 rho (float): Density
50 eta (float): Magnetic diffusivity
51 nu (float): Kinematic viscosity
52 """
53 self.Omega = Omega
54 self.q = q
55 self.B_z = B_z
56 self.rho = rho
57 self.eta = eta
58 self.nu = nu
59
60 # Derived quantities
61 self.kappa2 = 2 * Omega**2 * (2 - q) # Epicyclic frequency squared
62 self.v_A = B_z / np.sqrt(rho) # Alfvén velocity
63
64 # Magnetic Prandtl number
65 self.Pm = nu / eta if eta > 0 else np.inf
66
67 def dispersion_relation_ideal(self, gamma, k_z):
68 """
69 Ideal MHD dispersion relation for MRI.
70
71 (γ² + ν²_A k_z²)(γ² + κ²) - 4Ω²ν²_A k_z² = 0
72
73 Parameters:
74 gamma (complex): Growth rate
75 k_z (float): Vertical wavenumber
76
77 Returns:
78 complex: Residual
79 """
80 kappa2 = self.kappa2
81 Omega = self.Omega
82 v_A = self.v_A
83
84 term1 = gamma**2 + v_A**2 * k_z**2
85 term2 = gamma**2 + kappa2
86 coupling = 4 * Omega**2 * v_A**2 * k_z**2
87
88 return term1 * term2 - coupling
89
90 def dispersion_relation_resistive(self, gamma, k_z):
91 """
92 Resistive MHD dispersion relation.
93
94 Includes magnetic diffusion: adds -iηk² terms
95
96 Parameters:
97 gamma (complex): Growth rate
98 k_z (float): Vertical wavenumber
99
100 Returns:
101 complex: Residual
102 """
103 if self.eta == 0:
104 return self.dispersion_relation_ideal(gamma, k_z)
105
106 kappa2 = self.kappa2
107 Omega = self.Omega
108 v_A = self.v_A
109 eta = self.eta
110
111 # Modified dispersion relation with resistivity
112 gamma_eff = gamma + 1j * eta * k_z**2
113
114 term1 = gamma_eff**2 + v_A**2 * k_z**2
115 term2 = gamma_eff**2 + kappa2
116 coupling = 4 * Omega**2 * v_A**2 * k_z**2
117
118 return term1 * term2 - coupling
119
120 def find_growth_rate_ideal(self, k_z):
121 """
122 Find growth rate for given k_z (ideal MHD).
123
124 Parameters:
125 k_z (float): Vertical wavenumber
126
127 Returns:
128 float: Maximum growth rate (real part)
129 """
130 # Analytical solution exists for ideal MRI
131 kappa2 = self.kappa2
132 Omega = self.Omega
133 v_A = self.v_A
134
135 # MRI growth rate formula
136 omega_A2 = v_A**2 * k_z**2
137
138 discriminant = omega_A2 + kappa2
139 gamma2 = 0.5 * (omega_A2 + kappa2 -
140 np.sqrt((omega_A2 + kappa2)**2 - 16 * Omega**2 * omega_A2))
141
142 if gamma2 > 0:
143 return np.sqrt(gamma2)
144 else:
145 return 0.0
146
147 def find_growth_rate_resistive(self, k_z):
148 """
149 Find growth rate for given k_z (resistive MHD).
150
151 Uses numerical root finding.
152
153 Parameters:
154 k_z (float): Vertical wavenumber
155
156 Returns:
157 float: Maximum growth rate
158 """
159 # Initial guess from ideal solution
160 gamma_ideal = self.find_growth_rate_ideal(k_z)
161
162 # Solve dispersion relation numerically
163 def residual(gamma):
164 return abs(self.dispersion_relation_resistive(gamma, k_z))
165
166 # Search for maximum growth rate
167 gamma_range = np.linspace(0, gamma_ideal * 1.5, 100)
168 residuals = [residual(g) for g in gamma_range]
169 idx_min = np.argmin(residuals)
170
171 # Refine
172 try:
173 gamma_opt = fsolve(residual, gamma_range[idx_min])[0]
174 return max(gamma_opt, 0.0)
175 except:
176 return 0.0
177
178 def compute_growth_rate_spectrum(self, k_z_array):
179 """
180 Compute growth rate as function of k_z.
181
182 Parameters:
183 k_z_array (ndarray): Array of wavenumbers
184
185 Returns:
186 ndarray: Growth rates
187 """
188 gamma_array = np.zeros_like(k_z_array)
189
190 for i, k_z in enumerate(k_z_array):
191 if self.eta == 0:
192 gamma_array[i] = self.find_growth_rate_ideal(k_z)
193 else:
194 gamma_array[i] = self.find_growth_rate_resistive(k_z)
195
196 return gamma_array
197
198 def critical_wavenumber(self):
199 """
200 Compute critical wavenumber below which MRI is stable.
201
202 Returns:
203 float: k_crit
204 """
205 # For MRI, need B_z > 0 and k_z > 0
206 # Critical condition from dispersion relation
207 v_A = self.v_A
208 Omega = self.Omega
209
210 if v_A == 0:
211 return np.inf
212
213 # Rough estimate
214 k_crit = Omega / v_A
215
216 return k_crit
217
218 def plot_growth_rate_vs_k(self):
219 """
220 Plot growth rate spectrum.
221 """
222 k_z_array = np.logspace(-2, 2, 200) * self.Omega / self.v_A
223 gamma_array = self.compute_growth_rate_spectrum(k_z_array)
224
225 fig, ax = plt.subplots(figsize=(10, 6))
226
227 ax.plot(k_z_array * self.v_A / self.Omega, gamma_array / self.Omega,
228 'b-', linewidth=2, label='MRI growth rate')
229
230 # Maximum growth rate
231 gamma_max = np.max(gamma_array)
232 k_max = k_z_array[np.argmax(gamma_array)]
233
234 ax.axhline(y=gamma_max / self.Omega, color='r', linestyle='--',
235 alpha=0.5, label=f'γ_max/Ω = {gamma_max/self.Omega:.3f}')
236 ax.axvline(x=k_max * self.v_A / self.Omega, color='g', linestyle='--',
237 alpha=0.5, label=f'k_max v_A/Ω = {k_max*self.v_A/self.Omega:.2f}')
238
239 ax.set_xlabel(r'$k_z v_A / \Omega$', fontsize=12)
240 ax.set_ylabel(r'$\gamma / \Omega$', fontsize=12)
241 ax.set_title(f'MRI Growth Rate Spectrum (q={self.q}, Pm={self.Pm:.1e})',
242 fontsize=14)
243 ax.set_xscale('log')
244 ax.grid(True, alpha=0.3)
245 ax.legend(fontsize=11)
246
247 # Add theoretical max for ideal MRI
248 if self.eta == 0:
249 ax.axhline(y=0.75, color='orange', linestyle=':', linewidth=2,
250 alpha=0.6, label='Ideal limit: 3Ω/4')
251 ax.legend(fontsize=11)
252
253 plt.tight_layout()
254 return fig
255
256 def plot_stability_boundary(self):
257 """
258 Plot stability boundary in (k, Pm) space.
259 """
260 fig, ax = plt.subplots(figsize=(10, 7))
261
262 # Array of Pm values
263 Pm_array = np.logspace(-2, 2, 50)
264 k_array = np.logspace(-2, 2, 100) * self.Omega / self.v_A
265
266 # Compute growth rate for each (k, Pm)
267 gamma_grid = np.zeros((len(Pm_array), len(k_array)))
268
269 for i, Pm in enumerate(Pm_array):
270 # Set resistivity
271 eta_temp = self.nu / Pm if Pm > 0 else 0
272 eta_original = self.eta
273 self.eta = eta_temp
274
275 for j, k_z in enumerate(k_array):
276 gamma_grid[i, j] = self.find_growth_rate_resistive(k_z)
277
278 # Restore
279 self.eta = eta_original
280
281 # Normalize
282 gamma_grid /= self.Omega
283
284 # Contour plot
285 K, P = np.meshgrid(k_array * self.v_A / self.Omega, Pm_array)
286 levels = np.linspace(0, 0.75, 16)
287 contour = ax.contourf(K, P, gamma_grid, levels=levels, cmap='viridis')
288 plt.colorbar(contour, ax=ax, label=r'$\gamma / \Omega$')
289
290 # Stability boundary (γ = 0)
291 ax.contour(K, P, gamma_grid, levels=[0], colors='red',
292 linewidths=2, linestyles='--')
293
294 ax.set_xlabel(r'$k_z v_A / \Omega$', fontsize=12)
295 ax.set_ylabel(r'Magnetic Prandtl Number $Pm$', fontsize=12)
296 ax.set_title('MRI Stability Boundary', fontsize=14)
297 ax.set_xscale('log')
298 ax.set_yscale('log')
299 ax.grid(True, alpha=0.3)
300
301 plt.tight_layout()
302 return fig
303
304
305def main():
306 """
307 Main function demonstrating MRI analysis.
308 """
309 print("=" * 60)
310 print("Magnetorotational Instability (MRI) Analysis")
311 print("=" * 60)
312
313 # Create MRI model - ideal case
314 print("\n1. Ideal MHD (no resistivity):")
315 mri_ideal = MRIShearingBox(
316 Omega=1.0,
317 q=1.5, # Keplerian
318 B_z=0.1,
319 rho=1.0,
320 eta=0.0,
321 nu=0.0
322 )
323
324 print(f" Alfvén velocity: {mri_ideal.v_A:.3f}")
325 print(f" Epicyclic frequency: κ = {np.sqrt(mri_ideal.kappa2):.3f}")
326
327 # Resistive case
328 print("\n2. Resistive MHD:")
329 mri_resistive = MRIShearingBox(
330 Omega=1.0,
331 q=1.5,
332 B_z=0.1,
333 rho=1.0,
334 eta=0.01,
335 nu=0.01
336 )
337 print(f" Magnetic Prandtl number: Pm = {mri_resistive.Pm:.2f}")
338
339 # Plot growth rate spectrum
340 fig1 = mri_ideal.plot_growth_rate_vs_k()
341
342 # Plot stability boundary
343 fig2 = mri_ideal.plot_stability_boundary()
344
345 plt.savefig('/tmp/mri_shearing_box.png', dpi=150, bbox_inches='tight')
346 print("\nPlots saved to /tmp/mri_shearing_box.png")
347
348 plt.show()
349
350
351if __name__ == "__main__":
352 main()