태양 질량의 1.4~2.1배에 달하는 거대한 물질이 지름 약 20km의 초고밀도 구체에 응축된 중성자별(Neutron Star)은 원자핵 포화 밀도($\sim 2.8 \times 10^{14}\,\text{g/cm}^3$)를 초과하는 우주 유일의 극한 물질 실험실입니다. 일반상대성이론의 아인슈타인 중력장 방정식을 구대칭 정역학 평형 상태로 유도한 톨만-오펜하이머-볼코프(Tolman-Oppenheimer-Volkoff, TOV) 연립 비선형 미분방정식을 핵물질 상태방정식(EoS)과 결합하여 중심 밀도별 질량-반경(M-R) 궤적을 1,000Hz 주기로 0-Alloc C# 적분하는 상대론적 시뮬레이터를 공개합니다.

1. 2022 개정 교육과정 교과 연계 분석

  • 물리학 (일반상대성이론 중력장 방정식, 시공간 곡률 및 상대론적 유체역학): 슈바르츠실트 내부 시공간 계량 텐서 $ds^2 = -e^{2\Phi} c^2 dt^2 + e^{2\Lambda} dr^2 + r^2 d\Omega^2$, 압력 $P$에 의한 시공간 곡률 기여 및 상대론적 중력 강화 계수.
  • 핵물리학/천문학 (페르미 축퇴압, 파울리 배타 원리 및 핵자 상호작용): 비상대론적/상대론적 중성자 페르미 에너지 축퇴압 $P \propto \rho^{5/3}$ 및 $P \propto \rho^{4/3}$, 중간자 교환 상호작용(RMF 모델)에 의한 최대 허용 질량(오펜하이머-볼코프 한계 $\sim 2.2 M_\odot$).
  • 수학 (비선형 연립 상미분방정식 경계치 문제 및 4차 룬게-쿠타 수치 적분): 중심점($r=0, P=P_c, M=0$)에서 별 표면($P(R)=0$)까지의 TOV 2원 1계 비선형 ODE 4차 룬게-쿠타(RK4) $O(N)$ 적분.
  • 정보 (초저지연 1,000Hz 0-Alloc 천체물리 수치적분 코어): 가비지 컬렉터 부하 없이 스택 메모리(`ref struct`, `Span`)를 활용하여 NASA NICER X선 관측 데이터와 중력파 조석 변형률($\Lambda$)을 1ms 내에 실시간 상호 검증하는 아키텍처.

2. 톨만-오펜하이머-볼코프(TOV) 상대론적 정역학 수식 모델

일반상대성이론에서 정역학적 평형을 이루는 완전유체 구체의 압력 구배 $\frac{dP}{dr}$와 봉입 질량 $M(r)$의 시간 비의존 TOV 연립 미분방정식은 다음과 같습니다:

$$\frac{dP(r)}{dr} = -\frac{G M(r) \rho(r)}{r^2} \left[ 1 + \frac{P(r)}{\rho(r) c^2} \right] \left[ 1 + \frac{4\pi r^3 P(r)}{M(r) c^2} \right] \left[ 1 - \frac{2GM(r)}{c^2 r} \right]^{-1}$$ $$\frac{dM(r)}{dr} = 4\pi r^2 \rho(r)$$

뉴턴의 고전 중력 정역학식 $\frac{dP}{dr} = -\frac{G M \rho}{r^2}$과 비교할 때 대괄호 안의 세 가지 상대론적 보정 인자는 중심부의 극단적 압력과 에너지 밀도가 추가적인 인력 질량으로 작용함을 명확히 보여줍니다.

다방성(Polytropic) 핵물질 상태방정식 $P = K \rho^\Gamma$ 및 고밀도 쿼크 물질 전이를 반영하여 중심 압력 $P_c$를 변화시키며 적분을 수행하면, 별의 총 질량 $M_{total}$과 반경 $R_{star}$ 간의 비선형 안정성 한계 곡선이 도출됩니다.

3. 1,000Hz 0-Alloc C# TOV 수치적분 엔진

다음은 1ms 제어 주기 내에서 가비지 컬렉터 없이 스택 메모리를 사용하여 TOV 연립 미분방정식을 4차 룬게-쿠타로 적분하고 최대 질량 한계를 평가하는 C# 엔진입니다.

TovRelativisticEngine.cs (Zero-Allocation C# Relativistic Core) 
using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Astrophysics;

public readonly ref struct NeutronStarProfile
{
    public readonly double TotalMassSolarMasses;
    public readonly double RadiusKilometers;
    public readonly double CompactnessRatio;
    public readonly double SurfaceRedshift;
    public readonly bool IsRelativisticallyStable;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public NeutronStarProfile(double massSolar, double radiusKm, double compactness, double redshift, bool stable)
    {
        TotalMassSolarMasses = massSolar;
        RadiusKilometers = radiusKm;
        CompactnessRatio = compactness;
        SurfaceRedshift = redshift;
        IsRelativisticallyStable = stable;
    }
}

public ref struct TovRelativisticSimulator
{
    private const double G = 6.67430e-11;
    private const double C = 299792458.0;
    private const double SolarMassKg = 1.98847e30;
    private const double Dr = 10.0; // 10미터 공간 적분 스텝

    // 폴리트로프 상태방정식 상수 (핵물질 근사 P = K * rho^Gamma)
    private const double Gamma = 2.35;
    private const double K = 1.2e-4;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public NeutronStarProfile IntegrateStar(double centralDensityGcm3)
    {
        double rhoC = centralDensityGcm3 * 1000.0; // kg/m^3
        double pC = K * Math.Pow(rhoC, Gamma);

        double currentR = 10.0; // r=0 특이점 회피
        double currentM = (4.0 / 3.0) * Math.PI * currentR * currentR * currentR * rhoC;
        double currentP = pC;

        // 표면(P <= 0)까지 최대 3,000 스텝 적분 (최대 30km)
        for (int step = 0; step < 3000; step++)
        {
            if (currentP <= 1.0e-5) break;

            double rho = Math.Pow(currentP / K, 1.0 / Gamma);

            // TOV 미분계수 산출
            double cSq = C * C;
            double factor1 = 1.0 + currentP / (rho * cSq);
            double factor2 = 1.0 + (4.0 * Math.PI * currentR * currentR * currentR * currentP) / (currentM * cSq);
            double metricTerm = 1.0 - (2.0 * G * currentM) / (cSq * currentR);

            if (metricTerm <= 0.05) break; // 사건의 지평선 형성 (블랙홀 붕괴 한계)

            double factor3 = 1.0 / metricTerm;
            double dPdr = -((G * currentM * rho) / (currentR * currentR)) * factor1 * factor2 * factor3;
            double dMdr = 4.0 * Math.PI * currentR * currentR * rho;

            currentP += dPdr * Dr;
            currentM += dMdr * Dr;
            currentR += Dr;
        }

        double totalMassSolar = currentM / SolarMassKg;
        double radiusKm = currentR / 1000.0;
        double compactness = (G * currentM) / (currentR * C * C);
        double surfaceRedshift = 1.0 / Math.Sqrt(Math.Max(0.01, 1.0 - 2.0 * compactness)) - 1.0;

        // 오펜하이머-볼코프 안정성 한계 판정 (질량 < 2.3 M_sun)
        bool stable = totalMassSolar > 1.0 && totalMassSolar < 2.35 && compactness < 0.35;

        return new NeutronStarProfile(totalMassSolar, radiusKm, compactness, surfaceRedshift, stable);
    }
}

4. 우주천문학 및 첨단 물리 교육 응용 가이드

  • 한국천문연구원(KASI) 및 고에너지 천체물리 연구진: NICER X선 망원경 펄사 프로파일 모델링 및 LVK 중력파 관측소 킬로노바 잔해물 상태방정식 제약 검증.
  • 방통대 우호성 교수 랩실: 2022 개정 물리학Ⅱ 일반상대성이론 및 만유인력 단원을 초고밀도 중성자별 시뮬레이션과 결합한 현대 천체물리 가상 실습 교구.
  • 국방 위성항법 및 양자중력 계측 분야: 극한 중력장 시간 지연(중력 적색편이) 효과 보정 알고리즘 및 펄사 기반 심우주 독자 항법 시스템 기초 설계.