[ABP: 중력파 관측을 위한 정밀도는 얼마인가]
초정밀의 세계
학부 때 자아가 하늘을 찌르던 어떤 시절, 원주율을 소수점 아래 100자리까지 외우고 다니곤 했습니다. 아마 나보다 더 많은 자리수까지 원주율을 외우고 다닌 또래는 없었을 것이라 자신하면서 언젠가 꼭 뽐낼 자리가 있으려니 생각하기도 했죠. 하지만 정작 한국에서는 쓸 일이 없다가 미국에서 박사 과정 시작할 때 기회가 왔습니다. 100자리 까지는 다 못 외우고 있었지만 대략 30자리까지는 외웠던 것으로 기억하는데, 같이 수업을 듣던 물리학과 친구가 그냥 아무 표정 없이 150자리까지인가 제 앞에서 배틀하듯 외우는 모습을 보고 무릎을 꿇었습니다. 그렇습니다. 세상에는 고수가 늘 있기 마련입니다. 오늘 3월 14일은 파이 (pi) 데이입니다 (화이트 데이가 아닙니다.). 예전에 미국에 있을 때 수업듣던 수업의 물리학과 선생님이 해 주셨던 말씀이 생각납니다.
"수치를 계산할 때 학생들은 그냥 수치해석 엔진을 아무 생각없이 믿는 경향이 있다. 편하니까 그런거 이해하는데, 정밀한 실험물리학에서는 유효숫자 자리수가 정말 중요해지는 때가 있다. 원주율만 해도 그렇다. 너희 수치해석 엔진에서 원주율을 어떻게 처리한다고 생각하냐? 대충 3.141592... 뭐 이런 숫자를 쓴다고 생각하냐? 원주율의 유효숫자가 몇 자리수만 앞으로 땡겨오면 인류가 지금 쏘아올리는 우주선이 태양계를 벗어날 때쯤에는 오차가 어마어마해진다. 예를 들어 가장 가까운 항성계인 켄타우르스 알파별로 항로를 정했다 치자. 그런데 3.141592와 3.14159265358979정도의 격차는 그 별에 도달할 수 있는지 없는지의 차이를 만들어낼 수도 있다."
저는 이것이 도대체 무슨 이야기인지 잘 몰랐는데, 나중에 고급 수치해석 강의를 듣다가 행렬의 컨디션수 (condition number) 개념을 다시 배우며 비로소 조금씩 이해하기 시작했습니다. 예를 들어 어떤 변환 (projection or transformation)을 하는 행렬 A를 생각해 봅시다. 이 행렬의 역행렬을 수치해석방법으로 구하는 과정에서 실제 역행렬과 행렬 E만큼의 에러가 생긴다고 생각해 봅시다. 이를 P = inv(A), Q = inv(A+E) 라고 표현할 수 있을 것입니다. 그러면 행렬의 노름 (norm)을 이용하여 norm(Q-P)/norm(P) 값을 계산할 수 있을 것입니다. 이 값의 상한선은 행렬 A의 condition number kappa(A)와 norm(E)/norm(A)의 곱으로 정해집니다. 흥미롭게도 수치해석 엔진 처리 과정에서 발생하는 에러 행렬의 norm (즉, norm(E))은 미리 정해진 계산 오차 정밀도 (precision)에 따르게 되는데, 예를 들어 부동소숫점 (floating point) 정확도를 2^-52 정도로 설정한 상황이라면 norm(E)는 행렬 크기*floating point precision*norm(A)로 주어집니다. 예를 들어 2^-52의 floating point precision은 십진법에서는 2.2*10^-16 정도의 계산 오차 정밀도로 환원됩니다. 즉, 수치해석 엔진에서 얻어지는 수치들의 정밀도는 최대 소숫점 16자리까지 얻어질 수 있는 셈이 됩니다.
이해를 돕기 위해 A행렬을 2*2 행렬인 [4.1 2.8; 9.7 6.6]으로 생각해 봅시다. 이 경우 cond(A) = 2.2*10^3 이므로 norm(Q-P)/norm(P) <=2*2.2*10^3*2.22*10^-16 = 9.98*10^-13 으로 얻어집니다. 즉, 이 수치해석 엔진에서 얻어지는 행렬 A의 역행렬을 계산할 때, 그 허용 오차 정밀도 (계산 정확도)는 대략 적어도 소숫점 아래 13자리까지는 믿을 수 있는 셈이라는 뜻입니다. 실제로 norm(Q-P)/norm(P) = 8.7*10^-16 으로 나오긴 하는데, 어쨌든 안전한 한계를 알아낼 수 있다는 것은 수치해석 계산에서는 큰 이득입니다. 왜냐하면 어디까지 계산된 수치의 정밀도를 '일단' 믿고 진행할 것인지를 결정할 수 있기 때문입니다.
간혹 이러한 instrumental or algorithmic precision의 한계를 고려하지 않고 계산된 수치를 그냥 쓰는 경우가 생기는데, 대부분의 경우는 별탈이 나지 않지만, 극도의 정밀도를 요구하는 경우에는 이것이 문제가 될 수 있습니다. 예를 들어 최첨단 반도체 제조 공정에서 핵심 공정인 초미세 패터닝 공정의 총아 EUV lithography에서는 EUV 전자기파를 한 군데로 모으기 위해 전용 반사경이 필요한데, 적절한 광경로를 형성하기 위해서는 이들의 곡률 profile은 물론, bragg mirror로 쓰기 위한 소재 geometry가 극도의 정밀도로 가공되어야 합니다. 옹스트롬 (0.1 nm) 이하는 물론, 믿을 수 없을 정도의 정밀도인 10^-3 옹스트롬 (0.0001 nm)까지도 내려가야 하는 수준입니다. 보통 원자 한 개의 지름이 수 옹스트롬이니까 원자를 적어도 만 번 쪼개는 수준의 정밀도가 필요한 셈입니다. 이 정밀도의 제어는 당연히(?) 수치해석적으로 상황에 맞게 조정해야 (adaptive 하게) 해야 하는데, 그렇다면 그 제어 프로그램의 정밀도는 이보다 훨씬 높아야 할 것입니다. 그런데 상용 adaptive optics, frontline optics 광경로 제어 알고리즘은 이 정밀도를 상정하지 않기 때문에 이를 다시 처음부터 floating point 고려해가면서 짜야 합니다. 짜이스는 이 때문에 fine optics 제어 소프트웨어를 처음부터 다시 만들었습니다.
중력 간섭계 LIGO도 마찬가지 케이스입니다. 기본적으로 과거 빛의 속도를 정확하게 측정하기 위해 고안된 마이컬슨 간섭계 (Michelson interferometer)를 거대하게 확장한 구조를 갖는 LIGO는 블랙홀-블랙홀, 혹은 블랙홀-중성자별, 혹은 중성자별-중성자별 같은 거대한 천체들의 충돌과 합쳐짐에 따른 중력의 변화, 그리고 그에 따른 에너지 소멸과 시공간의 왜곡에서 필히 동반되는 중력파 (gravitational wave)를 아주 먼거리에 있는 검지기, 즉, 지구 상의 검지 장치에서 측정하기 위해 고안된 극도로 정밀한 실험 장치입니다. LIGO가 측정할 수 있는 공간 팽창 혹은 수축에 따른 변형률 (strain)의 측정 정밀도는 10^-21이나 됩니다 (더 정확히는 10^-21 1/sqrt(Hz)). 변형률이 10^-21이라는 뜻은 길이로 환산하면 얼마 정도 될까요? LIGO에 필요한 레이저 간섭계 한 축의 길이는 대략 4 km 정도 됩니다. 이 길이의 10^-21 변형률은 대략 원자 크기의 1/1000 정도 의 변형을 의미합니다. 신호 증폭을 위해 두 간섭계 사이에 레이저가 수천 번 왕복하는 상황을 고려한다고 해도 10^-21의 변형률은 원자 한 개 정도의 크기까지 정밀해져야 얻을 수 있는 극한 수준의 정밀도입니다. 이 정도 precision이라면 그에 걸맞는 데이터 처리 알고리즘도 개선되어야 할 것입니다. 기존의 수치해석 알고리즘에서 차용하던 디지털 신호 처리의 2^-52 수준의 precision은 어림도 없고, 적어도 2^-67 수준까지는 올라가야 할 것입니다. 노이즈가 잔뜩 끼어 있는 신호를 처리하기 위해서는 다양한 필터 (filter)를 활용해야 하는데, 이 필터들의 모든 수치해석 코드들도 다시 손봐야 함을 의미하기도 합니다.
LIGO에서 관측한 실험 측정값이 실제로 공간의 수축이라는 것을 확인하기 위해, 아인슈타인의 일반상대성 이론에 따라 모델링한 예측치도 필요합니다. 물리학자 Pretorius는 근사적 방법의 도움 없이 중력파의 특성을 수학적으로 계산할 수 있는 방법론을 개발하여 지난 2005년에 Phys. Rev. Lett.에 보고하기도 했습니다. (https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.95.121101?fbclid=IwAR1namzuxMRCOZl8Y5DhuGVHV2YyzkLBloBDkKjcI2zktS9rcZv0WRqFzeI). 이 논문은 중력파의 대략적인 특성을 알려주는데는 큰 도움이 되었지만 정확히 강력한 중력을 가진 천체의 충돌 과정과 에너지 변환, dynamics 등에 대한 계산은 보다 정밀한 수치해석적 방법을 필요로 했습니다. 특히 중력이 매우 강한 상황에서는 시공간이 휘게 되므로 예측 모델링 역시 이전의 수치해석 알고리즘 precision으로는 원하는 수준의 정밀도를 얻어낼 수 없습니다.
그러면 어떻게 해야 할까요? 일반상대성 이론에 따른 블랙홀 충돌-병합 과정의 동역학은 미국 NASA 고다드 우주센터 (NASA Goddard Space Flight Center)에서 근무하던 최대일 박사 (Dr. Dale Choi)라는 분이 개발한 수치해석 알고리즘 (Hahndol code)으로 계산될 수 있습니다. 시공간이 일정하다고 가정하고 계산하는 유한차분법과는 달리, 이 현상은 시공간이 일정하지 않고 거대한 중력을 가진 천체에 의해 '휘는' 것도 고려해야 합니다. 이를 위해 Mesh-adapted difference 방법이 필요한데, 이는 기존의 크랭크-니콜슨 (Crank-Nicholson) 시간 적분기와 유사하지만 비선형 고차 미분항 적분에 특화된 방식을 필요로 합니다. 특히 천체 주변의 공간을 나누는 메쉬는 대략 1억 2천만 개의 finite element로 나뉘어 추적되었는데, 이 과정에서 특히 천체 사이의 좁은 영역의 경계를 어떻게 처리할 것인지가 수치해석 적분의 관건이 됩니다. 한돌 코드에서는 블랙홀 한 쌍으로 인한 공간 왜곡 과정에서의 경계 문제를 풀기 위해 Bowen-York 알고리즘을 활용했고, 중력파를 waveform 여러 개로 쪼개어 이들의 수렴성 (convergence)를 테스트했습니다. 수치해석 결과의 정확도를 위해 그 수렴판단 기준 (convergence threshold)은 컴퓨터 자원이 허용하는 한도 내에서 계속 증가했고, 자체적으로 수렴성이 확인되는 수치에서 시뮬레이션 예측이 이뤄졌습니다.
이렇게 예측한 두 블랙홀의 충돌 시뮬레이션 결과는 LIGO 관측 결과와 각각의 precision 안에서 정확히 일치했고, 이는 LIGO 관측이 2년 후 노벨 물리학상을 수상하는 결정적인 근거를 제공했습니다 (2017년 노벨 물리학상 수상자는 Thorne, Weiss, Barish).
LIGO 외에도 작년부터 어마어마한 데이터를 쏟아내고 있는 제임스웹우주망원경 (JWST)이나 유럽의 거대강입자가속기 (LHC) 같은 현대 물리학 연구의 최전선에는 결국 실험적으로 측정한 데이터의 정밀도 한계와의 전쟁이 벌어지고 있다고 해도 과언이 아닙니다. 정밀도가 한 자릿수 올라갈 때마다 인류가 자연을 이해할 수 있는 범위는 넓어지면서 불확설성의 범위는 좁아진다고 해도 틀린 말은 아닐 것입니다. 진리로 가는 길이 보이기 시작하며, 곁가지 길로 빠질 위험도 줄어들기 때문입니다.
오늘 pi 데이를 맞아 매년 원주율 계산 결과를 끝없이 경신하는 수 많은 수학자들, 아마추어들, 그리고 precision과 매일 싸움을 이어가는 수 많은 과학자들과 연구자들, 엔지니어와 기술자들에게 응원과 감사의 말씀을 드립니다. 문자 그대로 precision의 문제는 인류 문명의 최전선의 문제나 마찬가지이고 이들은 그 최전선에서 인류를 대신하여 인류의 영역을 넓혀가고 있다고 해도 과언이 아닐 것입니다.
