수치 상대론

계산물리 시뮬레이션 수치해석 마지막 수정: 2026-09-01 04:13:41

1. 개요[편집]

수치 상대론 Numerical Relativity
푸는 것	아인슈타인 장방정식의 초기값 문제
분해	3+1 (ADM) — 랩스 α, 시프트 βⁱ, 3-계량 γ_ij, 외재곡률 K_ij
두 갈래	구속조건(타원형 4개) · 발전방정식(쌍곡형 12개)
안정 형식	BSSN · 일반화 조화 · Z4c/CCZ4
게이지 국룰	1+log 슬라이싱 + 감마-드라이버 시프트
돌파	2005 — 이진 블랙홀 병합 완주
출력물	Ψ₄ → 파형 카탈로그 → 검출기 템플릿

방정식은 1915년에 나왔고, 그걸 컴퓨터에 태워 블랙홀 두 개를 한 바퀴 돌리는 데 90년이 걸렸다. 사이 90년의 대부분은 “수식을 어떻게 다시 쓸 것인가”에 쓰였다.

수치 상대론은 아인슈타인 장방정식을 해석적으로 풀 수 없는 상황에서 초기값 문제로 재구성해 컴퓨터로 시간 적분하는 분야다. 대상은 이진 블랙홀 병합, 중성자별 병합, 중력붕괴, 임계현상 같은 강한 장·동역학적 시공간이다. 약한 장은 후-뉴턴 전개로, 정상 상태는 섭동론으로 충분하므로, 수치 상대론이 필요한 곳은 정확히 그 둘이 모두 무너지는 구간이다.

이 분야가 특이한 이유는 방정식이 어렵기 때문이 아니다. 방정식을 그대로 이산화하면 무조건 죽기 때문이다. 물리적으로 동일한 방정식을 대수적으로 다르게 적어 놓으면 어떤 것은 안정하고 어떤 것은 폭발하는데, 이 차이를 만드는 것이 쌍곡성(hyperbolicity) 이라는 순수 편미분방정식론의 개념이다. 수치 상대론의 역사는 사실상 “장방정식의 좋은 형식(formulation)을 찾는 역사”이고, 코드 최적화나 격자 해상도는 그다음 문제였다.

블랙홀 시공간의 기하학적 배경(지평선, 커 해, 좌표특이점)은 블랙홀 문서에 정리돼 있다. 이 문서는 수치해석 쪽 절반 — 형식 선택, 초기 데이터의 타원 문제, 이산화, 파형 추출 — 을 다룬다.

2. 3+1 분해와 두 종류의 방정식[편집]

$G_{\mu\nu}=8\pi T_{\mu\nu}$ 는 4차원 시공간 전체에 대한 관계식이라 “시간을 전진시킨다”는 개념이 없다. 시공간을 공간꼴 초곡면 $\Sigma_t$ 로 얇게 썰고 좌표를 그 위에 얹는 절차가 ADM 3+1 분해(아노윗–데저–미즈너, 1959~62)다.

ds^2 = -\alpha^2 dt^2 + \gamma_{ij}\left(dx^i+\beta^i dt\right)\left(dx^j+\beta^j dt\right)

$\alpha$ — 랩스. 이웃한 두 슬라이스 사이에 고유시간이 얼마나 흐르는가.
$\beta^i$ — 시프트. 같은 좌표값을 가진 격자점이 슬라이스를 옮겨 가며 옆으로 얼마나 미끄러지는가.
$\gamma_{ij}$ — 슬라이스에 유도된 3-계량. $K_{ij} = -\tfrac12 \mathcal L_n \gamma_{ij}$ 가 외재곡률, 즉 슬라이스가 4차원 안에서 어떻게 휘어 박혀 있는가.

$\alpha,\beta^i$ 는 물리가 아니라 좌표 선택이다. 4개의 게이지 자유도가 여기 전부 들어 있고, 방정식이 이들의 시간 발전을 정해 주지 않으므로 사용자가 직접 처방을 써넣어야 한다. 이것이 다른 어떤 편미분방정식 분야에도 없는 부담이다.

10개의 아인슈타인 방정식은 두 갈래로 갈린다.

구속조건(4개, 시간 미분이 없다):

R + K^2 - K_{ij}K^{ij} = 16\pi\rho, \qquad D_j\!\left(K^{ij}-\gamma^{ij}K\right) = 8\pi S^i

발전방정식(12개): $\partial_t \gamma_{ij}$ 와 $\partial_t K_{ij}$ 를 준다.

구속조건은 초기 슬라이스에서 반드시 성립해야 하고, 성립하면 발전에 의해 해석적으로는 계속 성립한다(비앙키 항등식). 문제는 수치적으로는 그렇지 않다는 것. 반올림 오차가 구속조건을 벗어난 성분을 만들고, 표준 ADM 형식에서 그 성분은 지수적으로 자란다. 1990년대 내내 이진 블랙홀 코드가 궤도 한 바퀴를 못 채우고 죽은 원인이다.

3. 쌍곡성 — 왜 ADM 은 안 되는가[편집]

핵심 진단은 1990년대 후반에 나왔다. 표준 ADM 형식은 약쌍곡형(weakly hyperbolic) 이고, 약쌍곡형 시스템은 잘 정립된 초기값 문제가 아니다. 주기호(principal symbol)가 대각화되지 않아 조던 블록이 남고, 그 결과 해가 $\|u(t)\|\lesssim e^{Kt}\|u(0)\|$ 형태가 아니라 시간의 거듭제곱 인자를 달고 자란다. 이산화 방식을 아무리 바꿔도 이 성질은 고쳐지지 않는다 — 수치 기법의 문제가 아니라 연속 방정식의 문제이기 때문이다.

해결책은 전부 “같은 물리, 다른 대수”였다.

BSSN(나카무라·오오하라·코지마 1987 → 시바타–나카무라 1995 → 바움가르테–샤피로 1999). 계량을 등각 분해해 $\gamma_{ij}=e^{4\phi}\tilde\gamma_{ij}$ ( $\det\tilde\gamma=1$ ), 외재곡률을 자취 $K$ 와 무자취 $\tilde A_{ij}$ 로 쪼갠 뒤, 결정적으로 등각 접속함수 $\tilde\Gamma^i \equiv -\partial_j\tilde\gamma^{ij}$ 를 독립 변수로 승격한다. 이 승격 덕에 리치 텐서에 있던 계량의 혼합 2계 미분이 사라지고 라플라시안 꼴만 남아 시스템이 강쌍곡형이 된다. 승격된 변수가 원래 정의와 어긋나지 않도록 하는 것 자체가 새로운 구속조건이 되는데, 실무적으로는 이걸 주기적으로 강제하는 편이 안정하다.
일반화 조화 형식(GHG). 좌표에 직접 파동방정식 $\Box x^\mu = H^\mu$ 를 부과한다. 그러면 장방정식이 각 계량 성분에 대한 준선형 파동방정식 꼴이 되어 쌍곡성이 자명해진다. 대가는 게이지를 자유롭게 고를 수 없다는 것 — 소스항 $H^\mu$ 를 어떻게 굴릴지가 새로운 설계 문제가 된다.
Z4 계열(보나 외 2003 → Z4c, CCZ4). 구속조건 위반을 숨기지 않고 아예 전파하는 장 $Z_\mu$ 로 승격한 뒤 감쇠항을 붙여 지수적으로 0으로 끌어내린다. “위반은 어차피 생기니 사라지게 만들자”는 발상이고, 요즘 중성자별 코드에서 널리 쓰인다.

구속조건 감쇠는 GHG 쪽에서도 결정적이었다. 위반 모드에 비례하는 항을 인위적으로 추가해 위반의 부분공간에서만 지수 감쇠를 걸면, 물리적 해는 건드리지 않으면서 표류만 억제할 수 있다. 2005년 돌파의 절반은 이 아이디어였다.

4. 초기 데이터 — 타원형 문제[편집]

시간 적분을 시작하려면 구속조건을 만족하는 $(\gamma_{ij},K_{ij})$ 가 필요하다. 이건 발전과 완전히 성격이 다른, 타원형 경계값 문제다.

표준 도구는 요크–리흐네로비치의 등각 횡무자취(conformal TT) 분해다. 계량을 $\gamma_{ij}=\psi^4\bar\gamma_{ij}$ 로, 외재곡률의 무자취 부분을 종방향(벡터 퍼텐셜의 대칭 미분)과 횡무자취 부분으로 쪼개면, 구속조건이 등각인자 $\psi$ 와 벡터 $W^i$ 에 대한 결합 타원 방정식계로 정리된다. 자유롭게 고르는 것(배경 계량, 횡무자취 성분, $K$ )과 방정식이 정해 주는 것이 깔끔하게 갈린다는 점이 이 분해의 미덕이다.

$\bar\gamma_{ij}=\delta_{ij}$ (등각평탄) + $K=0$ (최대 슬라이싱)을 택하면 운동량 구속이 선형이 되고, 놀랍게도 해석해가 존재한다 — 보웬–요크 해(1980). 임의의 운동량 $P^i$ 와 스핀 $S^i$ 를 가진 블랙홀에 대해

\bar A^{ij} = \frac{3}{2r^2}\Big(P^i n^j + P^j n^i - (\delta^{ij}-n^i n^j)P^k n_k\Big) + \frac{3}{r^3}\Big(\epsilon^{ikl}S_k n_l n^j + \epsilon^{jkl}S_k n_l n^i\Big)

같은 닫힌 형태가 나오고, 중첩도 된다. 남는 것은 해밀토니안 구속뿐인데 이건 $\psi$ 에 대한 비선형 포아송 방정식 꼴이라 수치로 푼다. 펑처 데이터(브란트–브뤼그만 1997)는 $\psi = 1 + \sum_a m_a/(2|\vec r - \vec r_a|) + u$ 로 특이 부분을 해석적으로 떼어내고 매끄러운 $u$ 만 격자에서 푸는 트릭이다. 요즘은 이 타원 문제를 스펙트럼 방법으로 기계 정밀도까지 푸는 코드(TwoPunctures)가 표준이다.

등각평탄 가정에는 대가가 있다. 커 시공간은 등각평탄이 아니다. 그래서 보웬–요크 데이터는 애초에 원하는 물리 상태가 아니고, 적분을 시작하면 그 차이가 정크 복사(junk radiation) — 초기 수십 $M$ 동안 격자를 훑고 지나가는 비물리적 파동 — 로 빠져나간다. 실무에서는 앞부분을 그냥 버린다. 스핀이 큰 경우( $a/M\gtrsim0.93$ ) 보웬–요크로는 도달할 수 없는 영역이 생기는 것도 알려진 한계라, 고스핀 파형에는 커 배경을 쓰는 초기 데이터가 필요하다.

이진계 초기 데이터의 다른 갈래는 확장 등각 얇은 샌드위치(XCTS) 다. 두 개의 이웃한 슬라이스를 함께 지정하는 형식이라 “준평형 궤도”를 요구하기 쉽고, 헬리컬 킬링 벡터 근사와 궁합이 좋다. 중성자별 이진계 초기 데이터는 사실상 전부 이쪽이다.

마지막으로 이심률 제거. 준원 궤도를 원하지만 초기 데이터의 접선속도를 잘못 주면 파형에 이심률 진동이 실린다. 표준 절차는 몇 궤도를 짧게 돌려 궤도반경 진동을 보고 초기 $(\dot r, \Omega)$ 를 뉴턴법 비슷한 반복으로 보정하는 것 — 즉 비싼 함수 평가를 쓰는 저차원 근 찾기 문제다. 보통 서너 번 반복으로 $e\lesssim10^{-4}$ 까지 내린다.

5. 게이지와 특이점 — 움직이는 펑처[편집]

발전 중 좌표를 어떻게 굴릴지가 남았다. 사실상 국룰이 된 조합이 움직이는 펑처 게이지다.

\partial_t \alpha = \beta^i\partial_i\alpha - 2\alpha K \qquad\text{(1+log 슬라이싱)}

\partial_t \beta^i = \tfrac34 B^i, \qquad \partial_t B^i = \partial_t\tilde\Gamma^i - \eta B^i \qquad\text{(감마-드라이버)}

1+log 슬라이싱은 특이점 회피 성질을 갖는다. 특이점 근처에서 $K$ 가 커지면 랩스가 0으로 붕괴해 슬라이스의 전진이 그 지역에서 멈춘다. 최대 슬라이싱( $\nabla^2\alpha = \alpha K_{ij}K^{ij}$ , 타원)도 같은 성질을 주지만 매 스텝 타원 문제를 풀어야 해서 비싸다. 1+log 는 그걸 쌍곡형 처방으로 바꿔 놓은 값싼 근사인 셈. 감마-드라이버의 $\eta$ 는 시프트의 감쇠율이고, 질량에 따라 조정해야 하는 사실상의 튜닝 노브다(질량비가 큰 계에서 이 값을 잘못 잡으면 코드가 죽는다).

특이점 자체는 두 방식으로 처리한다.

엑시전(excision) — 겉보기 지평선 안쪽 영역을 격자에서 통째로 도려낸다. 안쪽은 인과적으로 바깥에 영향을 주지 못하므로 경계조건을 아예 걸지 않아도 되지만, 도려낸 경계가 블랙홀을 따라 움직이므로 격자 관리가 지저분하다. SpEC 같은 유사스펙트럼 코드는 영역 분해와 궁합이 좋아 지금도 이 방식이다.
펑처(puncture) — 특이 부분을 해석적으로 분리하고 나머지만 발전시킨다. 2005~06년의 발견은 움직이는 펑처 게이지를 쓰면 펑처를 격자에 고정할 필요가 없다는 것이었다. 게이지가 슬라이스를 특이점에서 알아서 밀어내므로 펑처는 격자 위를 그냥 흘러 다니고, 코드는 그 위를 아무렇지 않게 유한차분한다.¹ 유한차분 + AMR 코드는 대부분 이쪽이다.

6. 이산화 — 유한차분/AMR 대 유사스펙트럼[편집]

수치 상대론에는 두 개의 학파가 있고, 둘 다 같은 파형을 낸다.

고차 유한차분 + Berger–Oliger AMR. 공간은 6차 또는 8차 중심차분, 시간은 4차 룽게–쿠타(룽게-쿠타법), 여기에 크라이스–올리거 소산을 얹어 격자 진동수 노이즈를 죽인다. 소산 차수는 차분 차수보다 높게 잡아 수렴차수를 깎지 않도록 한다. 블랙홀 주변만 촘촘히 봐야 하므로 적응 격자 세분화가 필수이고, 실무에서는 보통 “고정 격자 세분화” — 블랙홀을 따라다니는 중첩된 정육면체 상자들 — 로 충분하다. 파동을 손실 없이 멀리 보내야 하므로 파장당 격자점 수가 정확도를 지배하며, 이게 계산량을 잡아먹는 1순위다.

다영역 유사스펙트럼. 도메인을 구각·원통·정육면체 등으로 쪼개고 각 영역에서 체비쇼프 다항식이나 구면조화 기저로 전개한다. 매끄러운 해에 대해 지수 수렴하므로 같은 정확도를 훨씬 적은 자유도로 얻는다. 대신 격자 관리가 복잡하고 충격파가 있으면 무력하다. SXS 협력단의 SpEC 이 대표 주자이며, 긴 파형(수십~수백 궤도)을 낮은 비용으로 뽑는 데 강점이 있다. 자세한 이론은 스펙트럴 방법·의사스펙트럼법 참고.

바깥 경계도 골칫거리다. 이상적으로는 구속조건을 보존하면서 들어오는 중력파를 막는 경계조건이 필요한데, 순진한 유출 조건은 구속조건 위반을 안쪽으로 반사한다. 대안은 세 가지 — 경계를 아주 멀리 두고 파형 추출 시각 전에 반사가 도달하지 못하게 하기(무식하지만 널리 쓰임), 구속조건 보존 경계조건, 그리고 쌍곡형 슬라이싱/특성곡면 발전으로 널 무한대를 격자에 포함시키기.

7. 파형 추출과 카탈로그[편집]

시뮬레이션의 최종 산출물은 중력파형이다. 표준 경로는 뉴먼–펜로즈 스칼라 $\Psi_4$ 를 추출 구면에서 계산하고 스핀가중 $-2$ 구면조화 함수로 다극 전개하는 것:

\Psi_4(t,\theta,\varphi) = \sum_{\ell m} \psi_{\ell m}(t)\; {}_{-2}Y_{\ell m}(\theta,\varphi), \qquad \Psi_4 = \ddot h_+ - i\,\ddot h_\times

여기서 두 번의 시간적분이 필요한데, 이게 순진하게 하면 저주파 표류로 파형이 위아래로 휜다. 표준 해법이 고정주파수 적분 — 주파수 영역에서 $1/(2\pi i f)^2$ 를 곱되 궤도 최저 주파수 아래는 잘라내는 것 — 이다. 신호처리 쪽에서 보면 그냥 정규화된 이중적분기이고, 실제로 푸리에 변환 기반 필터링 문제와 다르지 않다.

두 번째 문제는 유한 반경이다. 계량 섭동으로서의 파형은 널 무한대에서만 게이지 불변인데, 격자는 유한하다. 처리는 여러 반경에서 뽑아 $1/r$ 로 외삽하거나, 코시–특성곡면 추출(CCE) 로 코시 발전의 결과를 특성곡면 코드에 넘겨 널 무한대까지 실제로 전파시킨다. 후자가 원리적으로 옳고, 요즘 카탈로그는 CCE 파형을 함께 제공한다. 메모리 효과처럼 진짜로 무한대에서만 정의되는 양은 이 경로가 아니면 얻을 수 없다.

이렇게 만든 파형은 SXS, RIT, 조지아텍 등의 공개 카탈로그에 수천 편 규모로 쌓여 있다. 이것이 검출기 데이터 분석으로 흘러 들어가는 방식은 블랙홀 문서에 정리돼 있는데, 요점만 다시 적으면 이렇다 — 파형 하나에 수만~수백만 코어시간이 들기 때문에 탐색용 템플릿을 매번 시뮬레이션할 수 없다. 그래서 수치 파형은 EOB·현상론 모형·축소기저 대리모형을 교정하는 데 쓰이고, 실제 탐색은 그 값싼 모형이 만든 템플릿 뱅크로 돌린다. 수치 상대론은 데이터 분석의 실행 엔진이 아니라 눈금자다.

8. 중성자별 — 물질이 붙으면[편집]

블랙홀 진공 문제에 물질이 들어오면 GRMHD 가 된다. 시공간 발전(보통 BSSN 또는 Z4c)과 상대론적 자기유체역학을 매 스텝 결합해 푸는데, 유체 쪽은 유한체적법 + 고분해능 충격파 포착(리만 문제 근사 해법)을 쓴다. 즉 한 코드 안에 매끄러운 장을 위한 고차 중심차분과 불연속을 위한 풍상 리만 솔버가 공존한다.

고통 지점은 유체 쪽에 몰려 있다.

보존량 → 원시변수 역변환이 닫힌 형태가 아니다. 매 격자점마다 비선형 방정식을 푸는데, 강자기장·저밀도에서 이 반복이 실패하는 것이 코드 사망 1순위.
상태방정식이 물리 입력이자 최대 불확실성이다. 유한 온도·조성 의존 표를 넣어야 하고, 표 보간의 열역학적 일관성이 깨지면 조용히 틀린 답이 나온다.
중성미자 수송이 킬로노바 광도와 원소 합성을 좌우한다. 완전한 볼츠만 수송은 감당이 안 되므로 누출(leakage) 근사나 M1 모멘트 닫힘을 쓴다 — 닫힘 관계가 곧 체계 오차다.

GW170817 이후 이 계열의 수요가 폭발했다. 파형에서 조석 변형률 $\Lambda$ 를 읽어 상태방정식을 제약하고, 병합 잔해에서 나오는 방출물의 조성과 광도곡선을 킬로노바 관측과 맞춘다. 중력파와 전자기 관측이 같은 시뮬레이션 출력을 공유하는 유일한 분야라, 다중신호 천문학의 계산 허브 역할을 한다.

9. 여담[편집]

이 분야의 표준 검증법은 수렴 테스트다. 해상도를 1.5배씩 올리며 구속조건 위반의 노름이 이론 차수대로 줄어드는지 보고, 줄어들지 않으면 어딘가 차수가 깨진 것이다. 해석해가 없는 문제에서 검증 및 확인을 하는 교과서적 사례라, 다른 분야에서도 그대로 베껴 쓸 만하다.
아인슈타인 툴킷, GRChombo, SpECTRE, AthenaK 등이 오픈소스로 공개돼 있다. 다만 “빌드가 된다”와 “수렴한 파형이 나온다”는 아주 다른 사건이다.
이진 블랙홀 병합의 최종 스핀·최종 질량은 수천 편의 시뮬레이션에 대한 피팅 공식으로 정리돼 있다. 90년을 갈아 넣은 결과물의 상당 부분이 결국 몇 줄짜리 대수식으로 요약되어 유통되는 셈인데, 이건 낭비가 아니라 그 공식이 어디까지 믿을 만한지를 아는 것이 시뮬레이션의 진짜 산출물이기 때문이다.²
극단 질량비( $q \gtrsim 10^3$ , LISA 표적)는 아직도 정면 돌파가 안 된다. 작은 천체를 격자에서 분해하려면 해상도 요구가 폭증하고, 궤도 수는 반대로 수만 바퀴가 된다. 여기는 자기력(self-force) 섭동론의 영역이고, 두 방법이 겹치는 중간 질량비 구간에서 서로를 검증하는 중이다.³

10. 관련 문서[편집]

“특이점 위를 유한차분한다”는 말이 처음 들으면 미친 소리 같은데, 실제로 펑처 격자점에는 값을 아예 두지 않거나(격자를 반칸 어긋나게 배치) 게이지가 만든 이른바 “트럼펫” 기하 덕에 그 지점의 장이 유한해진다. 2005년 이전 10년간 모두가 펑처를 얌전히 고정해 두려고 애썼는데, 정답은 놔주는 것이었다. 분야 전체가 뒤통수를 맞은 사건. ↩
이 피팅 공식들이 실제로 일하는 곳은 데이터 분석이다. 사후분포에서 최종 질량·스핀을 계산할 때마다 시뮬레이션을 돌릴 수는 없으니까. 대신 공식의 유효 범위 밖(초고스핀, 극단 질량비)에서 태연히 외삽하는 논문이 주기적으로 나오고, 그때마다 리뷰어가 잡는다. ↩
자기력 계산은 발산하는 자기장을 정규화해서 유한한 힘을 뽑아내는 작업이라 이름부터 불온한데, 정규화 처방(모드합, 유효장)이 여럿 있고 서로 일치한다는 것이 확인되기까지도 한참 걸렸다. 큰 그림에서는 후-뉴턴·자기력·수치 상대론이 각자의 정의역에서 계산한 값이 겹치는 구간에서 만나야 한다는 요구가 이 분야 전체의 상호 검증 구조다. ↩