12. 강착 원반 MHD¶

학습 목표¶

이 강의를 마치면 다음을 할 수 있어야 합니다:

강착 원반에서 각운동량 수송의 근본적인 문제 이해
자기 회전 불안정성(MRI) 분산 관계 유도 및 분석
천체물리 원반 강착에 MRI가 필수적인 이유 설명
MRI 성장률 및 특성 파장 계산
맥스웰 및 레이놀즈 응력을 통한 각운동량 수송 이해
MRI 난류를 α-원반 모델에 연결
원반 바람 및 제트 형성 메커니즘 설명 (Blandford-Payne, 자기 타워)
MRI 및 원반 물리학의 수치 모델 구현

1. 강착 원반 기초¶

1.1 각운동량 문제¶

강착 원반(Accretion disk)은 낙하하는 물질이 상당한 각운동량을 가질 때 컴팩트 천체 (블랙홀, 중성자별, 백색왜성) 및 원시별 주위에 형성됩니다. 물질은 회전하는 원반 구성으로 정착합니다.

케플러 회전:

중심 질량 M 주위의 반경 r에서 원 궤도를 도는 시험 입자의 경우:

원심력 = 중력
v²/r = GM/r²
v = √(GM/r)

각속도:

Ω(r) = v/r = √(GM/r³) ∝ r^{-3/2}

이것이 케플러 회전(Keplerian rotation)입니다: 각속도가 반경에 따라 감소합니다.

문제:

물질이 강착(안쪽으로 이동)하려면 각운동량을 잃어야 합니다. 하지만 어떻게?

비각운동량(Specific angular momentum):

ℓ = r² Ω = √(GM r)

슈바르츠실트 블랙홀의 최내곽 안정 원궤도(ISCO)에서:

r_ISCO = 6 GM/c²
ℓ_ISCO = √(6 GM² / c²)

큰 r의 물질은 훨씬 큰 ℓ를 가집니다. ISCO에 도달하려면 각운동량을 벗어야 합니다.

가능한 메커니즘:

분자 점성: 너무 작음 (천체물리 원반에서 Re ~ 10¹⁴)
중력 토크: 쌍성계 또는 자가 중력 원반에서 (한계적)
자기장 + 난류: 이것이 핵심입니다!

1.2 α-원반 모델 (Shakura & Sunyaev 1973)¶

각운동량 수송에 대한 상세한 이해가 없는 상태에서, Shakura & Sunyaev는 이것을 매개변수화했습니다:

점성 응력:

τ_rφ = -ρ ν_eff r dΩ/dr

여기서 ν_eff는 유효 점성(Effective viscosity)입니다 (분자가 아닙니다!).

매개변수화:

ν_eff = α c_s H

여기서: - α는 무차원 매개변수 (0 < α < 1) - c_s는 음속 - H는 원반 스케일 높이

물리적 해석:

크기 ~ H의 난류 와류가 속도 ~ c_s로 움직임
혼합 길이 이론: ν_eff ~ v_turb × ℓ_mix ~ c_s H
α는 각운동량 수송의 효율을 측정

전형적인 값:

관측으로부터 (X선 쌍성계, AGN에 맞춤):

α ~ 0.01 - 0.1

핵심 질문: 어떤 물리적 과정이 α를 설정합니까? 수십 년 동안 이것은 알려지지 않았습니다. 답: MRI 구동 난류.

1.3 원반 구조 방정식¶

질량 보존:

∂Σ/∂t + (1/r) ∂(r Σ v_r)/∂r = 0

여기서 Σ는 표면 밀도, v_r은 반경 속도입니다.

각운동량 보존:

∂(Σ r² Ω)/∂t + (1/r) ∂(r² Σ v_r r² Ω)/∂r = (1/r) ∂(r² τ_rφ)/∂r

에너지 방정식:

점성 소산이 원반을 가열합니다:

Q_vis = τ_rφ r dΩ/dr

이 에너지는 복사로 방출됩니다:

Q_rad = σ T_eff⁴

여기서 T_eff는 유효 표면 온도입니다.

정상 상태 강착:

정상 상태에서 ∂/∂t = 0이고 일정한 강착률 Ṁ이 있습니다:

Σ v_r × 2πr = -Ṁ

(유입이므로 음수.)

온도 프로파일:

기하학적으로 얇은 원반 (H ≪ r)의 경우:

T_eff ∝ (Ṁ / r³)^{1/4}

블랙홀 강착 원반의 경우:

T_eff ~ 10⁶ (Ṁ / Ṁ_Edd)^{1/4} (M / 10 M_☉)^{-1/4} (r / 10 R_s)^{-3/4} K

이것은 내부 영역에서 X선을 방출하기에 충분히 뜨겁습니다!

2. 자기 회전 불안정성 (MRI)¶

2.1 발견 및 중요성¶

Balbus & Hawley (1991)는 차등 회전 원반의 약한 자기장이 자기 회전 불안정성(Magnetorotational Instability, MRI)에 대해 선형 불안정하다는 것을 재발견했습니다 (Velikhov 1959, Chandrasekhar 1960 이후).

중요성:

MRI는 강착 원반 이론에서 가장 중요한 불안정성
난류 생성 → 유효 점성 → 강착
MRI 난류 원반에서 α ~ 0.01 설정
원시행성 원반, X선 쌍성계, AGN, 조석 파괴 사건에서 작동

2.2 물리적 메커니즘: 스프링 유추¶

설정:

반경 r 및 r + δr에 있는 두 유체 요소
자기장선으로 연결됨 (스프링처럼 작용)
원반이 케플러 회전: Ω(r) ∝ r^{-3/2}

교란되지 않은 상태:

두 요소 모두 국소 Ω(r)에서 회전.

교란:

내부 요소를 바깥쪽으로, 외부 요소를 안쪽으로 변위 (약간의 반경 교란).

진화:

자기장 없이:
내부 요소가 바깥쪽으로 이동 → 각운동량 보존 → 국소 케플러보다 느리게 회전 → 뒤처짐
외부 요소가 안쪽으로 이동 → 더 빠르게 회전 → 앞서 나감
방위각으로 떨어져 표류 → 안정 (Rayleigh 기준)
자기장과 함께:
자기 장력이 두 요소를 연결하는 스프링처럼 작용
내부 요소 (이제 바깥쪽, 느림)가 장에 의해 앞으로 당겨짐 → 각운동량 획득 → 더 바깥쪽으로 이동
외부 요소 (이제 안쪽, 빠름)가 뒤로 당겨짐 → 각운동량 손실 → 더 안쪽으로 이동
양의 피드백 → 불안정성!

핵심 통찰:

자기장이 각운동량을 바깥쪽으로 수송하여 (내부에서 외부 요소로), 내부 요소가 바깥쪽으로, 외부가 안쪽으로 이동하도록 허용 → Rayleigh 안정성에도 불구하고 불안정성.

2.3 국소 선형 분석: 전단 박스¶

MRI를 분석하기 위해 전단 박스(Shearing box) 근사를 사용합니다 (Goldreich & Lynden-Bell 1965):

좌표:

반경 r_0에서 국소 데카르트 프레임
x = r - r_0 (반경, 바깥쪽)
y = r_0 (φ - Ω_0 t) (방위각, 회전 프레임에서)
z (수직)

전단 흐름:

회전 프레임에서 배경 속도는:

v_y = -q Ω_0 x

여기서 q = -d ln Ω / d ln r입니다. 케플러의 경우: q = 3/2.

선형화된 MHD 방정식:

균일한 수직 장 B_0 = B_z ẑ 주위에서 교란:

∂v'/∂t + (v·∇)v' = -(1/ρ)∇p' + (1/ρμ_0)(∇×B')×B_0 + 2q Ω_0² x x̂
∂B'/∂t = ∇×(v'×B_0)
∇·v' = 0, ∇·B' = 0

2q Ω_0² x x̂ 항은 조석력(Tidal force)입니다 (중력의 국소 전개로부터).

평면파 해 찾기:

v' ~ exp(i k·x + γ t)
B' ~ exp(i k·x + γ t)

2.4 MRI 분산 관계¶

단순화를 위해, 수직 장 B_0 = B_z ẑ 및 파수 k = k_z를 가진 축대칭 모드 (k_y = 0 처음)를 고려합니다.

비압축성:

∇·v' = 0 가정 (아음속 교란에 유효).

분산 관계:

약간의 대수학 후 (유도는 Balbus & Hawley 1991, 1998 참조), 성장률 γ는 다음을 만족합니다:

γ⁴ + γ² [2κ² - (2 - q) Ω_0²] + κ² [κ² - (k v_A)²] = 0

여기서: - κ² = 2 Ω_0 (Ω_0 + r dΩ/dr) = (2 - q) Ω_0²는 주전원 주파수(Epicyclic frequency) 제곱 - v_A = B_0 / √(μ_0 ρ)는 알벤 속도

케플러 회전 (q = 3/2)의 경우:

κ² = Ω_0²

분산 관계는 다음이 됩니다:

γ⁴ + γ² (2 Ω_0² - Ω_0²/2) + Ω_0² [Ω_0² - (k v_A)²] = 0
γ⁴ + (3/2) Ω_0² γ² + Ω_0² [Ω_0² - (k v_A)²] = 0

불안정성 조건:

γ²가 양수 (지수 성장)이려면:

(γ²에 대한) 이차 방정식의 판별식이 양수이고 적어도 하나의 근 γ² > 0이어야 합니다.

풀면:

γ² = -(3/4) Ω_0² ± √[(9/16) Ω_0⁴ - Ω_0² (Ω_0² - (k v_A)²)]
    = -(3/4) Ω_0² ± Ω_0² √[(9/16) - 1 + (k v_A / Ω_0)²]
    = -(3/4) Ω_0² ± Ω_0² √[(k v_A / Ω_0)² - 7/16]

불안정성을 위해:

(k v_A / Ω_0)² < 7/16  (긴 파장)

그러면 + 근이 어떤 k v_A / Ω_0 범위에서 γ² > 0을 줍니다.

표준 결과 (Balbus & Hawley 1998):

약한 수직 장을 가진 케플러 원반의 경우:

γ² = (1/2) [(k v_A)² - Ω_0²] + (1/2) √[(k v_A)² + Ω_0²)² - 16 q Ω_0² (k v_A)²]

q = 3/2 (케플러)의 경우:

최대 성장률:

γ_max = (√(3)/4) Ω_0 ≈ 0.433 Ω_0

이것은 k v_A → 0 (긴 파장)에서 발생합니다.

핵심 결과:

불안정성 기준:
케플러 회전의 경우: 모든 자기장 (아무리 약해도!)이 불안정
조건: dΩ²/d ln r < 0 (바깥쪽으로 각속도 감소)
성장률:
궤도 주파수와 비슷: γ ~ Ω_0 (매우 빠름!)
E-폴딩 시간: τ = 1/γ ~ 1/(Ω_0) ~ 궤도 주기
가장 빠르게 성장하는 모드:
파장: λ_MRI ~ 2π v_A / Ω_0
약한 장의 경우, λ_MRI는 원반 높이보다 훨씬 작을 수 있음

2.5 비축대칭 모드¶

k_y ≠ 0 (방위각 구조)를 가진 모드의 경우, 불안정성이 지속됩니다. 반경 및 방위각 파수를 포함하면:

전체 분산 관계 (일반):

γ⁴ + γ² [κ² + (k·v_A)² - 2 Ω_0 k_y v_{Ay}]
    + κ² [(k·v_A)² - 4 Ω_0 k_y v_{Ay}] = 0

여기서 v_A = B_0 / √(μ_0 ρ)는 알벤 속도 벡터입니다.

핵심 포인트:

약한 장의 모든 방향 (수직, toroidal 등)에 대해 불안정성 존재
성장률은 케플러 원반에서 항상 ~ Ω_0
MRI는 자화된 원반에서 편재

2.6 케플러 원반이 Rayleigh 안정하지만 MRI 불안정한 이유¶

Rayleigh 기준:

회전 유체가 유체역학적 (비자기) 교란에 안정하려면:

d(r² Ω)² / dr > 0

케플러의 경우: Ω ∝ r^{-3/2} → r² Ω ∝ r^{1/2} → 도함수 > 0 → 안정.

이것이 케플러 원반이 전단만으로는 유체역학적 난류를 가질 수 없는 이유입니다.

MRI가 게임을 바꿉니다:

자기장은 Rayleigh 기준을 우회하는 각운동량 수송 채널을 제공합니다. 자기 장력은 안정화 원심 효과를 극복하는 음의 유효 점성 또는 불안정화 토크를 생성합니다.

3. MRI 난류에서의 각운동량 수송¶

3.1 맥스웰 응력¶

MHD에서 응력 텐서는 유체 운동 (레이놀즈 응력) 및 자기장 (맥스웰 응력) 모두로부터 기여를 받습니다.

총 응력 텐서:

T_{ij} = ρ v_i v_j + p δ_{ij} + (B_i B_j / μ_0 - B² δ_{ij} / (2μ_0))

각운동량 수송:

각운동량의 반경 수송 (단위 면적당)은:

τ_rφ = ρ v_r v_φ - B_r B_φ / μ_0

여기서: - ρ v_r v_φ는 레이놀즈 응력(Reynolds stress) (유체역학) - -B_r B_φ / μ_0는 맥스웰 응력(Maxwell stress) (자기)

부호 규약:

양의 τ_rφ → 각운동량의 바깥쪽 수송.

3.2 MRI 난류 시뮬레이션¶

전단 박스 MHD 코드를 사용한 수치 시뮬레이션 (Hawley, Balbus, Stone 등)은 다음을 발견합니다:

맥스웰 응력이 지배적:

⟨B_r B_φ⟩ / μ_0  ≫  ⟨ρ v_r v_φ⟩

일반적으로 맥스웰 응력이 레이놀즈 응력보다 ~10배 더 큽니다.

α 매개변수:

정의:

α = ⟨τ_rφ⟩ / ⟨p⟩

여기서 ⟨·⟩는 시간 및 부피 평균을 나타냅니다.

시뮬레이션으로부터 (예: Hawley, Gammie, Balbus 1995):

순 수직 플럭스 없음: α ~ 0.01-0.02
순 수직 플럭스와 함께: α ~ 0.05-0.5 (더 강한 수송)

장 구성이 중요:

순 플럭스 (정렬된 장): 지속된 채널 모드, 높은 α
순 플럭스 없음 (난류 장): 지속된 난류이지만 낮은 α

3.3 유효 점성¶

α 매개변수로부터:

ν_eff = α c_s H

c_s ~ 10 km/s, H ~ 0.1 R ~ 10⁹ cm, α ~ 0.01인 전형적인 원반의 경우:

ν_eff ~ 10¹⁵ cm²/s

이것은 분자 점성 (ν_mol ~ 1 cm²/s)에 비해 엄청나지만, MRI 난류로부터 자연스럽게 발생합니다.

강착 시간 척도:

τ_acc ~ R² / ν_eff ~ (10¹⁰ cm)² / (10¹⁵ cm²/s) ~ 10⁵ s ~ 1일

이것은 컴팩트 천체로의 급속한 강착을 허용하여 관측된 X선 변동성을 설명합니다.

4. 비선형 MRI 및 난류 포화¶

4.1 채널 해¶

채널 모드는 전단 박스에서 비압축성 MHD 방정식의 정확한 비선형 해입니다 (Goodman & Xu 1994).

구조:

방위각 방향의 이동 파동
선형 영역에서 지수적으로 성장
비선형 영역에서 지속됨 (순 수직 플럭스를 가진 원반의 경우)

에너지:

채널 모드는 포화 상태에서 대부분의 자기 에너지와 응력을 운반합니다 (순 플럭스 경우).

기생 불안정성:

채널 모드 자체가 그것들을 분해하는 기생 모드(Parasitic modes) (2차 불안정성)에 불안정하여 난류로 이어집니다.

4.2 포화 메커니즘¶

로렌츠 힘 역반응:

MRI가 성장함에 따라 자기장 강도가 증가합니다. 로렌츠 힘이 흐름을 변경합니다:

J × B ~ B² / L

자기 압력이 열 압력과 비슷해질 때:

B² / (2μ_0) ~ p

또는

β = 2μ_0 p / B² ~ 1

성장이 포화됩니다.

전형적인 포화:

시뮬레이션은 다음을 발견합니다:

⟨B²⟩ / (2μ_0 ⟨p⟩) ~ 1-10  (단위 차수)

자기 에너지는 장 구성에 따라 아열(Sub-thermal)에서 초열(Super-thermal)입니다.

4.3 층화된 원반에서의 MRI¶

(z 방향 중력을 가진) 층화된 전단 박스 시뮬레이션에서:

나비 효과:

자기장이 나비 패턴을 발전시킵니다: 중간면 위와 아래에서 교대 toroidal 장 극성, z에서 전파.

대역 흐름(Zonal flows):

레이놀즈 응력으로 인해 대규모 방위각 흐름 (φ와 독립)이 발생합니다.

강착 응력:

수직 범위에 걸쳐 평균한 유효 α는 비층화 경우와 유사하지만 (~0.01-0.1), 상당한 수직 변화가 있습니다.

5. 사영역 및 비이상 MHD 효과¶

5.1 원시행성 원반에서의 MRI¶

문제:

MRI는 이온화 (기체와 자기장 사이의 결합)를 필요로 합니다. 원시행성 원반에서:

내부 영역 (< 0.1 AU): 뜨겁고, 열적으로 이온화됨 → MRI 활성
외부 영역 (> 1 AU) 중간면: 차갑고, 밀도 높음 → 불량하게 이온화됨 → MRI 억제됨

사영역(Dead zone):

MRI에 대해 이온화가 너무 낮은 영역. 기준:

자기 레이놀즈 수: Rm = v L / η_Ohm > 10⁴ (MRI의 경우)

여기서 η_Ohm = c² / (4π σ)는 옴 저항입니다.

젊은 별 주위 1 AU의 원반 중간면에서: - 온도: T ~ 100-300 K - 이온화 분율: x_e ~ 10^{-13} (우주선, 방사성 붕괴로부터) - Rm ~ 10-100 → MRI 억제됨

층상 강착(Layered accretion):

활성 영역: 표면 근처 (UV, X선, 우주선에 의해 이온화됨)
사영역: 중간면 (중성, MRI 없음)
강착은 주로 활성층에서 진행됨

5.2 비이상 MHD 효과¶

약하게 이온화된 플라즈마에서 세 가지 비이상 효과가 중요해집니다:

1. 옴 확산(Ohmic diffusion):

∂B/∂t = ∇×(v×B) + η_Ohm ∇²B

저항 η_Ohm이 소규모 장 변동을 감쇠시킵니다.

2. 양극성 확산(Ambipolar diffusion):

이온과 중성자가 완벽하게 결합되지 않습니다. 자기장이 중성자를 통해 미끄러집니다:

∂B/∂t = ∇×(v×B) + ∇×(η_AD (∇×B)×B / B²)

여기서 η_AD ~ B / (ρ_n γ_in), γ_in은 이온-중성자 충돌률입니다.

3. 홀 효과(Hall effect):

전류 존재 하에서 전자가 이온에 대해 표류합니다:

∂B/∂t = ∇×(v×B) - ∇×(η_H (∇×B)×B / |B|)

여기서 η_H ~ B / (e n_e).

MRI에 대한 영향:

옴 확산: 작은 척도에서 MRI 억제 → 임계 파장 증가
양극성 확산: 밀도 높고 약하게 이온화된 영역에서 MRI 억제
홀 효과: MRI를 수정 가능 (장 방향 및 홀 항의 부호에 따라 향상 또는 억제)

결과:

원시행성 원반에서 MRI는 다음에서만 활성일 수 있습니다: - 표면층 - 내부 뜨거운 영역 - 향상된 이온화 영역 (별 근처, 바람 쓸린 영역)

대안 메커니즘 (예: 수직 전단 불안정성, 중력 불안정성)이 사영역에서 작동할 수 있습니다.

6. 원반 바람 및 제트¶

6.1 Blandford-Payne 메커니즘 (1982)¶

원반을 관통하는 자기장이 원심 구동 바람(Centrifugally driven wind)을 발사할 수 있습니다.

설정:

원반을 관통하는 poloidal 자기장 B_p
원반에 고정된 장선, 바깥쪽으로 굽힘
기체가 장선을 따라 흐름

원심 발사 조건:

장선이 원반 법선 (수직)으로부터 각도 θ로 기울어져야 합니다:

θ > 30°  (임계 각도)

물리적 그림:

반경 r의 기체가 각속도 Ω(r)을 가짐
바깥쪽으로 기울어진 장선을 따라 이동하면 각운동량을 보존: L = r v_φ = 상수
더 큰 원통 반경 R > r에서 원심력 L²/(R³)이 중력 + 자기 장력을 초과 가능
기체가 바깥쪽으로 던져짐 (원심 가속)

가속:

바람이 다음으로 가속됩니다:

v_∞² ~ 2 G M / r_0

여기서 r_0는 발사 반경입니다. 이것은 탈출 속도와 비슷합니다.

질량 손실률:

바람으로 방출되는 강착 흐름의 분율:

Ṁ_wind / Ṁ_acc ~ (B_p / √(4π ρ v_K²)) ~ β_p^{-1/2}

여기서 β_p = 8π p / B_p²는 플라즈마 베타입니다.

6.2 자기 타워¶

강한 toroidal 장 (B_φ)의 경우, 자기 압력 기울기가 집속된 유출을 구동할 수 있습니다.

메커니즘:

차등 회전이 poloidal 장을 감음 → 강한 toroidal 장 B_φ
Toroidal 장이 압력 B_φ² / (2μ_0)을 가짐
z (수직)의 압력 기울기 → 위쪽 힘
기체가 축을 따라 위로 밀림

후프 응력(Hoop stress):

Toroidal 장도 핀치 힘 (후프 응력)을 가합니다:

F_pinch ~ -∂(B_φ² / (2μ_0)) / ∂R

이것이 흐름을 축을 향해 집속 → 제트 형성.

시뮬레이션:

MHD 시뮬레이션 (예: Lynden-Bell, Kato, Kudoh)은 다음의 조합이: - Toroidal 장 축적 - 자기 타워 형성 - 후프 응력에 의한 집속

젊은 항성 천체 및 AGN에서 관측된 것과 유사한 쌍극 제트를 생성할 수 있음을 보여줍니다.

6.3 AGN 제트에의 응용¶

활동 은하핵(Active Galactic Nuclei, AGN):

은하 중심의 초대질량 블랙홀 (10⁶-10⁹ M☉)이 기체를 강착시키고 강력한 상대론적 제트를 발사합니다: - 길이: Mpc 척도까지 - 속도: v ~ 0.1-0.99 c - 전력: L_jet ~ 10⁴²-10⁴⁷ erg/s

자기 발사:

주요 모델: Blandford-Znajek 메커니즘 (1977) - 블랙홀 지평선을 관통하는 자기장 (원반만이 아님) - 블랙홀 회전이 장선을 비틂 → 전자기 에너지 추출 - 전력: P ~ B² a² c 여기서 a는 블랙홀 스핀

대안으로, 원반으로부터의 Blandford-Payne이 기여할 수 있습니다.

집속:

Toroidal 장 (회전으로부터) + 외부 압력 (원반 바람, 고치) → 집속된 제트.

관측 증거:

VLBI 영상: ~10 슈바르츠실트 반경까지 해상된 제트
Event Horizon Telescope (EHT): 지평선 척도의 M87* 제트 기저
편광: 정렬된 poloidal+toroidal 장과 일치

7. Python 구현¶

7.1 MRI 분산 관계¶

import numpy as np
import matplotlib.pyplot as plt

def mri_dispersion():
    """
    케플러 원반에 대한 MRI 분산 관계를 풉니다.

    분산 관계 (수직 장, 비압축성):
      γ⁴ + γ² [κ² + (k v_A)²] - 3 Ω² (k v_A)² = 0

    케플러의 경우: κ² = Ω²
    """
    Omega = 1.0  # 궤도 주파수 (정규화)

    # k v_A / Omega 범위
    k_vA_over_Omega = np.linspace(0, 2, 200)

    # 분산 관계: γ⁴ + γ² [Ω² + (k v_A)²] - 3 Ω² (k v_A)² = 0
    # X = γ² / Ω², Y = (k v_A / Ω)²로 두면
    # X² + X [1 + Y] - 3 Y = 0

    Y = k_vA_over_Omega**2

    # X = γ² / Ω²에 대한 이차 방정식 풀기
    # X = -(1+Y)/2 ± √[(1+Y)²/4 + 3Y]

    discriminant = (1 + Y)**2 / 4 + 3 * Y
    X_plus = -(1 + Y)/2 + np.sqrt(discriminant)
    X_minus = -(1 + Y)/2 - np.sqrt(discriminant)

    # 성장률 γ / Ω (양의 근 선택)
    gamma_over_Omega = np.where(X_plus > 0, np.sqrt(X_plus), 0)

    # 최대 성장률 (k v_A → 0에서)
    gamma_max = np.sqrt(X_plus[0]) * Omega
    print(f"최대 성장률: γ_max / Ω = {gamma_max:.4f}")
    print(f"             γ_max = {gamma_max:.4f} Ω")

    # 플롯
    plt.figure(figsize=(10, 6))
    plt.plot(k_vA_over_Omega, gamma_over_Omega, 'b-', linewidth=2.5)
    plt.axhline(gamma_max, color='r', linestyle='--', linewidth=1.5, label=f'최대: γ/Ω = {gamma_max:.3f}')
    plt.axvline(0, color='k', linestyle='-', linewidth=0.5)
    plt.axhline(0, color='k', linestyle='-', linewidth=0.5)

    plt.xlabel('$k v_A / \\Omega$', fontsize=14)
    plt.ylabel('$\\gamma / \\Omega$', fontsize=14)
    plt.title('MRI 분산 관계 (케플러 원반, 수직 장)', fontsize=16)
    plt.legend(fontsize=12)
    plt.grid(True, alpha=0.3)
    plt.xlim(0, 2)
    plt.ylim(0, 0.5)
    plt.savefig('mri_dispersion.png', dpi=150)
    plt.show()

    # 가장 불안정한 파장
    idx_max = np.argmax(gamma_over_Omega)
    k_vA_max = k_vA_over_Omega[idx_max] * Omega
    print(f"\n가장 불안정한 모드: k v_A = {k_vA_max:.4f} Ω")
    print(f"파장: λ = 2π / k = {2*np.pi / k_vA_max:.2f} (v_A / Ω)")

mri_dispersion()

7.2 장 강도 대 MRI 성장률¶

import numpy as np
import matplotlib.pyplot as plt

def mri_growth_vs_field():
    """
    자기장 강도의 함수로 MRI 성장률을 플롯합니다.
    """
    Omega = 1.0  # 궤도 주파수
    rho = 1.0    # 밀도 (정규화)
    mu_0 = 1.0   # 투자율 (정규화)

    # 자기장 강도 범위
    B = np.logspace(-3, 1, 100)

    # 알벤 속도
    v_A = B / np.sqrt(mu_0 * rho)

    # 최대 성장 (k v_A → 0)의 경우, γ_max ≈ 0.75 Ω (케플러)
    # 하지만 성장률은 k v_A / Ω에 의존

    # 고정 파장 선택: k = Ω / (H), 여기서 H ~ c_s / Ω
    # 그러면 k v_A / Ω = v_A / c_s
    c_s = 1.0  # 음속 (정규화)

    k_vA_over_Omega = v_A / c_s

    # 분산 관계로부터 성장률 계산
    Y = k_vA_over_Omega**2
    discriminant = (1 + Y)**2 / 4 + 3 * Y
    X_plus = -(1 + Y)/2 + np.sqrt(discriminant)
    gamma_over_Omega = np.where(X_plus > 0, np.sqrt(X_plus), 0)

    # 플롯
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

    # B 대 성장률
    ax1.semilogx(B, gamma_over_Omega, 'b-', linewidth=2.5)
    ax1.axhline(0.75, color='r', linestyle='--', linewidth=1, label='최대 (k→0): γ/Ω ≈ 0.75')
    ax1.set_xlabel('자기장 $B$ (정규화)', fontsize=14)
    ax1.set_ylabel('성장률 $\\gamma / \\Omega$', fontsize=14)
    ax1.set_title('자기장 강도 대 MRI 성장률', fontsize=16)
    ax1.legend(fontsize=12)
    ax1.grid(True, alpha=0.3)

    # 가장 빠르게 성장하는 모드의 파장
    lambda_MRI = 2 * np.pi * v_A / (Omega * gamma_over_Omega)
    lambda_MRI = np.where(gamma_over_Omega > 0, lambda_MRI, np.nan)

    ax2.loglog(B, lambda_MRI, 'b-', linewidth=2.5)
    ax2.set_xlabel('자기장 $B$ (정규화)', fontsize=14)
    ax2.set_ylabel('MRI 파장 $\\lambda_{MRI} / H$', fontsize=14)
    ax2.set_title('가장 빠르게 성장하는 MRI 파장', fontsize=16)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('mri_growth_vs_field.png', dpi=150)
    plt.show()

mri_growth_vs_field()

7.3 전단 박스 국소 분석¶

import numpy as np
import matplotlib.pyplot as plt

def shearing_box_trajectories():
    """
    전단 박스에서 유체 요소 궤적을 시각화합니다.

    전단 흐름: v_y = -q Ω x  (케플러: q = 3/2)
    """
    Omega = 1.0
    q = 1.5  # 케플러

    # 초기 위치
    N_particles = 20
    x0 = np.random.uniform(-2, 2, N_particles)
    y0 = np.random.uniform(-2, 2, N_particles)

    # 시간 진화
    t_max = 10.0
    dt = 0.1
    Nt = int(t_max / dt)

    # 궤적
    fig, ax = plt.subplots(figsize=(10, 10))

    for i in range(N_particles):
        x = np.zeros(Nt)
        y = np.zeros(Nt)

        x[0] = x0[i]
        y[0] = y0[i]

        for n in range(Nt - 1):
            # 전단 흐름
            v_y = -q * Omega * x[n]

            # 업데이트 (Euler)
            x[n+1] = x[n]
            y[n+1] = y[n] + dt * v_y

        ax.plot(x, y, alpha=0.7, linewidth=1.5)
        ax.plot(x[0], y[0], 'go', markersize=5)
        ax.plot(x[-1], y[-1], 'ro', markersize=5)

    ax.set_xlabel('반경 $x$', fontsize=14)
    ax.set_ylabel('방위각 $y$', fontsize=14)
    ax.set_title('전단 박스에서 유체 요소 궤적 (케플러)', fontsize=16)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    ax.legend(['시작', '끝'], fontsize=12)
    plt.savefig('shearing_box_trajectories.png', dpi=150)
    plt.show()

shearing_box_trajectories()

7.4 Blandford-Payne 바람 해 (단순화)¶

import numpy as np
import matplotlib.pyplot as plt

def blandford_payne_wind():
    """
    Blandford-Payne 원심 바람의 단순화된 모델.

    가정:
      - 장선 형태: z = r tan(θ)
      - 각운동량 보존: r v_φ = r₀² Ω₀
      - 원심 대 중력 균형
    """
    # 매개변수
    GM = 1.0       # 중력 매개변수 (정규화)
    r_0 = 1.0      # 발사 반경
    Omega_0 = np.sqrt(GM / r_0**3)  # r_0에서의 케플러 각속도

    # 장선 기울기
    theta_deg = np.array([30, 45, 60, 75])  # 도
    theta = np.radians(theta_deg)

    # 장선을 따른 원통 반경
    R = np.linspace(r_0, 5*r_0, 100)

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

    for th, th_d in zip(theta, theta_deg):
        # 높이 z = (R - r_0) tan(θ)
        z = (R - r_0) * np.tan(th)

        # 구 반경
        r_sph = np.sqrt(R**2 + z**2)

        # 각운동량 (보존)
        L = r_0**2 * Omega_0

        # 방위각 속도
        v_phi = L / R

        # 원심력 (단위 질량당)
        F_cent = v_phi**2 / R

        # 중력 (반경 성분)
        F_grav = GM / r_sph**2

        # 유효 퍼텐셜
        Phi_eff = -GM / r_sph + L**2 / (2 * R**2)

        # 장선 플롯
        ax1.plot(R, z, label=f'θ = {th_d}°', linewidth=2)

        # 유효 퍼텐셜 플롯
        ax2.plot(r_sph, Phi_eff, label=f'θ = {th_d}°', linewidth=2)

    ax1.set_xlabel('원통 반경 $R$ ($r_0$ 단위)', fontsize=14)
    ax1.set_ylabel('높이 $z$ ($r_0$ 단위)', fontsize=14)
    ax1.set_title('Blandford-Payne 장선 기하학', fontsize=16)
    ax1.legend(fontsize=12)
    ax1.grid(True, alpha=0.3)

    ax2.set_xlabel('구 반경 $r$ ($r_0$ 단위)', fontsize=14)
    ax2.set_ylabel('유효 퍼텐셜 $\\Phi_{eff}$', fontsize=14)
    ax2.set_title('바람 가속을 위한 유효 퍼텐셜', fontsize=16)
    ax2.axhline(0, color='k', linestyle='--', linewidth=0.5)
    ax2.legend(fontsize=12)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('blandford_payne_wind.png', dpi=150)
    plt.show()

    print("원심 발사 임계 각도: θ > 30°")
    print("θ < 30°의 경우: 유효 퍼텐셜에 장벽이 없음 → 바람이 원심적으로 구동될 수 없음")

blandford_payne_wind()

7.5 MRI 난류 에너지 진화 (토이 모델)¶

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def mri_turbulence_energy():
    """
    MRI 난류 에너지 진화를 위한 토이 모델.

    방정식:
      dE_K/dt = -α Ω E_K + β E_M - ε_K
      dE_M/dt = γ_MRI E_M - β E_M - ε_M

    여기서:
      E_K: 운동 에너지
      E_M: 자기 에너지
      γ_MRI: MRI 성장률
      β: 에너지 교환 (로렌츠 힘)
      ε: 소산
    """
    # 매개변수
    Omega = 1.0
    gamma_MRI = 0.75 * Omega  # MRI 성장률
    alpha_visc = 0.01         # 유효 점성 매개변수
    beta_exchange = 0.1       # 에너지 교환률
    epsilon_K = 0.05          # 운동 소산
    epsilon_M = 0.05          # 자기 소산

    def dE_dt(E, t):
        E_K, E_M = E

        dE_K_dt = -alpha_visc * Omega * E_K + beta_exchange * E_M - epsilon_K * E_K
        dE_M_dt = gamma_MRI * E_M - beta_exchange * E_M - epsilon_M * E_M

        return [dE_K_dt, dE_M_dt]

    # 초기 조건
    E0 = [1.0, 0.01]  # 초기 운동 에너지, 작은 자기 씨앗

    # 시간
    t = np.linspace(0, 100, 1000)

    # 풀기
    E = odeint(dE_dt, E0, t)

    E_K = E[:, 0]
    E_M = E[:, 1]
    E_total = E_K + E_M

    # 플롯
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

    ax1.plot(t, E_K, 'b-', linewidth=2, label='운동 $E_K$')
    ax1.plot(t, E_M, 'r-', linewidth=2, label='자기 $E_M$')
    ax1.plot(t, E_total, 'k--', linewidth=2, label='총 $E_K + E_M$')
    ax1.set_xlabel('시간 ($\\Omega^{-1}$ 단위)', fontsize=14)
    ax1.set_ylabel('에너지', fontsize=14)
    ax1.set_title('MRI 난류: 에너지 진화', fontsize=16)
    ax1.legend(fontsize=12)
    ax1.grid(True, alpha=0.3)

    # 반로그
    ax2.semilogy(t, E_K, 'b-', linewidth=2, label='운동 $E_K$')
    ax2.semilogy(t, E_M, 'r-', linewidth=2, label='자기 $E_M$')
    ax2.set_xlabel('시간 ($\\Omega^{-1}$ 단위)', fontsize=14)
    ax2.set_ylabel('에너지 (로그 척도)', fontsize=14)
    ax2.set_title('MRI 난류: 에너지 진화 (로그 척도)', fontsize=16)
    ax2.legend(fontsize=12)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('mri_turbulence_energy.png', dpi=150)
    plt.show()

    # 포화 값
    print(f"포화:")
    print(f"  운동 에너지: {E_K[-1]:.4f}")
    print(f"  자기 에너지: {E_M[-1]:.4f}")
    print(f"  비율 E_M / E_K: {E_M[-1] / E_K[-1]:.4f}")

mri_turbulence_energy()

8. 요약¶

강착 원반 MHD는 회전, 중력, 자기장의 상호작용에 의해 지배됩니다:

각운동량 문제:
케플러 원반은 Rayleigh 안정 → 유체역학적 난류 없음
물질이 강착하려면 각운동량을 벗어야 함 → 자기장 필요
자기 회전 불안정성 (MRI):
Balbus & Hawley (1991): 차등 회전 원반의 약한 자기장이 불안정
성장률: γ ~ Ω (매우 빠름!)
가장 빠르게 성장하는 파장: λ_MRI ~ v_A / Ω
불안정성 기준: dΩ²/d ln r < 0 (케플러 만족)
메커니즘: 자기 장력이 스프링처럼 작용하여 각운동량 수송을 바깥쪽으로 가능하게 함
각운동량 수송:
맥스웰 응력: -B_r B_φ / μ_0 (지배적)
레이놀즈 응력: ρ v_r v_φ (부차적)
α-매개변수: MRI 난류로부터 α ~ 0.01-0.1
X선 쌍성계, AGN에서 관측된 강착률 설명
비선형 포화:
MRI는 B² / (2μ_0) ~ p까지 성장
로렌츠 힘이 흐름을 수정 → 난류 포화
채널 모드, 기생 불안정성, 지속된 난류
사영역:
약하게 이온화된 영역 (원시행성 원반)에서 옴, 양극성, 홀 효과가 MRI 억제
층상 강착: 표면 근처 활성, 중간면 사영역
대안 불안정성 작동 가능
원반 바람 및 제트:
Blandford-Payne: 기울어진 자기장선을 따라 발사되는 원심 바람 (θ > 30°)
자기 타워: Toroidal 장 압력이 유출을 구동하고 집속
응용: YSO 제트, AGN 제트, 펄서 바람

MRI는 천체물리학에서 가장 중요한 발견 중 하나로, 강착 원반 난류에 대한 오랫동안 찾던 메커니즘을 제공하고 컴팩트 천체로의 급속한 강착을 가능하게 합니다.

연습 문제¶

MRI 성장률: M = 10 M☉ 블랙홀 주위 반경 r = 10¹⁰ cm의 원반에 대해 궤도 주파수 Ω 및 최대 MRI 성장률 γ_max를 계산하세요.
MRI 파장: v_A = 10 km/s 및 Ω = 10⁻³ s⁻¹에 대해 가장 빠르게 성장하는 MRI 모드의 파장을 추정하세요.
맥스웰 응력: ⟨B_r B_φ⟩ = 10⁻² × ⟨p⟩인 경우, 유효 α 매개변수는 얼마입니까?
강착 시간 척도: α = 0.01, c_s = 10 km/s, H = 10⁹ cm, R = 10¹¹ cm에 대해 강착 시간 척도 τ_acc ~ R² / ν_eff를 추정하세요.
등분배 장: ρ = 10⁻⁹ g/cm³, c_s = 10⁷ cm/s인 원반에서 β = B² / (8π p) = 1에서의 자기장 강도를 계산하세요.
사영역 기준: 이온화 분율 x_e = 10⁻¹³, 온도 T = 200 K, 밀도 ρ = 10⁻⁹ g/cm³인 1 AU의 원시행성 원반에 대해 옴 저항 및 자기 레이놀즈 수를 추정하세요. MRI가 활성입니까?
Blandford-Payne 각도: 왜 θ = 30°가 원심 바람 발사의 임계 각도입니까? (힌트: 장선을 따른 원심력과 중력의 균형을 고려하세요.)
제트 전력: 질량 M = 10⁹ M☉, 강착률 Ṁ = 0.1 Ṁ_Edd, 제트 효율 η_jet = 0.1인 블랙홀에 대해 제트 전력을 erg/s로 추정하세요.
Python 연습: 수직 장에 더해 toroidal 장 B_φ를 포함하도록 MRI 분산 관계 코드를 수정하세요. 성장률이 어떻게 변합니까?
심화: MRI 구동 α-점성을 가진 1D 수직 적분 원반 진화 모델을 구현하세요. 물질의 고리로 시작하여 시간에 따라 확산하고 강착하는 것을 관찰하세요.

이전: 태양 MHD | 다음: 핵융합 MHD