复合多方模型下致密星的随机分析

《Scientific Reports》：Stochastic analysis of compact stars under composite polytropes

　　

致密物质研究因理论高能模拟与现代观测天文学的进步而取得显著进展。中子星（Neutron Star, NS）作为探测极端密度与强引力场下物质行为的天然实验室，其结构建模依赖于物态方程（Equation of State, EoS）的合理选取。多方关系（Polytropic Relation）P=Kρ^1+1/n作为恒星天体物理学中的基本参数化EoS，因其简洁性与物理直观性而被广泛应用，其中n为多方指数，K为依赖于组分与熵的常数。然而，真实致密星体具有显著的内部分层结构，单一多方指数无法充分表征这种层状特征。为此，研究人员构建了复合相对论多方模型，使多方指数随半径r平滑变化，通过双曲正切函数n(r)=η₀?η₁tanh[(r?r_c)/ε]实现核心与包层之间的连续过渡，其中r_c为核心边界，ε为过渡宽度。

研究人员将爱因斯坦场方程（Einstein Field Equations）与广义复合多方EoS耦合，推导得到复合Tolman–Oppenheimer–Volkoff方程组（Composite TOV, CTOV）。该方程组为耦合的非线性常微分方程组，其解析求解通常不可行。传统确定性数值求解器（如Runge-Kutta方法）在面对此类刚性方程、参数敏感性及不确定性量化需求时存在局限。为此，研究人员采用基于蒙特卡洛（Monte Carlo, MC）的数值积分方法，将积分过程视为随机采样过程，通过 hit-or-miss 策略处理方程刚性，同时实现对参数空间的概率化探索与自然的不确定性量化。

该研究的核心目标包括：验证MC方法求解复杂非线性CTOV系统的有效性；量化复合模型参数对致密星质量-半径关系及内部结构的预测影响；提供能够自然输出预测不确定性的框架，以反映模型对输入假设的敏感性。

研究结果表明，相对论参数σ的增大显著降低Emden函数θ与质量函数ν的数值，形成更为致密的恒星构型。更尖锐的核心-包层过渡（ε=0.01）系统性地产生比平滑过渡（ε=0.03）更高的致密性。推导得到的质量-半径关系成功复现了观测到的中子星多样性，既匹配PSR J0030+0451等低质量、大半径系统，也吻合PSR J1614?2230等高质量致密脉冲星。特别地，刚性复合构型（n_c=1, n_e=2, x_c=0.7）的最大质量分析显示，ε=0.01时M_max≈3.47 M_⊙，ε=0.03时M_max≈2.71 M_⊙，对应最小半径R_min≈10.6–13.2 km，与当前观测约束一致。这些发现证实复合多方过程为分层致密星建模及约束高密度物态方程提供了灵活且物理动机充分的框架。

致密天体物理学汇聚了从最强引力吸引至近光速物质喷流等极端物理现象，为天体物理、核物理及引力理论检验提供了天然实验室。中子星、白矮星等致密遗迹的结构建模核心在于物态方程（Equation of State, EoS），其关联恒星内部压强与密度。多方关系（Polytropic Relation）P=Kρ^1+1/n作为基本参数化EoS被广泛应用，其中n为多方指数，K为依赖于组分与熵的常数。然而，真实中子星并非均匀球体，其内部存在受不同微观物理过程调控的分层结构，单一多方指数无法充分表征这一特征。因此，复合或多方模型成为必要，其中多方指数n被允许随径向坐标r变化，以有效拼接不同EoS段表征层状内部。对于此类恒星，相应的CTOV方程组为耦合非线性常微分方程组，其解析求解通常不可行，数值积分至关重要。

传统数值求解器（如Runge-Kutta方法）面临多重挑战：方程在中心与表面区域呈现刚性，解对中心密度ρ_c、内部过渡锐度及相对论参数σ等参数选择高度敏感；且经典确定性方法仅提供单一解轨迹，缺乏对模型预测不确定性的固量化框架。此外，致密星内部可能存在压强各向异性，源于高磁场、相变、固体壳层、超流性或π/K介子凝聚体等效应，导致径向压强P_r与切向压强P_t出现显著差异。尽管本研究聚焦于各向同性情形，所提出框架可扩展至包含此类效应。

为解决CTOV系统求解挑战并稳健探索参数空间、量化不确定性，研究人员转向随机数值方法。MC算法设计用于从复杂高维概率分布中采样，特别适用于参数退化、似然面复杂的问题。在天体物理学中，MC已成为参数估计与模型比较的标准工具。本研究将MC方法直接应用于CTOV方程组的数值求解，而非仅用于方程求解后的后验参数估计。通过将微分方程积分视为随机采样过程，研究人员获得了密度、压强与质量剖面的稳健数值解，并概率化地探索参数空间，系统研究了中心密度ρ_c、核心边界r_c、过渡宽度ε及相对论参数σ等关键物理参数对恒星结构及全局性质（总质量M与半径R）的影响。

研究采用的关键技术方法如下：基于光滑变化多方指数的复合CTOV理论构建；将MC hit-or-miss方法适配于常微分方程初始值问题的数值积分，以步长Δx=0.01离散化无量纲半径x，每步生成N=2×10⁶个均匀随机样本；将物理参数（ρ_c, n_c, n_e, x_c, ε, σ）视为具有物理动机先验分布的随机变量进行采样；执行m=10次独立MC运行生成解集合，通过线性插值至公共径向网格后计算集合均值与标准差；与四阶Runge-Kutta参考解对比验证，收敛测试确认误差按O(1/√N)缩放。

研究结果分为以下部分：

Emden与质量函数分析：牛顿极限（σ=0）下，Emden函数θ从中心θ=1单调递减至表面θ=0。相对论参数σ增大导致θ更快速下降，例如均匀指标n_c=n_e=2在x=2.0处θ从牛顿情形的0.60降至σ=0.5时的0.15，表明更强相对论压缩与更致密构型。集合解显示10次独立MC运行产生的θ(x)与ν(x)曲线在核心区域几乎不可区分，近表面处展开形成不确定性带。模型1（n_c=0.5, n_e=1, x_c=0.2）与模型2（n_c=1, n_e=2, x_c=0.7）的比较揭示：σ从0增至0.5时，模型1的x_s从3.14收缩至1.859（降幅40%），ν(x_s)从2.3734降至0.497（降幅79%）；模型2的x_s从3.074降至1.7657（降幅44%），ν(x_s)从2.558降至0.4567（降幅85%）。过渡宽度ε显著调制响应：ε=0.01的锐利过渡产生比ε=0.03平滑过渡系统性地更高的致密性。模型2的更大核心产生更陡峭的初始质量增长，随后包层中更渐进的增加。

稳定性分析：所有模型满足密度与压强的正性、有限性及单调递减性。绝热指数Γ(x)>4/3确保对径向振荡的稳定性。因果性条件v_sc=1, n_e=2）全参数范围内满足，但在软模型（n_c=0.5, n_e=1）σ>0.3时于内核区域出现超光速，限制其物理允许参数范围为σ≤0.3。Zeldovich条件0<ω_c=P_c/(ρ_cc²)=σ<1在所有模型中成立。零能量条件（NEC）、弱能量条件（WEC）、强能量条件（SEC）及主导能量条件（DEC）均被满足。

质量-半径关系与观测比较：理论M-R曲线与21颗中子星观测数据比较显示，软模型（n_c=0.5, n_e=1, x_c=0.2）产生的较大半径与PSR J0030+0451（M≈1.34 M_⊙, R≈12.71 km）及PSR J0437-4715（M=1.44 M_⊙, R=13.02 km）等低质量、大半径系统一致；刚性模型（n_c=1, n_e=2, x_c=0.7）的高致密性则匹配PSR J1614?2230（M≈1.91 M_⊙, R≈9.7 km）、PSR J0348+0432（M≈2.01 M_⊙）及PSR J0952?0607（M≈2.35 M_⊙）等高质量致密脉冲星。热核爆发源（4U 1608-52、4U 1820-30）等中介系统位于两曲线之间，表明其内部结构可能介于两组参数之间。

最大质量分析：刚性复合构型（n_c=1, n_e=2, x_c=0.7）的最大质量随σ强烈单调增加。锐利过渡模型（ε=0.01）的M_max从σ=0.05时的0.168 M_⊙增至σ=0.50时的3.47 M_⊙；平滑过渡模型（ε=0.03）从0.12 M_⊙增至2.71 M_⊙。最大质量对应的中心密度随σ增加而降低约一个数量级，从~2.2×10¹⁶ g cm^?3降至~2.2×10¹⁵ g cm^?3。该模型在σ≥0.35（ε=0.01）或σ≥0.40（ε=0.03）时达到M_max≥2.0 M_⊙，满足PSR J1614?2230、PSR J0348+0432、PSR J0740+6620及PSR J0952?0607等观测约束。最小半径R_min在10.6–13.2 km范围内，与NICER测量 favored 的11.9–13.0 km及GW170817潮汐形变率约束一致，不支持半径低于~10.5 km的构型。与现有理论比较：自由简并中子气的经典结果M_max≈0.71 M_⊙远低于观测；均匀多方n=1的结果M_max≈0.56–0.82 M_⊙亦显著低于复合模型预测；Rhoades-Ruffini因果上限M_max≤3.2 M_⊙被锐利过渡模型在σ=0.50时的3.47 M_⊙所接近。与核EoS（DD2：M_max≈2.42 M_⊙, R_1.4≈13.2 km；GM1：M_max≈2.37 M_⊙）的比较 confirms 模型参数空间的物理合理性。

软核心模型（n_c=0.5, n_e=1, x_c=0.2）呈现不同的结构行为：M_max从σ=0.05时的0.117 M_⊙增至σ=0.50时的3.107 M_⊙（ε=0.01），但因果性限制σ≤0.3的物理允许范围内仅达M_max≈0.54–1.43 M_⊙与观测阈值尚有差距。

研究结论：复合多方过程通过核心刚度（n_c）、包层刚度（n_e）及核心尺寸分数（x_c）的灵活组合，能够覆盖观测到的中子星质量-半径多样性。MC方法不仅作为有效的数值求解器，更通过集合解提供了固有的不确定性量化与参数退化分析能力。该框架为利用NICER数据等观测约束进行MCMC参数估计、进一步约束中子星物态方程奠定了方法论基础，未来可扩展至包含旋转、更现实核物理及更尖锐相变的模型。