用宇宙线透视金字塔：\(\mu\) 子成像
2026 年 10 月
1 导读
每分钟大约有一个宇宙线 \(\mu\) 子穿过你的指甲盖大小的面积。它们诞生在十几千米高的大气层里，本来平均只能活 2.2 微秒，却能以接近光速穿透整个大气，甚至穿透上百米厚的岩石。\(\mu\) 子穿过物质时几乎走直线，被吸收的多少取决于路上“有多少物质”；被偏折的多少取决于物质的原子序数。于是，把一台能记录 \(\mu\) 子方向的探测器放在金字塔里、火山脚下或集装箱旁边，只要耐心计数，就能给这些庞然大物拍“X 光片”。本讲从狭义相对论、带电粒子与物质的相互作用和计数统计三块基础出发，一步一步算清楚：为什么 \(\mu\) 子能到达地面，它能穿透多厚的石头，金字塔里的空腔会让计数多出多少，需要曝光多少天才能确信空腔存在，以及怎样用散射角把铀和铅区分开来。
\(\mu\) 子是 1936–1937 年在宇宙线云室照片中发现的带电粒子，质量约为电子的 207 倍[1], [2]。现代测得其质量 \(m_\mu c^2={105.66\,}\mathrm{MeV}\)，静止时的平均寿命 \(\tau_\mu={2.197\,}\mathrm{\mu s}\)[3]。1940 年，Rossi 与 Hall 在美国科罗拉多州的 Echo Lake（海拔 3240 m）和 Denver（海拔 1616 m）比较不同动量的宇宙线 \(\mu\) 子在大气中的衰减（1941 年发表）。他们用铅吸收体把 \(\mu\) 子按穿透本领分成两组，发现穿透本领较弱（较慢）的一组衰变速率约为较强一组的 3 倍[4]，并由此推出 \(\mu\) 子静止寿命约为 \((2.4\pm0.3)\,\mu\mathrm{s}\)，这是运动时钟变慢的第一批定量证据之一[4], [5]。1962–1963 年，Frisch 与 Smith 在新罕布什尔州华盛顿山顶和海平面的剑桥分别计数，只选取在山顶速度为 \(0.9950c\) 至 \(0.9954c\) 的 \(\mu\) 子，测得时间膨胀因子为 \(8.8\pm0.8\)，与按其探测几何计算的有效值 \(8.4\pm2\) 相符，并把实验过程拍成了教学影片 Time Dilation: An Experiment with Mu-Mesons（1962）[5], [6]。第 7.7 节和例题 1 会用本讲推出的公式重新分析这个实验。
1955 年，澳大利亚物理学家 George 用宇宙线计数测量了雪山水电工程一条隧道上方岩层的面密度，这是已知最早的 \(\mu\) 子“透射测量”[7]。1960 年代，Alvarez 领导的团队把火花室放进埃及哈夫拉（Khafre）金字塔底部的贝尔佐尼室，按方向统计穿过金字塔的 \(\mu\) 子数，并与计算机模拟比较，1970 年在 Science 上发表结论：在所观测的范围内，没有与国王墓室体积相当的未知空腔[8], [9]。2007 年，田中宏幸（H. K. M. Tanaka）等人用核乳胶片拍到了日本浅间山火山口下方的密度结构[10]；2014 年他们用 \(\mu\) 子“延时摄影”观察到萨摩硫黄岛火山通道中岩浆柱的升降[11]。
2015 年起，国际合作 ScanPyramids 用三种不同技术的探测器（第 8.1 节介绍）观测胡夫金字塔（今高 139 m、底宽 230 m）：名古屋大学的核乳胶片放在王后墓室（位置 NE2）及其东侧的壁龛（位置 NE1），日本高能加速器研究机构（KEK）的闪烁体径迹仪放在王后墓室内，法国原子能委员会（CEA）的 Micromegas 气体探测器放在金字塔北面外侧。2017 年他们在 Nature 上报告，在大走廊（高 8.6 m、长 46.7 m）上方发现一个长度至少 30 m、截面与大走廊相当的“大空洞”（ScanPyramids Big Void）。核乳胶两个位置的统计显著性分别为 13.7\(\sigma\) 和 12.7\(\sigma\)；KEK 径迹仪测得超额 5%–11%，多数切片超过 10\(\sigma\)；CEA 两台望远镜合并为 5.8\(\sigma\)[9]。2023 年同一合作又精确测定了北面“人字形”石块后方一条长约 9 m、截面约 \(2\,\mathrm{m}\times2\,\mathrm{m}\) 的走廊状结构[12]（图 1）。
另一类方法利用 \(\mu\) 子在物质中的偏折。2003 年，美国洛斯阿拉莫斯国家实验室的 Borozdin 等人在 Nature 上提出并演示：测量每个宇宙线 \(\mu\) 子穿过物体前后的方向，由多次库仑散射的强弱可以为高原子序数材料成像，用于在货柜中寻找铀、钚等核材料[13]。2011 年福岛第一核电站事故后，\(\mu\) 子成像被用来从反应堆厂房外部估计熔毁燃料的位置：2012 年的模拟研究论证了散射法的可行性[14]，KEK 等单位用透射法观测了 1 号机组，估计压力容器内残留燃料的多少[15]。
学习目标
学完本讲，你应当能够：
用固有时与指数衰变律计算运动 \(\mu\) 子的存活概率，并计入沿途能量损失对时间膨胀的影响；
由“一次擦肩而过的冲量”估计带电粒子的电离能量损失，理解 Bethe 公式的标度、最小电离值约 \(2\,\mathrm{MeV\,cm^2/g}\) 与“不透明度”\(X=\int\rho\,\dd l\) 的意义，并由 \(-\dd E/\dd X=a+bE\) 求出穿透厚度 \(X\) 所需的最小能量；
由积分能谱的局部斜率推出透射计数对不透明度的灵敏度 \(\dd\ln\Phi/\dd\ln X\)，估计空腔造成的超额与达到给定显著性所需的曝光时间，并理解系统误差设定的“显著性天花板”；
由随机行走与卢瑟福散射理解多次散射角 \(\psi_0\propto(1/p)\sqrt{x/X_0}\)，用它区分铁、铅、铀，并估计动量未知时所需的 \(\mu\) 子数；
完成“天顶角分布”实验数据的处理：几何接收度、线性化、加权最小二乘、不确定度与系统误差。
力学与电磁学：冲量、动量与能量、库仑定律、电场的高斯定理；
近代物理：知道狭义相对论的时间膨胀与 \(E=mc^2\) 即可，第 3 节会补齐本讲用到的全部公式；
数学：泰勒展开、定积分换元、对数求导、简单的概率（期望与方差）。
本讲新引入的工具（相对论运动学、固有时、衰变律、阻止本领与 Bethe 公式、连续减速近似、泊松分布与显著性、多次散射与辐射长度、强度与立体角、加权最小二乘）都在首次使用前给出定义和推导。
建议学时
共 2 学时（约 90 分钟）：前置模块 30 分钟；宇宙线 \(\mu\) 子与 Frisch–Smith 实验 15 分钟；透射成像 25 分钟；散射成像 15 分钟；小结 5 分钟。例题与实验题供课后练习，附录供深入阅读。
推导路线
时间膨胀 + 衰变律 \(\Rightarrow\) \(\mu\) 子能到达地面；玻尔冲量估计 \(\Rightarrow\) Bethe 公式 \(\Rightarrow\) \(-\dd E/\dd X=a+bE\) \(\Rightarrow\) 最小穿透能量 \(E_{\min}(X)\)；能谱 \(\Rightarrow\) 透射强度 \(\Phi(X)\) \(\Rightarrow\) 空腔超额 \(\varepsilon\approx\eta'\ell/L\)；泊松统计 \(\Rightarrow\) 显著性 \(n_\sigma=\varepsilon\sqrt N\) \(\Rightarrow\) 曝光时间。另一条线：卢瑟福散射 + 随机行走 \(\Rightarrow\) \(\psi_0\propto\sqrt{x/X_0}/p\) \(\Rightarrow\) 材料识别。
2 名词与符号
| 符号 | 含义 | 常用单位 |
|---|---|---|
| \(\beta=v/c\), \(\gamma\) | 速度与光速之比；洛伦兹因子 \(1/\sqrt{1-\beta^2}\) | 1 |
| \(E\), \(p\), \(m\) | 总能量（含静能）、动量、静质量 | GeV, GeV/\(c\), GeV/\(c^2\) |
| \(\tau_\mu\), \(\tau\) | \(\mu\) 子静止平均寿命 2.197 μs；运动粒子经历的固有时 | \(\mu\)s |
| \(P\) | 存活概率 | 1 |
| \(\rho\), \(X\) | 密度；不透明度（面密度）\(X=\int\rho\,\dd l\) | g/cm\(^3\), g/cm\(^2\) |
| \(-\dd E/\dd X\) | 质量阻止本领（单位面密度的平均能量损失） | MeV cm\(^2\)/g |
| \(Z\), \(A\) | 原子序数、摩尔质量（数值上等于质量数） | 1, g/mol |
| \(\bar I\) | 平均激发能（Bethe 公式中） | eV |
| \(b\) | 碰撞参数（第 4.2、6 节与附录 A、B） | m |
| \(K\), \(r_{\rm e}\) | Bethe 常数 \(4\pi N_{\rm A}r_{\rm e}^2m_{\rm e}c^2\)；经典电子半径 | MeV cm\(^2\)/mol, m |
| \(W_{\max}\), \(\delta\) | 单次碰撞最大能量转移；密度效应修正 | MeV, 1 |
| \(a\), \(b\) | 能损参数：\(-\dd E/\dd X=a+bE\)（此处 \(b\) 与碰撞参数无关） | MeV cm\(^2\)/g, cm\(^2\)/g |
| \(E_{\min}(X)\) | 恰好穿透不透明度 \(X\) 所需的最小能量 | GeV |
| \(\theta\), \(\varphi\), \(\Omega\) | 天顶角、方位角、立体角 | rad, sr |
| \(\tan\theta_x\), \(\tan\theta_y\) | 径迹在东西、南北竖直平面上投影倾角的正切 | 1 |
| \(\lambda_{\rm d}\), \(\lambda_I\), \(\Lambda\) | 衰变长度；核相互作用长度；核子流衰减长度 | m；g/cm\(^2\) |
| \(E_\pi^*\), \(E_K^*\) | \(\pi\)、K 介子“衰变/碰撞”分界能量 | GeV |
| \(I(p,\theta)\) | 微分强度：单位时间、单位垂直面积、单位立体角、单位动量的粒子数 | m\(^{-2}\)s\(^{-1}\)sr\(^{-1}\)(GeV/\(c\))\(^{-1}\) |
| \(I(>p,\theta)\) | 积分强度（动量大于 \(p\)） | m\(^{-2}\)s\(^{-1}\)sr\(^{-1}\) |
| \(\Phi(X,\theta)\) | 透射强度 \(I(>E_{\min}(X),\theta)\) | m\(^{-2}\)s\(^{-1}\)sr\(^{-1}\) |
| \(J\) | 穿过平面的通量 | m\(^{-2}\)s\(^{-1}\) |
| \(S_{\rm eff}\), \(\Delta\Omega\), \(G\) | 有效面积；立体角窗口；接收度 \(\int S_{\rm eff}\dd\Omega\) | m\(^2\), sr, m\(^2\) sr |
| \(g(u)\) | \(u\ee^u/(\ee^u-1)\)，\(u=bX\) | 1 |
| \(\eta\), \(\eta'\) | 积分谱局部斜率 \(-\dd\ln I/\dd\ln p\)；灵敏度 \(-\dd\ln\Phi/\dd\ln X\) | 1 |
| \(L\), \(\ell\) | 视线方向岩石总长度；视线穿过空腔的长度 | m |
| \(\varepsilon\) | 空腔造成的相对计数超额 | 1 |
| \(N\), \(\mu_N\) | 计数；计数的期望值 | 1 |
| \(n_\sigma\) | 显著性（偏离几个标准差） | 1 |
| \(\sigma_{\rm sys}\) | 无目标预言的相对系统误差 | 1 |
| \(w\), \(\Delta\) | 空洞角宽；分析窗口宽度 | 1（以 \(\tan\theta\) 计） |
| \(S\), \(T\) | 探测器面积、曝光时间（\(T_{\rm a}\) 为气温） | m\(^2\), s 或天 |
| \(X_0\), \(E_{\rm s}\) | 辐射长度；Rossi–Greisen 能量 21.2 MeV | g/cm\(^2\) 或 cm, MeV |
| \(\psi_x\), \(\psi_0\) | 投影散射角；其均方根（Highland 公式） | mrad |
| \(\kappa\) | 峰度比 \(\avg{\psi^4}/\avg{\psi^2}^2\) | 1 |
| \(p_{\rm eff}\), \(p_{\rm c}\) | 有效动量 \(\avg{(\beta p)^{-2}}^{-1/2}\)；截断动量 | GeV/\(c\) |
| \(Q_1\), \(Q_2\), \(R\) | 两种假设预言的均方投影角及其比值 \(Q_1/Q_2\) | rad\(^2\), 1 |
| \(M\) | 入射重粒子质量（Bethe 公式）；\(M_{\rm a}\)、\(R_{\rm g}\) 为空气摩尔质量与气体常数 | — |
| \(n\) | 天顶角分布指数，\(I\propto\cos^n\theta\) | 1 |
| \(h\), \(H\) | 海拔高度；大气标高 | km |
| \(d\), \(D\) | 实验题中探测平面边长、两平面间距 | m |
| —— | 例题中另有局部定义的符号以题面为准（如例题 2 的 \(t_{\rm v}\)、\(H_{\rm p}\)、\(B_{\rm p}\)、\(b_{\rm d}\)） | — |
能量单位统一用电子伏特：\({1\,}\mathrm{eV}={1.602\times 10^{-19}\,}\mathrm{J}\)，\({1\,}\mathrm{GeV}=10^9\,\mathrm{eV}\)。动量用 GeV/\(c\)，质量用 GeV/\(c^2\)，这样 \(pc\)、\(mc^2\) 与 \(E\) 都以 GeV 计。“米水当量”（m.w.e.）是不透明度的工程单位：1 m 厚的水的面密度，1 m.w.e. \(={100\,}\mathrm{g/cm^2}\)。
3 前置模块一：相对论运动学与衰变
3.1 洛伦兹因子、能量与动量
静质量为 \(m\)、速度为 \(v=\beta c\) 的粒子，其能量与动量为 \begin{align*} E&=\gamma mc^2,\qquad p=\gamma mv, \\ \gamma&=\frac{1}{\sqrt{1-\beta^2}} . \tag{3.1} \end{align*} 由两式消去 \(\gamma\) 得能量–动量关系 \[\tag{3.2} E^2=(pc)^2+(mc^2)^2 .\]
式 (3.1) 两式相除得 \(\beta=pc/E\)；相乘与相除可得 \[\tag{3.3} \beta\gamma=\frac{p}{mc},\qquad \gamma=\frac{E}{mc^2}.\] 验证式 (3.2)：\(E^2-(pc)^2=\gamma^2m^2c^4(1-\beta^2)=m^2c^4\)。当 \(pc\gg mc^2\)（“极端相对论”）时 \(E\approx pc\)，\(\beta\approx1-\tfrac12(mc^2/E)^2\)。本讲的 \(\mu\) 子能量多在 GeV 以上，\(m_\mu c^2={0.1057\,}\mathrm{GeV}\)，所以常可不区分 \(E\) 与 \(pc\)；只在 \(E\lesssim{1\,}\mathrm{GeV}\) 时需要用精确式。
海平面 \(\mu\) 子的平均能量约 4 GeV[16]。此时 \(\gamma=4/0.10566=37.9\)，\(\beta=0.99965\)。速度与光速只差万分之三点五，但 \(\gamma\) 已经是 38。
3.2 时间膨胀与固有时
粒子随身携带的钟所走过的时间叫固有时，记作 \(\tau\)。若在实验室（地面）参考系中粒子以速度 \(v(t)\) 运动，实验室时间 \(\dd t\) 内粒子的钟只走过 \[\tag{3.4} \dd\tau=\frac{\dd t}{\gamma(t)}=\sqrt{1-\beta^2}\,\dd t .\] 这就是时间膨胀：在实验室看来，运动的钟走得慢了 \(\gamma\) 倍。
狭义相对论中，时间膨胀是对“相对实验室匀速运动的钟”导出的。若钟在加速，可以把运动切成很多小段，每一小段内速度近似不变，用式 (3.4) 得到这一小段的固有时，再求和（积分）。这样做的前提是钟本身不因加速而损坏：\(\mu\) 子的“钟”是它的衰变概率，与加速度无关：在储存环中做高速圆周运动（加速度极大）的 \(\mu\) 子，其寿命仍只由 \(\gamma\) 决定，这已被实验证实[17]。本讲中 \(\mu\) 子在大气或岩石中逐渐减速，正需要这种逐段积分。
\(\mu\) 子沿竖直方向飞过距离 \(\dd x\) 用时 \(\dd t=\dd x/v\)，于是 \[\tag{3.5} \dd\tau=\frac{\dd x}{\gamma v}=\frac{\dd x}{\beta\gamma c}=\frac{m\,\dd x}{p}.\] 最后一步用了式 (3.3)。式 (3.5) 说明：每走单位路程，粒子的钟走过的时间与动量成反比；动量越大，“老得越慢”。
3.3 指数衰变律
一个 \(\mu\) 子在它自己的钟走过 \(\dd\tau\) 的时间内衰变的概率为 \(\dd\tau/\tau_\mu\)，与它已经“活了多久”无关，也与它的速度、所处环境无关。\(\tau_\mu\) 是一个常数，叫平均寿命。这是对大量粒子衰变实验的总结（量子力学给出其解释）。
设固有时 \(\tau\) 时仍存活的概率为 \(P(\tau)\)。在接下来的 \(\dd\tau\) 内，存活需要“到 \(\tau\) 时还活着”且“这一小段内不衰变”，两者独立，故 \begin{align*} &P(\tau+\dd\tau)=P(\tau)\Big(1-\frac{\dd\tau}{\tau_\mu}\Big)\\ \Longrightarrow\quad &\dv{P}{\tau}=-\frac{P}{\tau_\mu}. \end{align*} 以 \(P(0)=1\) 积分得 \[\tag{3.6} P=\ee^{-\tau/\tau_\mu}.\] 再用式 (3.5) 把固有时换成路程：若动量 \(p\) 在飞行中不变，飞过 \(L\) 后 \(\tau=mL/p\)， \begin{align*} P&=\exp\Big(-\frac{L}{\lambda_{\rm d}}\Big), \\ \lambda_{\rm d}&\equiv\beta\gamma c\tau_\mu=\frac{p}{mc}\,c\tau_\mu . \tag{3.7} \end{align*} \(\lambda_{\rm d}\) 叫衰变长度。对 \(\mu\) 子，\(c\tau_\mu={658.6\,}\mathrm{m}\)。若动量在飞行中变化，则用 \(\tau=\int m\,\dd x/p(x)\) 代入式 (3.6)。
设 \(\mu\) 子在 15 km 高处产生（第 7.2 节说明这个数从何而来），竖直向下，能量 4 GeV，先忽略能量损失。
不考虑时间膨胀：飞行时间 \(t=L/v\approx{15\,}\mathrm{km}/c={50.0\,}\mathrm{\mu s}=22.8\,\tau_\mu\)，存活概率 \(\ee^{-22.8}\approx1.3\times10^{-10}\)；
考虑时间膨胀：\(\lambda_{\rm d}=37.9\times{658.6\,}\mathrm{m}={24.9\,}\mathrm{km}\)，\(P=\ee^{-15/24.9}=0.55\)。
两者相差 9 个数量级。实际上地面每平方米每秒能接收到上百个 \(\mu\) 子（第 7.4 节），没有时间膨胀就无法解释。用 \(\mu\) 子自己的参考系看，同一件事表现为“大气层厚度被洛伦兹收缩为 \({15\,}\mathrm{km}/\gamma={0.40\,}\mathrm{km}\)”，\(\mu\) 子在寿命内足以飞过。图 2 给出不同到达能量下的三种计算结果。
4 前置模块二：带电粒子穿过物质
4.1 面密度与不透明度
带电粒子穿过物质时，主要和物质中的电子与原子核打交道。穿过的“物质多少”用沿路径的质量面密度衡量：
沿粒子路径对密度积分 \[\tag{4.1} X=\int_{\text{路径}}\rho(l)\,\dd l\] 叫不透明度（opacity）或面密度，单位 g/cm\(^2\)。它等于以路径为轴、截面为 1 cm2 的细柱体内的物质质量。例如 90 m 厚、密度 2.2 g/cm3 的石灰岩，\(X=2.2\times9000={1.98\times 10^{4}\,}\mathrm{g/cm^2}\)，即 198 米水当量。
地下与岩石中的 \(\mu\) 子计算常用 PDG 约定的“标准岩石”：\(Z=11\)，\(A=22\)（\(Z/A=0.5\)），\(\rho={2.65\,}\mathrm{g/cm^3}\)，平均激发能 \(\bar I={136.4\,}\mathrm{eV}\)（见第 4.3 节），辐射长度（第 6.3 节定义）26.54 g/cm2[3], [18]。真实岩石的成分与密度不同（金字塔的石灰岩约 2.2 g/cm3），但以面密度计的能损与标准岩石相差不大，所以本讲用石灰岩的密度算不透明度，能损参数借用标准岩石。
为什么用 \(X\) 而不用长度？下面会看到，单位路程的电离能量损失正比于电子数密度 \(n_{\rm e}=\rho N_{\rm A}Z/A\)（\(N_{\rm A}\) 为阿伏伽德罗常数），而绝大多数物质的 \(Z/A\) 都接近 \(1/2\)（氢除外）。于是“单位面密度”的能量损失几乎与材料无关，岩石、水、铁的差别都在 30% 以内。用 \(X\) 作自变量，可以把“穿过多少物质”和“物质是什么”大致分开。
4.2 一次擦肩而过：玻尔的冲量估计
一个电荷 \(ze\)、速度 \(v\) 的重粒子（质量 \(M\gg m_{\rm e}\)）从一个电子旁边飞过，最近距离（碰撞参数）为 \(b\)。在碰撞的短暂时间内，(i) 重粒子几乎不偏转，速度不变；(ii) 电子来不及移动，可视为静止的自由电子（原子束缚的作用在下文用 \(b_{\max}\) 计入）。
取重粒子沿 \(x\) 轴运动，\(t=0\) 时最靠近电子。电子受的力沿运动方向的分量在 \(t\) 与 \(-t\) 时大小相等、方向相反，积分为零；只剩垂直分量。这个分量由库仑定律给出（非相对论情形）： \[F_\perp(t)=\frac{ze^2}{4\pi\varepsilon_0}\,\frac{b}{(b^2+v^2t^2)^{3/2}} .\] 对时间积分（令 \(u=vt/b\)，\(\int_{-\infty}^{\infty}(1+u^2)^{-3/2}\dd u=2\)）得冲量 \[\tag{4.2} \Delta p_\perp=\int_{-\infty}^{\infty}F_\perp\,\dd t=\frac{ze^2}{4\pi\varepsilon_0}\,\frac{2}{bv}.\] 一个更快捷的看法分两步。(i) 在某一时刻，以重粒子轨迹为轴作一个半径 \(b\)、长度无限的圆柱面，由高斯定理，穿过侧面的电通量 \(\int E_\perp(x)\,2\pi b\,\dd x=ze/\varepsilon_0\)（\(x\) 为沿轴坐标）。(ii) 匀速运动电荷的场随电荷整体平移：\(E_\perp(x,t)=E_\perp(x-vt)\)。所以固定在 \(x=0\) 处的电子感受到的场随时间的变化 \(E_\perp(-vt)\)，与某一时刻场沿轴的分布形状相同，\(\int E_\perp\,\dd t=\int E_\perp\,\dd x/v\)。两式合起来，\(e\int E_\perp\dd t=2ze^2/(4\pi\varepsilon_0bv)\)。这个论证只用到高斯定理和“场随电荷平移”，不依赖场的具体形状，因此对高速（相对论）粒子也成立；相对论下垂直场的峰值增大、作用时间缩短，两者恰好抵消（运动电荷电场的这一结论本讲不推导）。
电子获得的动能为 \(\Delta E(b)=\Delta p_\perp^2/(2m_{\rm e})\)。引入经典电子半径 \(r_{\rm e}=e^2/(4\pi\varepsilon_0m_{\rm e}c^2)={2.818\times 10^{-15}\,}\mathrm{m}\)，写成 \[\tag{4.3} \Delta E(b)=\frac{2z^2r_{\rm e}^2\,m_{\rm e}c^2}{\beta^2}\,\frac{1}{b^2}.\]
重粒子前进 \(\dd x\) 时，碰撞参数在 \(b\) 到 \(b+\dd b\) 之间的电子位于体积为 \(2\pi b\,\dd b\,\dd x\) 的薄圆筒内，个数为 \(n_{\rm e}\,2\pi b\,\dd b\,\dd x\)。把每个电子获得的能量相加： \begin{align*} -\dv{E}{x}&=\int_{b_{\min}}^{b_{\max}}n_{\rm e}\,\Delta E(b)\,2\pi b\,\dd b\\ &=\frac{4\pi n_{\rm e}z^2r_{\rm e}^2m_{\rm e}c^2}{\beta^2}\ln\frac{b_{\max}}{b_{\min}} . \end{align*} 积分在 \(b\to0\) 和 \(b\to\infty\) 都发散（对数发散），所以必须有物理上的截断：
\(b_{\min}\)：\(b\) 太小时式 (4.3) 给出的能量超过运动学允许的最大值 \(W_{\max}\approx2m_{\rm e}c^2\beta^2\gamma^2\)（正碰时电子能得到的最大能量），冲量近似失效；
\(b_{\max}\)：碰撞持续时间约为 \(b/(\gamma v)\)。若它长于电子在原子中的典型运动周期 \(1/\bar\omega\)，电子会“绝热地”跟随，几乎得不到能量。所以 \(b_{\max}\approx\gamma v/\bar\omega\)。
除以 \(\rho\) 化为质量阻止本领，并代入 \(n_{\rm e}=\rho N_{\rm A}Z/A\)： \begin{align*} -\frac{\dd E}{\dd X}&=K\,z^2\,\frac{Z}{A}\,\frac{1}{\beta^2} \\ &\quad\times\ln\frac{b_{\max}}{b_{\min}}, \\ K&\equiv4\pi N_{\rm A}r_{\rm e}^2m_{\rm e}c^2 \\ &=0.307\ \mathrm{MeV\,cm^2/mol}. \tag{4.4} \end{align*}
式 (4.4) 是 1913 年玻尔给出的经典结果的核心[19]。它已经告诉我们三件事：能损只依赖于粒子的速度（而非质量）；在 \(\beta\to1\) 之前按 \(1/\beta^2\) 下降；对材料的依赖只通过 \(Z/A\)（以及对数里的弱依赖）。
4.3 Bethe 公式与最小电离
量子力学的完整计算（Bethe，1930 年代）给出重带电粒子的平均质量阻止本领[3]： \begin{align*} &-\frac{\dd E}{\dd X}=K z^2\frac{Z}{A}\frac{1}{\beta^2} \\ &\times\bigg[\frac12\ln\frac{2m_{\rm e}c^2\beta^2\gamma^2W_{\max}}{\bar I^2} \\ &\qquad-\beta^2-\frac{\delta(\beta\gamma)}{2}\bigg]. \tag{4.5} \end{align*} \(\bar I\) 是物质的平均激发能（原子中电子的平均“束缚能标”，约 \({10\,}\mathrm{eV}\times Z\)，例如标准岩石 136.4 eV、铅 823 eV）；\(W_{\max}\) 为单次碰撞的最大能量转移（精确式为 \begin{align*} W_{\max}&=\frac{2m_{\rm e}c^2\beta^2\gamma^2}{1+2\gamma m_{\rm e}/M+(m_{\rm e}/M)^2}, \end{align*} \(M\) 为入射粒子质量；当 \(2\gamma m_{\rm e}/M\ll1\)、对 \(\mu\) 子即 \(\gamma\ll100\) 时 \(\approx2m_{\rm e}c^2\beta^2\gamma^2\)。\(\gamma\sim100\) 时 \(W_{\max}\) 约减半，但它只在对数里，对能损的影响约 3%）；\(\delta\) 是密度效应修正：高速粒子的电场使介质极化，屏蔽了远处的电子，削弱了对数项的增长。适用范围约 \(0.1\lesssim\beta\gamma\lesssim1000\)，精度为百分之几[3]。
对照式 (4.4)：取 \(b_{\max}=\gamma v/\bar\omega\)，量子力学中 \(b_{\min}\) 由电子在粒子参考系中的德布罗意波长决定，\(b_{\min}\approx\hbar/(2\gamma m_{\rm e}v)\)，于是 \begin{align*} \ln(b_{\max}/b_{\min})&=\ln(2\gamma^2m_{\rm e}v^2/\hbar\bar\omega)\\&=\ln(2m_{\rm e}c^2\beta^2\gamma^2/\bar I) \end{align*} （\(\hbar\bar\omega\equiv\bar I\)），正是 \(W_{\max}\approx2m_{\rm e}c^2\beta^2\gamma^2\) 时式 (4.5) 方括号中的主项。附录 A 给出这一对应的细节。
\(\beta^{-2}\) 随速度增大而下降，对数项随 \(\gamma\) 缓慢上升，两者之间有一个宽而平的极小值，位置约在 \(\beta\gamma\approx3\)–\(3.5\)[3]。对标准岩石（\(Z/A=0.5\)，\(\bar I={136.4\,}\mathrm{eV}\)），取 \(\beta\gamma=3.5\)（\(\beta^2=0.925\)），忽略 \(\delta\)： \begin{align*} \ln\frac{2m_{\rm e}c^2\beta^2\gamma^2}{\bar I}&=\ln\frac{2\times0.511\times10^6\times12.25}{136.4}\\ &=11.4, \end{align*} \begin{align*} -\frac{\dd E}{\dd X}&\approx\frac{0.307\times0.5}{0.925}\,(11.4-0.93)\\ &=1.74\ \mathrm{MeV\,cm^2/g}. \end{align*} PDG 表值为 1.69 MeV cm2/g，位于 \(\beta\gamma=3.7\)[18]。手算误差 3%。
宇宙线 \(\mu\) 子的 \(\beta\gamma\) 多在 10–100 之间，处在极小值右侧的平坦区，能损率接近极小值：\(-\dd E/\dd X\approx{2\,}\mathrm{MeV\cdot cm^2/g}\)（准确值随材料在 1.1–2.0 之间，图 3）。这样的粒子叫最小电离粒子（minimum ionizing particle，MIP）。以后我们把岩石中的电离能损近似为常数 \(a\approx{2.2\,}\mathrm{MeV\cdot cm^2/g}\)（略高于极小值，计入了对数上升）。
式 (4.5) 给的是平均能损。单次碰撞偶尔把很大能量交给一个电子（产生“\(\delta\) 电子”），所以一段薄物质中的实际能损分布有一条长尾，最可几值低于平均值。穿过很厚的物质（如上百米岩石）时，大量碰撞的总和的相对涨落变小，用平均能损描述是合理的，但仍有百分之几到十几的“射程歧离”。本讲在需要精确数值时直接查 PDG 表。
4.4 高能 \(\mu\) 子的辐射损失
能量高于几百 GeV 时，三种“辐射”过程变得重要：轫致辐射（\(\mu\) 子在原子核电场中被急剧偏转而辐射光子）、直接电子对产生（在核电场中直接产生一对电子与正电子）、光核作用（\(\mu\) 子通过交换一个虚光子打碎原子核）。这三种过程损失的能量都与 \(\mu\) 子自身的能量成正比（粗略地说，每个过程都以一定概率带走 \(\mu\) 子能量的一个固定比例）。于是总能损可以写成[16] \[\tag{4.6} -\dv{E}{X}=a+bE ,\] \(a\) 是电离损失，\(b\) 是三种辐射过程的“分数能损率”之和。（注意：这里的 \(b\) 与第 4.2 节的碰撞参数 \(b\) 无关，本讲中二者不会在同一处出现。）两者都随能量缓慢变化：对标准岩石，\(E={10\,}\mathrm{GeV}\) 到 10 TeV，\(a\) 从 2.17 增至 2.93 MeV cm2/g，\(b\) 从 1.9 × 10−6 增至 4.35 × 10−6 cm2/g[16]。两项相等的能量 \[\tag{4.7} E_{\mu\rm c}=\frac{a}{b}\] 叫 \(\mu\) 子的临界能量，标准岩石约 500 GeV–700 GeV[16], [18]。本讲取常数 \(a={2.2\,}\mathrm{MeV\cdot cm^2/g}\)、\(b={4.0\times 10^{-6}\,}\mathrm{cm^2/g}\)（\(E_{\mu\rm c}={550\,}\mathrm{GeV}\)）。用这对常数由式 (4.8) 算出的穿透（动）能，与 PDG 射程表相比：10 GeV–1 TeV 内相差在 \(\pm9\%\) 以内，5 GeV 处偏高 14%，1 GeV 处偏高约 22%，2 TeV 处偏低 11%（图 4a）。本讲关心的金字塔（几十 GeV）与火山（几百 GeV）都在精度较好的区间。
4.5 连续减速近似与最小穿透能量
把 \(\mu\) 子的能量损失看成沿路径连续、确定的过程，每一段都按平均能损 \(a+bE\) 减速；忽略能损的涨落（射程歧离）、大角散射造成的路径弯曲，也忽略 \(\mu\) 子在岩石中衰变（第 3.2 节的衰变长度在 GeV 以上为几十千米，远大于岩石厚度）。
由式 (4.6) 分离变量：\(\dd X=-\dd E/(a+bE)\)。设 \(\mu\) 子以能量 \(E_0\) 进入岩石，恰好在走完 \(X\) 时停下。积分下限取 0 意味着这里的 \(E\) 应理解为动能；在 \(E\gg m_\mu c^2\) 时动能、总能量与 \(pc\) 不必区分。 \[X=\int_0^{E_0}\frac{\dd E}{a+bE}=\frac1b\ln\Big(1+\frac{bE_0}{a}\Big).\] 反解出“恰好穿透”的入射能量，记为 \(E_{\min}\)： \[\tag{4.8} E_{\min}(X)=\frac ab\left(\ee^{bX}-1\right).\] 能量高于 \(E_{\min}\) 的 \(\mu\) 子能穿透，低于的被“吃掉”。
\(bX\ll1\)（只有电离损失重要）：\(\ee^{bX}-1\approx bX+\tfrac12(bX)^2\)，\(E_{\min}\approx aX\,(1+\tfrac12bX)\)。金字塔 \(X\approx{2\times 10^{4}\,}\mathrm{g/cm^2}\)，\(bX=0.08\)， \begin{align*} E_{\min}&\approx2.2\ \mathrm{MeV\,cm^2/g}\\&\quad\times2\times10^{4}\ \mathrm{g/cm^2}\times1.04\\&=46\ \mathrm{GeV}. \end{align*}
\(bX\gtrsim1\)（辐射损失主导）：\(E_{\min}\) 随 \(X\) 指数增长。火山 \(X\approx{2\times 10^{5}\,}\mathrm{g/cm^2}\)，\(bX=0.8\)，\(E_{\min}\approx{550\,}\mathrm{GeV}\times1.23={670\,}\mathrm{GeV}\)；PDG 射程表给出 728 GeV。
量纲：\(a/b\) 是能量，\(bX\) 无量纲。
第 41 届复赛“质子断层扫描”利用的是质子射程对能量的依赖。本讲反过来用：入射粒子的能量不受我们控制（宇宙线谱是给定的），能量损失只用来决定“哪些 \(\mu\) 子能穿过来”，测量的是穿透数目，而不是单个粒子的射程。
5 前置模块三：计数统计
5.1 泊松分布
把曝光时间 \(T\) 切成 \(n_{\rm s}\) 个很短的小段，每段内到达一个 \(\mu\) 子的概率为 \(q\ll1\)，到达两个的概率可忽略，各段相互独立（宇宙线 \(\mu\) 子的到达彼此无关）。第 \(k\) 段的计数 \(N_k\) 只取 0 或 1： \begin{align*} &\avg{N_k}=q,\qquad \avg{N_k^2}=q,\\ &\operatorname{Var}N_k=\avg{N_k^2}-\avg{N_k}^2=q(1-q)\approx q . \end{align*} 总计数 \(N=\sum_kN_k\)。期望相加：\(\mu_N\equiv\avg{N}=n_{\rm s}q\)。独立随机量之和的方差等于方差之和（交叉项 \(\avg{(N_j-q)(N_k-q)}=0\)）： \begin{align*} \operatorname{Var}N&=n_{\rm s}q(1-q) \\ &\to n_{\rm s}q=\mu_N, \\ \sigma_N&=\sqrt{\mu_N}. \tag{5.1} \end{align*} 这里取了 \(n_{\rm s}\to\infty\)、\(q\to0\)、\(n_{\rm s}q\) 固定的极限。在同样的极限下，\(N\) 的分布是泊松分布 \(\Pr(N)=\mu_N^N\ee^{-\mu_N}/N!\)（附录 C）；\(\mu_N\gtrsim20\) 时它很接近均值 \(\mu_N\)、标准差 \(\sqrt{\mu_N}\) 的高斯分布。
计到 \(N\) 个事例，其统计误差约 \(\sqrt N\)，相对误差 \(1/\sqrt N\)。想把相对误差减小一半，计数（即曝光时间或面积）要增大 4 倍。
5.2 显著性
设某方向“没有空腔”时期望计数为 \(\mu_0\)（由模拟或参照区域预言，假定其误差可忽略），实测为 \(N\)。定义 \[\tag{5.2} n_\sigma=\frac{N-\mu_0}{\sqrt{\mu_0}} ,\] 即实测偏离预言几个标准差。若没有空腔，\(N\) 近似为均值 \(\mu_0\)、标准差 \(\sqrt{\mu_0}\) 的高斯分布，纯属涨落而超过 \(n_\sigma\) 的概率（单侧）为 \(n_\sigma=2\)：2.3%；\(n_\sigma=3\)：\(1.3\times10^{-3}\)；\(n_\sigma=5\)：\(2.9\times10^{-7}\)。粒子物理习惯以 5\(\sigma\) 作为“发现”的门槛。
(i) 若 \(\mu_0\) 不是可靠的理论预言，而来自另一个计数 \(N_0\)（如参照方向），则差值 \(N-N_0\) 的方差为 \(N+N_0\)（独立量方差相加），显著性变为 \((N-N_0)/\sqrt{N+N_0}\)。(ii) 在一张图像的成百上千个像素里寻找“最显著”的一个，即使没有空腔，也很可能遇到 3\(\sigma\) 的涨落，这叫“到处找”效应（look-elsewhere effect）；这正是门槛定在 5\(\sigma\) 而不是 3\(\sigma\) 的原因之一。
5.3 用样本估计方差及其误差
散射成像要由 \(N\) 个散射角 \(\psi_1,\dots,\psi_N\)（均值为零）估计分布的方差 \(\avg{\psi^2}\)。最自然的估计量是 \(\hat s^2=\frac1N\sum\psi_i^2\)。
\(\hat s^2\) 是 \(N\) 个独立随机量 \(\psi_i^2\) 的平均，其期望 \(\avg{\psi^2}\)，方差为单个 \(\psi^2\) 的方差除以 \(N\)： \begin{align*} \operatorname{Var}\hat s^2&=\frac{\avg{\psi^4}-\avg{\psi^2}^2}{N}=\frac{(\kappa-1)\avg{\psi^2}^2}{N},\\ \kappa&\equiv\frac{\avg{\psi^4}}{\avg{\psi^2}^2}. \end{align*} 所以 \[\tag{5.3} \frac{\sigma(\hat s^2)}{\avg{\psi^2}}=\sqrt{\frac{\kappa-1}{N}} .\] 高斯分布 \(\kappa=3\)，相对误差 \(\sqrt{2/N}\)。若分布有长尾（\(\kappa\) 大），同样的 \(N\) 给出的估计更差。附录 C 证明：宽度不同的高斯分布按权重混合后，\(\kappa=3\avg{\sigma^4}/\avg{\sigma^2}^2\ge3\)。
6 前置模块四：多次库仑散射
6.1 单次小角散射
\(\mu\) 子从原子核旁飞过时也受到库仑力。核的质量远大于 \(\mu\) 子，核几乎不动，\(\mu\) 子被偏折。仿照式 (4.2)，把电子换成电荷 \(Ze\) 的原子核、把“受冲量者”换成 \(\mu\) 子：\(\mu\) 子获得横向冲量 \(\Delta p_\perp=2zZe^2/(4\pi\varepsilon_0bv)\)，偏转角 \begin{align*} \psi(b)&=\frac{\Delta p_\perp}{p}=\frac{2zZ\,r_{\rm e}m_{\rm e}c^2}{pv\,b} \\ &=\frac{2zZ\,r_{\rm e}}{b}\cdot\frac{m_{\rm e}c^2}{\beta\,pc}. \tag{6.1} \end{align*} 这是卢瑟福散射的小角极限。偏转角与 \(pv=\beta pc\) 成反比（非相对论时 \(pv\) 等于两倍动能，极端相对论时 \(pv\approx E\)）：能量越高的粒子越“硬”，越难被偏折。
6.2 随机行走与高斯分布
厚度 \(x\) 的物质中，\(\mu\) 子会与成千上万个原子核发生小角碰撞，每次偏折的方向（在垂直于运动方向的平面内）是随机的。把总偏转角在 \(x\) 方向的投影写成 \(\psi_x=\sum_j\psi_{x,j}\)。各次碰撞相互独立且 \(\avg{\psi_{x,j}}=0\)，所以和的方差等于方差之和，正如式 (5.1) 的推导： \[\avg{\psi_x^2}=\sum_j\avg{\psi_{x,j}^2}\propto\text{碰撞次数}\propto x .\] 这就是随机行走：均方偏离随步数线性增长，均方根随 \(\sqrt x\) 增长。由中心极限定理（大量独立小量之和近似服从高斯分布），\(\psi_x\) 近似为高斯分布，其标准差记为 \(\psi_0\)。少数近距离的“硬”碰撞会在分布外侧形成比高斯更长的尾巴。
原子数密度 \(n_{\rm A}=\rho N_{\rm A}/A\)。走过 \(\dd x\) 时碰撞参数在 \(b\)–\(b+\dd b\) 的核数为 \(n_{\rm A}2\pi b\,\dd b\,\dd x\)，每次贡献空间偏转角平方 \(\psi(b)^2\)。求和： \begin{align*} \frac{\dd\avg{\psi_{\rm sp}^2}}{\dd x}&=n_{\rm A}\int\psi(b)^2\,2\pi b\,\dd b\\ &=8\pi n_{\rm A}z^2Z^2r_{\rm e}^2\Big(\frac{m_{\rm e}c^2}{\beta pc}\Big)^2\\ &\quad\times\ln\frac{b_{\max}}{b_{\min}} . \end{align*} 积分与能量损失的积分（第 4.2 节）完全同形，也是对数发散：\(b_{\max}\) 取原子电子云对核电荷的屏蔽半径，\(b_{\min}\) 取原子核半径。空间角是两个垂直投影的合成，\(\psi_{\rm sp}^2=\psi_x^2+\psi_y^2\)，故投影角的均方为其一半。原子中 \(Z\) 个电子也会使 \(\mu\) 子偏转，粗略计入的办法是把 \(Z^2\) 换成 \(Z(Z+1)\)。于是 \[\tag{6.2} \avg{\psi_x^2}\propto\frac{Z(Z+1)\rho}{A}\cdot\frac{x}{(\beta pc)^2}.\]
6.3 辐射长度与 Highland 公式
高能电子在物质中主要通过轫致辐射损失能量，平均每走过一个辐射长度 \(X_0\)，能量降为原来的 \(1/\ee\)。PDG 给出[3] \begin{align*} \frac{1}{X_0}\approx{}&4\alpha r_{\rm e}^2\frac{N_{\rm A}}{A}Z(Z+1) \\ &\times\ln\big(184\,Z^{-1/3}\big) \tag{6.3} \end{align*} （\(\alpha=1/137\) 为精细结构常数。此式略去了库仑修正等项，对铅、铀偏小约 10%，对铁小于 1%；本讲的数值一律用 PDG 表值）。表值：水 36.08 g/cm2（36.1 cm），标准岩石 26.54 g/cm2，铁 13.84 g/cm2（1.76 cm），铅 6.37 g/cm2（0.56 cm），铀 6.00 g/cm2（0.32 cm）[3]。
为什么一个描述“电子辐射”的长度会出现在 \(\mu\) 子的散射里？比较式 (6.2) 与式 (6.3)：两者都正比于 \(Z(Z+1)\rho/A\)——轫致辐射本质上也是电子在原子核库仑场中的“偏转”。所以 \(x/X_0\) 恰好就是散射公式里与材料有关的那个组合。附录 B 计算两个对数之比（约为 2），得到 \begin{align*} \avg{\psi_{\rm sp}^2}&\approx\Big(\frac{E_{\rm s}}{\beta pc}\Big)^2\frac{x}{X_0}, \\ E_{\rm s}&=m_{\rm e}c^2\sqrt{\frac{4\pi}{\alpha}}=21.2\ \mathrm{MeV}, \tag{6.4} \end{align*} 即 Rossi 与 Greisen 1941 年的结果[20]。
实用中，投影角分布中央 98% 用高斯近似，其标准差为[3], [21], [22] \begin{align*} \psi_0={}&\frac{13.6\ \mathrm{MeV}}{\beta pc}\,z\sqrt{\frac{x}{X_0}} \\ &\times\left[1+0.038\ln\frac{xz^2}{X_0\beta^2}\right] . \tag{6.5} \end{align*} 适用于 \(10^{-3}<x/X_0<100\)，精度约 11%。与式 (6.4) 的投影版本 \(E_{\rm s}/\sqrt2={15.0\,}\mathrm{MeV}\) 相比，系数 13.6 与对数修正来自对精确的 Molière 理论分布的拟合：中央部分比“全部散射角的均方根”窄，因为后者被少数大角单次散射拉大。多层材料要先把 \(x/X_0\) 相加，再代入一次式 (6.5)（不能把各层的 \(\psi_0\) 平方相加）[3]。
\(p={3\,}\mathrm{GeV/c}\) 的 \(\mu\) 子（\(\beta\approx1\)）穿过 10 cm 铁：\(x/X_0=10/1.757=5.69\)， \begin{align*} \psi_0&=\frac{13.6}{3000}\sqrt{5.69}\,(1+0.038\ln5.69)\\&=11.5\ \mathrm{mrad}. \end{align*} 换成 10 cm 铅：\(x/X_0=17.8\)，\(\psi_0={21.2\,}\mathrm{mrad}\)；铀：\(x/X_0=31.6\)，\(\psi_0={28.8\,}\mathrm{mrad}\)。动量加倍，\(\psi_0\) 减半；厚度加倍，\(\psi_0\) 增大 \(\sqrt2\) 倍（略多，来自对数项）。
7 宇宙线 \(\mu\) 子
7.1 立体角与强度
从一点 O 看一块小面元 \(\dd A\)，面元到 O 的距离为 \(r\)，面元法线与视线夹角为 \(\chi\)。面元在垂直于视线方向上的投影为 \(\dd A\cos\chi\)，它对 O 所张的立体角定义为 \[\tag{7.1} \dd\Omega=\frac{\dd A\cos\chi}{r^2},\] 单位为球面度（sr），整个球面对球心张 \(4\pi\) sr。取以 O 为心的单位球面上的面元（\(\chi=0\)、\(r=1\)），用天顶角 \(\theta\) 与方位角 \(\varphi\) 表示，面元为 \((\dd\theta)(\sin\theta\,\dd\varphi)\)，所以 \(\dd\Omega=\sin\theta\,\dd\theta\,\dd\varphi\)。
（积分）强度 \(I(\theta,\varphi)\)：单位时间内、穿过垂直于方向 \((\theta,\varphi)\) 的单位面积、方向落在单位立体角内的粒子数，单位 \(\mathrm{m^{-2}s^{-1}sr^{-1}}\)。微分强度 \(I(p,\theta)\) 再“按动量分档”，多一个 \((\mathrm{GeV}/c)^{-1}\)；对动量从 \(p\) 积到无穷，得积分强度 \(I(>p,\theta)\)。
穿过平面的通量 \(J\)：单位时间穿过某个平面单位面积的粒子数。方向与平面法线成 \(\theta\) 角的粒子“看到”的平面面积缩小为 \(\cos\theta\) 倍（投影），所以 \[\tag{7.2} J=\int I(\theta,\varphi)\cos\theta\,\dd\Omega .\]
7.2 从初级宇宙线到 \(\mu\) 子
从宇宙空间射到大气顶部的高能粒子叫初级宇宙线，主要是质子，其次是氦核和少量更重的原子核。在几 GeV 到约 \(10^{6}\) GeV 之间，其强度随能量近似按 \(E^{-2.7}\) 下降[16]。它们在大气中与氮、氧原子核碰撞，产生大量次级粒子，统称空气簇射（图 5）。
质子在空气中的核相互作用长度 \(\lambda_I={90.1\,}\mathrm{g/cm^2}\)[3]：平均每穿过 90.1 g/cm2 的空气发生一次核碰撞。首次碰撞所在的深度服从均值为 \(\lambda_I\) 的指数分布，所以首次碰撞的平均深度为 \(X=\lambda_I\)。 竖直大气深度 \(X(h)\) 是高度 \(h\) 以上单位面积空气柱的质量，由流体静力学平衡 \(\dd p/\dd h=-\rho g\) 得 \(X(h)=p(h)/g\)；海平面 \(X={1033\,}\mathrm{g/cm^2}\)。若某一层大气等温（温度 \(T_{\rm a}\)），由理想气体 \(\rho=pM_{\rm a}/(R_{\rm g}T_{\rm a})\)（\(M_{\rm a}\) 为空气摩尔质量，\(R_{\rm g}\) 为气体常数）可得 \(\dd p/p=-\dd h/H\)，即在这一层内 \begin{align*} X(h)&=X(h_1)\,\ee^{-(h-h_1)/H}, \\ H&=\frac{R_{\rm g}T_{\rm a}}{M_{\rm a}g}, \\ \rho(h)&=\frac{X(h)}{H}. \tag{7.3} \end{align*} 这里 \(h_1\) 是该层的下边界。真实大气在 11 km 以下（对流层）温度随高度近似线性下降，在 11–20 km（平流层下部）近似等温，\(T_{\rm a}\approx{217\,}\mathrm{K}\)，\(H={6.34\,}\mathrm{km}\)（1976 年美国标准大气）。对流层顶 \(h_1={11\,}\mathrm{km}\) 处的深度由标准大气算得 \(X(h_1)={231\,}\mathrm{g/cm^2}\)。令 \(X=\lambda_I\)：\(h={11\,}\mathrm{km}+H\ln(231/90.1)={17.0\,}\mathrm{km}\)；对标准大气直接数值求解也得 17.0 km。（若把 217 K 的等温大气一直用到海平面，会得到 15.5 km；但对流层并不等温，这个数与下面的 15 km 接近只是巧合。）
PDG 说“\(\mu\) 子大多产生在约 15 km 高处”[16]，依据的是另一个量：\(\mu\) 子的母粒子 \(\pi\) 介子（下一小节介绍）由首次及后续几代碰撞不断产生，同时母核子流按 \(\ee^{-X/\Lambda}\) 衰减（\(\Lambda\approx{121\,}\mathrm{g/cm^2}\) 为核子流的衰减长度，比 \(\lambda_I\) 长，因为碰撞后的领头核子仍带走相当一部分能量），低能（\(E_\pi\ll E_\pi^*\)，几乎全部衰变）的 \(\pi\) 的强度在 \(X=\Lambda\) 处达到极大[16]。按标准大气，\(X={121\,}\mathrm{g/cm^2}\) 对应 15.1 km。
7.3 \(\pi\) 介子：衰变还是再碰撞
\(\pi\) 介子是最轻的一类强相互作用粒子（介子），有 \(\pi^+\)、\(\pi^-\)、\(\pi^0\) 三种；中微子是几乎不与物质发生作用的中性轻子；K 介子是比 \(\pi\) 重的另一类介子。碰撞产生的带电 \(\pi\) 介子（\(m_\pi c^2={139.6\,}\mathrm{MeV}\)，静止寿命 26.0 ns，\(c\tau_\pi={7.80\,}\mathrm{m}\)[3]）有两种命运：衰变成 \(\mu\) 子和中微子（\(\pi^+\to\mu^+\nu_\mu\)，\(\pi^-\to\mu^-\bar\nu_\mu\)），或者在衰变前再次撞上原子核。
考虑竖直向下、能量 \(E\) 的 \(\pi\) 介子位于平流层下部、大气深度 \(X\) 处（该层近似等温）。它走过面密度 \(\dd X\)，对应路程 \(\dd l=\dd X/\rho=H\,\dd X/X\)（式 (7.3)）。
衰变概率：由式 (3.7)，每走 \(\dd l\) 衰变的概率为 \(\dd l/(\gamma c\tau_\pi)\)（\(\beta\approx1\)），即 \(\dfrac{m_\pi c^2}{E}\dfrac{H}{c\tau_\pi}\dfrac{\dd X}{X}\)；
相互作用概率：\(\dd X/\lambda\)（\(\lambda\) 为 \(\pi\) 的相互作用长度，与 \(\lambda_I\) 同量级）。
两者之比为 \(\dfrac{E_\pi^*}{E}\cdot\dfrac{\lambda}{X}\)，其中 \begin{align*} E_\pi^*&\equiv\frac{m_\pi c^2\,H}{c\tau_\pi} \\ &=\frac{0.1396\ \mathrm{GeV}\times6340\ \mathrm{m}}{7.80\ \mathrm{m}} \\ &=113\ \mathrm{GeV}. \tag{7.4} \end{align*} 在 \(\pi\) 产生的深度 \(X\sim\lambda\) 处，\(E\ll E_\pi^*\) 的 \(\pi\) 几乎都衰变，\(E\gg E_\pi^*\) 的 \(\pi\) 多半先碰撞。对带电 K 介子（\(m_Kc^2={493.7\,}\mathrm{MeV}\)，\(c\tau_K={3.71\,}\mathrm{m}\)）同样得到 \(E_K^*={844\,}\mathrm{GeV}\)。
\(E_\mu\ll E_\pi^*\)：\(\pi\) 衰变概率接近 1，\(\mu\) 子谱基本“继承”初级谱的形状，\(\dd I/\dd E\propto E^{-2.7}\)（在 10–100 GeV 范围成立；更低能量处还受 \(\mu\) 子自身衰变和能量损失的影响而变平）；
\(E_\mu\gg E_\pi^*\)：衰变概率 \(\propto E_\pi^*/E\)，\(\mu\) 子谱比初级谱陡一个幂次，\(\propto E^{-3.7}\)。
PDG 给出的 Gaisser 公式正体现了这一点（适用于 \(E_\mu>100/\cos\theta\) GeV、\(\theta<70^\circ\)）[16]： \begin{align*} \frac{\dd N_\mu}{\dd E_\mu\,\dd\Omega}\approx{}&\frac{0.14\,E_\mu^{-2.7}}{\mathrm{cm^2\,s\,sr\,GeV}} \\ &\times\bigg[\frac{1}{1+\frac{1.1E_\mu\cos\theta}{115\ \mathrm{GeV}}} \\ &\quad+\frac{0.054}{1+\frac{1.1E_\mu\cos\theta}{850\ \mathrm{GeV}}}\bigg], \tag{7.5} \end{align*} 式中 \(E_\mu\) 以 GeV 为单位。两个分母中的 115 GeV 与 850 GeV 正是式 (7.4) 估计的 \(E_\pi^*\) 与 \(E_K^*\)。因子 1.1 反映 \(\mu\) 子只继承母介子的一部分能量（由衰变运动学，\(E_\mu\) 在 \(0.57E_\pi\) 到 \(E_\pi\) 之间），按陡峭的能谱加权后近似为 \(E_\pi\approx1.1E_\mu\)；这是 Gaisser 解析近似中的有效系数，本讲不推导。\(\cos\theta\) 的出现是因为斜着穿过大气时，同样的深度处空气更稀薄（\(\rho=X\cos\theta/H\)），\(\pi\) 更容易衰变，分界能量变为 \(E_\pi^*/\cos\theta\)。
7.4 通量与角分布
强度与通量的定义见第 7.1 节。
设 \(I=I_0\cos^n\theta\)（\(0\le\theta\le\pi/2\)，与 \(\varphi\) 无关）。代入式 (7.2)，令 \(u=\cos\theta\)： \begin{align*} J&=2\pi I_0\int_0^{\pi/2}\cos^{n+1}\theta\sin\theta\,\dd\theta \\ &=2\pi I_0\int_0^1u^{n+1}\dd u=\frac{2\pi I_0}{n+2}. \tag{7.6} \end{align*} 海平面实测：动量大于 1 GeV/c 的竖直强度 \(I_0\approx{70\,}\mathrm{m^{-2}\cdot s^{-1}\cdot sr^{-1}}\)，角分布 \(n\approx2\)[16]。于是 \begin{align*} J&=\pi I_0/2\\&=110\ \mathrm{m^{-2}.s^{-1}}\\&=0.66\ \mathrm{cm^{-2}.min^{-1}}; \end{align*} 计入 1 GeV/c 以下的 \(\mu\) 子，实验者常记的数是“水平探测器约 1 cm−2 min−1”[16]。注意：若忘了式 (7.2) 中的 \(\cos\theta\)，结果会变成 \(2\pi I_0/(n+1)\)，偏大 33%。
天顶角 \(\theta\) 方向的 \(\mu\) 子要穿过 \(\sec\theta\) 倍的大气：一方面路程更长、飞行时间更长，衰变得更多；另一方面能量损失更多（\(\sec\theta\times{2\,}\mathrm{GeV}\)），低能 \(\mu\) 子到不了地面。两者都使斜向强度降低。\(\cos^2\theta\) 是 \(E_\mu\sim{3\,}\mathrm{GeV}\) 附近的典型行为；能量更低时分布更陡，能量更高（\(E_\mu\gg E_\pi^*\)）时反而接近 \(\sec\theta\)（上一小节的 \(\cos\theta\) 效应占优）[16]。所以“\(n\)”与探测器的能量阈值有关，这正是实验题（第 11 节）要测的量。
7.5 能谱
本讲需要在 1 GeV 到几 TeV 的整个范围内使用海平面 \(\mu\) 子谱。我们采用 Reyna 对多组实验数据的拟合[23]：先给出竖直方向的微分强度 \[\tag{7.7} I_{\rm v}(p)=c_1\,p^{-(c_2+c_3\lg p+c_4\lg^2p+c_5\lg^3p)},\] \(\lg\) 表示以 10 为底的对数，\(p\) 以 GeV/\(c\) 计，\(I_{\rm v}\) 以 \(\mathrm{cm^{-2}s^{-1}sr^{-1}(GeV/c)^{-1}}\) 计，\(c_1=0.00253\)，\(c_2=0.2455\)，\(c_3=1.288\)，\(c_4=-0.2555\)，\(c_5=0.0209\)；再用经验标度关系推到任意天顶角： \[\tag{7.8} I(p,\theta)=\cos^3\theta\;I_{\rm v}(p\cos\theta).\] 适用范围 \(p>{1\,}\mathrm{GeV/c}\)、\(p<2000/\cos\theta\) GeV/\(c\)。这是对数据的经验拟合，没有简单的理论推导，我们把它当作“测得的谱”使用。
对式 (7.7) 从 1 GeV/c 积分到无穷，得竖直积分强度 70.2 m−2 s−1 sr−1，与 PDG 的 \(\approx70\) 一致（脚本 fig_spectrum.py 用 assert 核对）。与式 (7.5)（竖直）相比：100–300 GeV/c 内两者相差不到 13%，但 Gaisser 公式在更高能量处偏低，500 GeV/c 处约低 20%，1 TeV/c 处约低 30%（图 6a）。不同参数化在高能段的这种差别，就是火山成像（\(E_{\min}\) 为几百 GeV）中谱形不确定度的量级（第 8.5 节）。
\[\tag{7.9} \eta(p)\equiv-\frac{\dd\ln I(>p)}{\dd\ln p}.\] 若 \(I(>p)\propto p^{-\eta}\) 是严格幂律，\(\eta\) 就是那个幂指数；对真实谱，\(\eta\) 随 \(p\) 缓慢变化。由图 6b：\(p={10\,}\mathrm{GeV/c}\) 时 \(\eta\approx1.4\)，45 GeV/c 时 \(\eta\approx1.9\)，700 GeV/c 时 \(\eta\approx2.3\)。\(\eta\) 将决定透射成像的灵敏度。
7.6 大气中的能量损失与存活
整层大气的面密度 1033 g/cm2，按最小电离约 2 MeV cm2/g 计，竖直穿过大气的 \(\mu\) 子损失约 2 GeV，这与 PDG 的说法一致[16]。能量变化使式 (3.7) 中的衰变长度沿途变化，固有时要逐段积分： \[\tag{7.10} \tau=\int\frac{m\,\dd x}{p(x)},\qquad P=\ee^{-\tau/\tau_\mu}.\] \(p(x)\) 由能损方程 \(\dd E/\dd x=-\rho(x)\,(\dd E/\dd X)\) 与 \(E^2=(pc)^2+(mc^2)^2\) 给出。
PDG 写道：“能量 2.4 GeV 的 \(\mu\) 子衰变长度为 15 km，计入能量损失后缩短为 8.7 km”[16]。由式 (3.7)，\(\lambda_{\rm d}=(2.4/0.1057)\times{658.6\,}\mathrm{m}\times\beta={14.9\,}\mathrm{km}\)。本讲的脚本（atmos.py）从 15 km 高处以 2.4 GeV 出发，按标准大气和 PDG 空气能损表积分式 (7.10)，\(\mu\) 子到达地面时剩 0.51 GeV，\(\tau=1.67\tau_\mu\)；定义“有效衰变长度”\({15\,}\mathrm{km}/(\tau/\tau_\mu)={9.0\,}\mathrm{km}\)，与 PDG 的 8.7 km 相差 3%。
7.7 Frisch–Smith 实验
Frisch 与 Smith 的实验设计很巧妙[6]。探测器是一大块塑料闪烁体，上方放置吸收体。只有那些恰好停在闪烁体里的 \(\mu\) 子才被记录：\(\mu\) 子进入闪烁体时发出一次闪光，停下来后衰变出的电子发出第二次闪光，两次闪光的时间差就是这个 \(\mu\) 子静止时的寿命样本。这一设计一举两得：
能停在闪烁体里，意味着 \(\mu\) 子在进入吸收体前的能量落在一个窄范围内，即选出了一个确定的速度窗口（山顶处 \(0.9950c\)–\(0.9954c\)）；
同一台仪器测出的衰变时间分布给出静止 \(\mu\) 子的寿命，就是“\(\mu\) 子钟”的定标。
在山顶每小时约记录 563 个，在高度低约 1.9 km 的海平面每小时约 412 个（据文献转述[5]）。若没有时间膨胀，\(\mu\) 子从山顶飞到海平面用时 \(t\approx{1.9\,}\mathrm{km}/(0.995c)={6.4\,}\mathrm{\mu s}\)，约 2.9 个寿命，海平面只应剩 \(563\times\ee^{-2.9}\approx30\) 个/小时（作者用实测衰变曲线算得约 27 个）。实测 412 个，意味着 \(\mu\) 子在这段路上只经历了约 0.7 μs 的固有时，时间膨胀因子约 9。
还有一个细节：\(\mu\) 子在 1.9 km 空气中损失几百 MeV，速度从 \(0.995c\) 降到约 \(0.988c\)，\(\gamma\) 从约 10 降到约 6.3，所以“膨胀因子”沿途变化。例题 1 将定量处理，并说明为什么海平面的吸收体必须比山顶的薄。（据转述，Frisch 与 Smith 给出的海平面速度为 \(0.9881c\)–\(0.9897c\)、\(\gamma\approx6.8\)[5]，比本讲按 PDG 空气能损与标准大气算出的 6.3 大约 8%。差别可能来自他们所用的能损数据或空气柱面密度，原文正文未能核实。若取 \(\gamma_2=6.8\)，固有时约减小 4%，有效时间膨胀因子正好是他们给出的 8.4。）
交互演示加载中……（需启用 JavaScript）
8 透射成像
8.1 测量量：沿视线的不透明度
能测出带电粒子穿过位置的探测平面叫位置灵敏探测器；两层以上这样的平面上下排列，由各层的击中点连线得到粒子方向，整套装置叫径迹探测器（或“望远镜”）。只有在很短时间窗口内各层同时有信号才记为一个事例，叫符合，用以排除噪声。径迹方向常用它在两个竖直平面（东西、南北）上投影的倾角正切表示：若方向单位矢量为 \(\uvec n=(n_x,n_y,n_z)\)（\(x\) 向东、\(y\) 向北、\(z\) 向上），则 \(\tan\theta_x=n_x/n_z\)，\(\tan\theta_y=n_y/n_z\)（附录 D）。本讲提到的几种：
把能测出 \(\mu\) 子入射方向的探测器（径迹探测器）放在待测物体下方或内部。对方向 \((\theta,\varphi)\) 附近的一个小立体角 \(\Delta\Omega\)，曝光时间 \(T\) 内的期望计数为 \begin{align*} \mu_N&=\Phi\big(X(\theta,\varphi),\theta\big)\,S_{\rm eff}\,\Delta\Omega\,T, \\ \Phi(X,\theta)&\equiv I\big(>E_{\min}(X),\theta\big), \tag{8.1} \end{align*} 其中 \(X(\theta,\varphi)=\int\rho\,\dd l\) 沿这条视线从探测器积分到物体表面，\(S_{\rm eff}\) 是探测器对该方向的有效面积（平板探测器为 \(S\cos\theta_{\rm d}\)，\(\theta_{\rm d}\) 是视线与探测器法线的夹角），\(E_{\min}\) 由式 (4.8) 给出。测出 \(\mu_N\) 的估计值 \(N\)，反解出 \(X\)，再除以由几何测量得到的路径长度 \(L\)，就得到这条视线上的平均密度 \(\bar\rho=X/L\)。一张“\(\mu\) 子照片”的每个像素就是一条视线上的 \(X\) 或 \(\bar\rho\)（图 7）。
探测器对某个方向范围的接收度定义为 \(G=\int S_{\rm eff}(\theta,\varphi)\,\dd\Omega\)，单位 \(\mathrm{m^2\,sr}\)，即“有效面积 × 立体角”。若强度在该范围内近似不变，计数 \(=I\times G\times T\)。式 (8.1) 中的 \(S_{\rm eff}\Delta\Omega\) 就是小立体角 \(\Delta\Omega\) 的接收度。
8.2 透射强度对不透明度的灵敏度
由链式法则， \[\frac{\dd\ln\Phi}{\dd\ln X}=\frac{\dd\ln I(>E)}{\dd\ln E}\bigg|_{E=E_{\min}}\cdot\frac{\dd\ln E_{\min}}{\dd\ln X}.\] 第一个因子是 \(-\eta(E_{\min})\)（式 (7.9)，在 \(E\gg mc^2\) 时不区分 \(E\) 与 \(pc\)）。第二个因子由式 (4.8) 求导： \begin{align*} \frac{\dd\ln E_{\min}}{\dd\ln X}&=\frac{X}{E_{\min}}\frac{\dd E_{\min}}{\dd X}=\frac{X\cdot a\,\ee^{bX}}{(a/b)(\ee^{bX}-1)}\\ &=\frac{bX\,\ee^{bX}}{\ee^{bX}-1}\equiv g(bX). \end{align*} 于是 \begin{align*} \eta'&\equiv-\frac{\dd\ln\Phi}{\dd\ln X} \\ &=\eta(E_{\min})\,g(bX), \\ g(u)&=\frac{u\,\ee^{u}}{\ee^{u}-1} \\ &=1+\frac u2+\order{u^2}. \tag{8.2} \end{align*}
不透明度相对变化 \(\delta X/X\) 引起透射强度相对变化 \(\delta\Phi/\Phi\approx-\eta'\,\delta X/X\)。金字塔（\(X\approx{2\times 10^{4}\,}\mathrm{g/cm^2}\)，\(\eta\approx1.9\)，\(bX=0.08\)）：\(\eta'\approx2.0\)；火山（近水平视线 \(\theta\approx80^\circ\)，\(X\approx{1.25\times 10^{5}\,}\mathrm{g/cm^2}\)，\(bX=0.5\)，该天顶角下 \(\eta\approx2.0\)，\(g=1.27\)）：\(\eta'\approx2.5\)。（若误用竖直谱的斜率，\(E_{\min}\approx{360\,}\mathrm{GeV}\) 处 \(\eta=2.27\)，会得到 \(\eta'\approx2.9\)。）反过来，用 \(N\) 个计数测定的平均密度的统计精度为 \[\tag{8.3} \frac{\delta\bar\rho}{\bar\rho}=\frac{\delta X}{X}=\frac{1}{\eta'}\frac{\delta N}{N}=\frac{1}{\eta'\sqrt N}.\] 因为 \(\eta'\approx2\)，测 1% 的密度只需要约 \(2500\) 个计数。
8.3 空腔造成的超额
设视线上岩石密度均匀为 \(\rho\)，总长 \(L\)，其中一段长 \(\ell\) 是空腔（充满空气，密度可忽略）。不透明度从 \(X=\rho L\) 变为 \(X'=\rho(L-\ell)=X(1-\ell/L)\)。相对计数超额定义为 \[\varepsilon\equiv\frac{\Phi(X')-\Phi(X)}{\Phi(X)} .\] 若 \(\ell\ll L\)，对 \(\ln\Phi\) 在 \(\ln X\) 处作一阶泰勒展开： \begin{align*} &\ln\Phi(X')-\ln\Phi(X)\\&\approx-\eta'\,[\ln X'-\ln X]\\&=-\eta'\ln(1-\ell/L)\approx\eta'\ell/L. \end{align*} 于是 \[\tag{8.4} \varepsilon\approx\eta'\,\frac{\ell}{L}\qquad(\ell\ll L).\] 保留的是 \(\ell/L\) 的一阶项；略去的二阶项来自 \(\ln(1-\ell/L)\) 的展开、\(\ee^{(\cdot)}-1\) 的展开和 \(\eta'\) 随 \(X\) 的变化；前两项使超额增大，第三项使之略减小（\(\eta'\) 随 \(X\) 增大），合计使真实超额略大于一阶估计。
取王后墓室中的乳胶片看向大空洞的方向：天顶角约 \(29^\circ\)，岩石长度 \(L\approx{90\,}\mathrm{m}\)（文献给出该视场内岩石厚度 65–115 m[9]），石灰岩 \(\rho={2.2\,}\mathrm{g/cm^3}\)。\(X={1.98\times 10^{4}\,}\mathrm{g/cm^2}\)，\(E_{\min}={45\,}\mathrm{GeV}\)，在 \(\theta=29^\circ\) 的谱上 \(\eta=1.88\)，\(\eta'=1.96\)。
\(\ell={3\,}\mathrm{m}\)：式 (8.4) 给出 6.5%，用完整谱精确计算为 6.8%；
\(\ell={5\,}\mathrm{m}\)：10.9%，精确 11.8%。
KEK 径迹仪在不同切片中测得的超额为 5.2%–11%[9]，按式 (8.4) 相当于视线穿过约 2.5–5 m 的空腔，与“截面与大走廊相当”的描述相容。
交互演示加载中……（需启用 JavaScript）
8.4 统计：需要多少 \(\mu\) 子
记“无空腔”时窗口内的期望计数为 \(\mu_0\)，有空腔时为 \((1+\varepsilon)\mu_0\)。按式 (5.2)，期望的显著性 \[\tag{8.5} n_\sigma=\frac{(1+\varepsilon)\mu_0-\mu_0}{\sqrt{\mu_0}}=\varepsilon\sqrt{\mu_0}.\] 要达到 \(n_\sigma\)，需要 \(\mu_0=(n_\sigma/\varepsilon)^2\) 个本底计数；结合式 (8.1)，曝光时间 \[\tag{8.6} T=\frac{(n_\sigma/\varepsilon)^2}{\Phi\,S_{\rm eff}\,\Delta\Omega}.\] 超额减半，所需时间变为 4 倍；这是计数实验的“平方律”。地球物理中评估 \(\mu\) 子成像可行性时用的正是这一类估计[24]。
名古屋大学的乳胶片有效面积 0.45 m2，在位置 NE1 曝光 98 天[9]。取分析窗口 \(\Delta(\tan\theta_x)=0.05\)、\(0.4\le\tan\theta_y<0.7\)（\(\tan\theta_y\) 的范围是文献中用于展示超额的切片，\(\tan\theta_x\) 的窗宽为本讲的假设），由附录 D 的式 (D.1)， \begin{align*} \Delta\Omega&=0.05\times0.3\times(1+0.55^2)^{-3/2}\\&=0.0101\ \mathrm{sr}. \end{align*} 乳胶片平放，\(S_{\rm eff}=0.45\cos29^\circ={0.39\,}\mathrm{m^2}\)。由上一小节，\(\Phi=I(>{45\,}\mathrm{GeV},29^\circ)={0.68\,}\mathrm{m^{-2}\cdot s^{-1}\cdot sr^{-1}}\)。于是每天 \(\mu_0=0.68\times0.39\times0.0101\times86400=234\) 个，98 天 \(2.3\times10^4\) 个。取 \(\varepsilon=8\%\)：\(n_\sigma=0.08\times\sqrt{2.3\times10^4}=12.1\)。文献报告 13.7\(\sigma\)[9]。考虑到我们对窗口、岩石厚度和超额都只取了代表值，这个一致程度说明式 (8.5) 抓住了问题的主要部分。达到 5\(\sigma\) 只需 \((5/0.08)^2/234\approx17\) 天。
探测器若放在金字塔外面（如 CEA 的望远镜在北面外侧），看向空洞的视线要穿过更多石头，\(\Phi\) 迅速下降：\(X\) 增大一倍，\(\Phi\) 约降为 \(2^{-2}=1/4\)（\(\eta'\approx2\)）。同时 \(\ell/L\) 变小，超额也降低。CEA 用两台望远镜各两个月得到合并 5.8\(\sigma\)[9]，正反映了这一点。
交互演示加载中……（需启用 JavaScript）
8.5 系统误差与本底
若“无空腔”期望 \(\mu_0\) 本身有相对系统误差 \(\sigma_{\rm sys}\)（岩石密度与几何的不确定、谱的形状误差、探测效率的不均匀等），它与泊松涨落独立，方差相加：\(\sigma^2=\mu_0+\sigma_{\rm sys}^2\mu_0^2\)。于是 \begin{align*} n_\sigma&=\frac{\varepsilon\,\mu_0}{\sqrt{\mu_0+\sigma_{\rm sys}^2\mu_0^2}} \\ &=\frac{\varepsilon}{\sqrt{1/\mu_0+\sigma_{\rm sys}^2}} \\ &\xrightarrow{\ \mu_0\to\infty\ }\ \frac{\varepsilon}{\sigma_{\rm sys}}. \tag{8.7} \end{align*} 无论曝光多久，显著性都不会超过 \(\varepsilon/\sigma_{\rm sys}\)。\(\varepsilon=8\%\)、\(\sigma_{\rm sys}=2\%\) 时上限只有 4\(\sigma\)（图 8a）。
实际分析靠两个办法压低 \(\sigma_{\rm sys}\)：(i) 用物体内没有目标的区域对模拟做归一化，让谱的整体幅度、探测效率等大尺度误差相消（谱的整体幅度本身就有 10%–15% 的不确定[16]）。但谱的形状误差不能这样消去：火山成像的 \(E_{\min}\) 在几百 GeV，那里不同参数化相差 20%–30%（第 7.5 节），由式 (8.3) 相当于密度的系统误差 8%–12%，必须用同一探测器在开阔天空方向实测谱来标定；(ii) 寻找的是局部的鼓包，与相邻方向比较，只有在小角度尺度上起伏的误差才进入 \(\sigma_{\rm sys}\)。ScanPyramids 估计石块尺度（约 1 m）的表面起伏造成的相对误差小于 2%[9]。
假信号本底：低能电子、质子、被散射偏转的 \(\mu\) 子，甚至从背后入射被误判方向的粒子，会混入“穿过岩石”的计数。它们与岩石厚度无关，使厚岩石方向的计数偏高，因而把密度测低（例题 4）。火山观测常在探测平面之间夹铅板以吸收低能粒子[11]。
岩石中的多次散射：若按入射能量 45 GeV 计算，穿过 90 m 石灰岩（约 750 个辐射长度，已超出式 (6.5) 的适用范围，只作量级估计）投影角的高斯宽度 \(\psi_0\approx{10\,}\mathrm{mrad}\)（相应地 \(\abs{\psi_x}\) 的中位数只有 \(0.67\psi_0\approx{7\,}\mathrm{mrad}\)）。但能透射过来的 \(\mu\) 子在路径末端——恰好是“空腔到探测器”这一段——能量最低，而散射正比于 \(\int\dd X/(pc)^2\)，主要发生在末段。按透射谱抽样、计入沿途能损（脚本
rock_ms.py）：每个 \(\mu\) 子的 \(\psi_0\) 的中位数约 14 mrad，比按入射能量算的 10 mrad 大；\(\abs{\psi_x}\) 的中位数约 8 mrad（45 m 外约 0.4 m）。约 1/5 的 \(\mu\) 子出射能量低于 5 GeV，散射大得多，10% 的 \(\mu\) 子偏角超过约 40 mrad（45 m 外约 1.8 m）。所以空腔的像有一个清晰的核心和一圈 1–2 m 的晕；ScanPyramids 的模拟必须计入散射，原因正在于此。
9 散射成像
9.1 原理
透射成像只对“挡住了多少 \(\mu\) 子”敏感，适合几十米到几千米的大物体。对集装箱这样只有几米大、但可能藏着几厘米重金属的物体，大部分 \(\mu\) 子都能穿过，透射几乎不变；有用的信息在偏折里。在物体上方和下方各放一组径迹探测器，测出每个 \(\mu\) 子的入射方向与出射方向，两者之差就是散射角 \((\psi_x,\psi_y)\)；把入射与出射径迹外推，求两条直线的最近点，可以估计散射主要发生在哪里（“最近点法”）。对很多 \(\mu\) 子统计，就得到“散射强度”的三维分布。
由式 (6.5)，散射角的宽度由 \(x/X_0\) 决定，而 \(X_0\) 随原子序数急剧下降（式 (6.3) 中的 \(Z(Z+1)/A\)）。表 2 列出几种材料的数值。
| 材料 | \(X_0\)/cm | 厚度/cm | \(x/X_0\) | \(\psi_0\)/mrad |
|---|---|---|---|---|
| 水 | 36.1 | 100 | 2.77 | 7.8 |
| 铁 | 1.76 | 10 | 5.69 | 11.5 |
| 铅 | 0.561 | 10 | 17.8 | 21.2 |
| 铀 | 0.317 | 10 | 31.6 | 28.8 |
9.2 动量未知的代价
宇宙线 \(\mu\) 子的动量分布很宽，而简单的径迹探测器不测动量。每个 \(\mu\) 子的散射角来自宽度 \(\psi_0(p)\propto1/(\beta p)\) 的高斯分布，\(p\) 因 \(\mu\) 子而异，所以实测分布是许多不同宽度的高斯分布的混合。
对给定的 \(p\)，\(\avg{\psi_x^2\,|\,p}=\psi_0(p)^2\)。对 \(p\) 的分布求平均（忽略 Highland 公式中对数项随 \(p\) 的缓慢变化）： \begin{align*} \avg{\psi_x^2}&=\psi_0(p_{\rm ref})^2\,(\beta p)_{\rm ref}^2 \\ &\quad\times\avg{\frac{1}{(\beta p)^2}} \\ &\equiv\psi_0(p_{\rm ref})^2\,\frac{(\beta p)_{\rm ref}^2}{p_{\rm eff}^2}, \\ p_{\rm eff}&\equiv\avg{(\beta p)^{-2}}^{-1/2}. \tag{9.1} \end{align*} 高斯分布的四阶矩为 \(3\sigma^4\)，所以混合后 \(\avg{\psi_x^4}=3\avg{\psi_0^4}\)，峰度比 \(\kappa=3\avg{(\beta p)^{-4}}/\avg{(\beta p)^{-2}}^2\ge3\)（附录 C）。
\(p>{1\,}\mathrm{GeV/c}\) 用式 (7.7) 的竖直谱；\(p<{1\,}\mathrm{GeV/c}\) 取平坦谱（PDG：“1 GeV 以下谱几乎平坦”[16]）。这是一个模型假设。
只有穿透整个物体（以及下方探测器中的吸收层）的 \(\mu\) 子才会被下方径迹探测器记录，所以被记录的 \(\mu\) 子谱在某个截断动量 \(p_{\rm c}\) 以下为零。\(p_{\rm c}\) 由沿途的总能损决定：物体越厚、越重，\(p_{\rm c}\) 越高；同一台装置看“铀块”和“铅块”时截断不同（例题 3）。\(p_{\rm eff}\) 与 \(\kappa\) 对 \(p_{\rm c}\) 很敏感（表 3，按入射动量计算、不计物体内的减速）：
| \(p_{\rm c}\)/(GeV/\(c\)) | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 1.0 |
|---|---|---|---|---|---|---|---|
| \(p_{\rm eff}\)/(GeV/\(c\)) | 1.07 | 1.26 | 1.43 | 1.59 | 1.74 | 1.89 | 2.18 |
| \(\kappa\) | 15.9 | 12.4 | 10.5 | 9.28 | 8.46 | 7.87 | 7.09 |
以 \(p_{\rm c}={0.5\,}\mathrm{GeV/c}\) 为例：\(p_{\rm eff}={1.43\,}\mathrm{GeV/c}\)，\(\kappa=10.5\)。不测动量时，实测均方根散射角约为“3 GeV/c 单一动量”值的 \(3/1.43=2.1\) 倍，分布有很长的尾巴（图 10b）。\(\mu\) 子在物体内减速会使最慢的那部分散射得更厉害，实际的 \(\kappa\) 还要更大；单次大角散射（Highland 公式只描述中央 98%）也会进一步加长尾巴。
9.3 所需 \(\mu\) 子数
每个 \(\mu\) 子测得两个投影角 \(\psi_x,\psi_y\)。给定动量 \(p\) 时，二者是宽度同为 \(\sigma(p)=\psi_0(p)\) 的独立高斯变量；但 \(p\) 未知且对两者相同：动量小的 \(\mu\) 子两个投影都大，动量大的两个都小。所以 \(\psi_x^2\) 与 \(\psi_y^2\) 正相关，不能当作 \(2N\) 个独立样本。正确的做法是把每个 \(\mu\) 子当作一个样本，取 \(m_i=\tfrac12(\psi_{x,i}^2+\psi_{y,i}^2)\)，\(N\) 个 \(\mu\) 子的 \(m_i\) 相互独立。先对给定 \(p\) 求条件矩，再对 \(p\) 平均（附录 C）： \begin{align*}
\avg{m}&=\avg{\sigma^2},\\
\avg{m^2}&=\tfrac14\avg{3\sigma^4+2\sigma^4+3\sigma^4}=2\avg{\sigma^4},
\end{align*} \begin{align*}
\frac{\operatorname{Var}m}{\avg{m}^2}&=\frac{2\avg{\sigma^4}-\avg{\sigma^2}^2}{\avg{\sigma^2}^2} \\
&=\frac{2\kappa}{3}-1 . \tag{9.2}
\end{align*} 单一动量时 \(\kappa=3\)，结果为 1，与“\(2N\) 个独立高斯样本”的 \(\sqrt{2/(2N)}\) 一致；\(\kappa=10.5\) 时为 6.0，而错误地把两投影当独立样本会得到 \((\kappa-1)/2=4.75\)，低估 26%。蒙特卡罗（ex3_numbers.py，按例题 3 中铅假设的截断抽样，\(\kappa=10.0\)；每组 300 个 \(\mu\) 子、重复 \(2\times10^4\) 组）测得的相对误差与式 (9.2) 相符到 0.1%。
设两个假设（如“铀块”与“铅块”）预言的均方投影角为 \(Q_1>Q_2\)，比值 \(R=Q_1/Q_2\)。由式 (9.2)，\(N\) 个 \(\mu\) 子的均方估计 \(\bar m\) 的相对误差为 \(\sqrt{(2\kappa/3-1)/N}\)。要求两种预言之差 \(Q_1-Q_2=(R-1)Q_2\) 达到 \(n_\sigma\) 个标准差： \begin{align*} &(R-1)\ge n_\sigma\sqrt{\frac{2\kappa/3-1}{N}} \\ \Rightarrow\ &N\ge\Big(\frac{2\kappa}{3}-1\Big)\Big(\frac{n_\sigma}{R-1}\Big)^2 . \tag{9.3} \end{align*}
这里用了一个近似：误差按假设 2 的预言 \(Q_2\) 与峰度计算。例题 3 把能损截断、两投影相关、“走满厚度”的接收度都计入后，区分 10 cm 铀块与铅块（\(n_\sigma=5\)）需要约 700 个 \(\mu\) 子、约半小时；计入物体内减速的完整数值模型给出约 20–90 分钟（取决于下方是否有剔除慢 \(\mu\) 子的吸收层，也与块体在货物中的竖直位置有关）。实际系统常同时粗测 \(\mu\) 子动量（例如利用探测器自身各层中的散射），把慢 \(\mu\) 子单独处理，以缩短时间。
交互演示加载中……（需启用 JavaScript）
10 决赛风格例题
山顶探测器比海平面探测器高 \(h={1.91\,}\mathrm{km}\)。山顶选出的 \(\mu\) 子速度取 \(v_1=0.9952c\)。山顶计数率 \(N_1=563\) 个/小时，海平面计数率 \(N_2=412\) 个/小时。已知 \(\tau_\mu={2.197\,}\mathrm{\mu s}\)，\(m_\mu c^2={105.66\,}\mathrm{MeV}\)；两探测器之间的竖直空气柱面密度 \(\Delta X={213.5\,}\mathrm{g/cm^2}\)；在 \(\beta\gamma\approx6\)–10 时，\(\mu\) 子在空气中的能损取常数 \(a_{\rm air}={1.95\,}\mathrm{MeV\cdot cm^2/g}\)，在铁中取 \(a_{\rm Fe}={1.54\,}\mathrm{MeV\cdot cm^2/g}\)；铁的密度 7.87 g/cm3。
若没有时间膨胀，海平面每小时应有多少个？
若有时间膨胀但忽略能量损失，海平面每小时应有多少个？
计入能量损失。(a) 求 \(\mu\) 子到达海平面时的 \(\gamma_2\) 与速度。(b) 把空气密度视为均匀，证明 \(\mu\) 子飞行的固有时与实验室时间分别为 \begin{align*} \tau&=\frac{\operatorname{arcosh}\gamma_1-\operatorname{arcosh}\gamma_2}{kc},\\ t&=\frac{\sqrt{\gamma_1^2-1}-\sqrt{\gamma_2^2-1}}{kc},\\ k&\equiv\frac{\gamma_1-\gamma_2}{h}. \end{align*} （\(\operatorname{arcosh}x=\ln(x+\sqrt{x^2-1})\)。）(c) 求预期计数率与“有效时间膨胀因子”\(t/\tau\)。
由实测计数率反推时间膨胀因子。若山顶与海平面各计数 1 小时，统计误差多大？要使时间膨胀因子的统计误差降到 10%，两地各需计数多久？
海平面探测器上方的吸收体若与山顶相同，选出的是不是“同一批 \(\mu\) 子”？若要选出在山顶速度为 \(v_1\) 的那批 \(\mu\) 子，海平面的铁吸收体应比山顶薄多少？
在 \(\mu\) 子的瞬时静止参考系中，它能“活着”到达海平面的原因是迎面而来的大气层被洛伦兹收缩。写出该参考系中“收缩后的空气层逐段掠过 \(\mu\) 子所需时间”的积分，证明它等于第 (3)(b) 问的 \(\tau\)。
(1) 实验室飞行时间 \begin{align*} t_0&=h/v_1\\&=1910\ \mathrm{m}/(0.9952\times2.998\times10^{8}\ \mathrm{m/s})\\&=6.40\ \mu\mathrm{s}=2.91\tau_\mu. \end{align*} 由式 (3.6)（此时固有时等于实验室时），\(N_1\ee^{-2.91}=563\times0.0543=31\) 个/小时（作者按实测衰变曲线算得 27 个，差别来自衰变时间分布的细节）。
(2) \(\gamma_1=1/\sqrt{1-0.9952^2}=10.22\)，固有时 \(\tau=t_0/\gamma_1={0.627\,}\mathrm{\mu s}\)，预期 \(563\,\ee^{-0.627/2.197}=423\) 个/小时。
(3)(a) 能量 \(E=\gamma m_\mu c^2\)。穿过 \(\Delta X\) 损失 \(a_{\rm air}\Delta X=1.95\times213.5\ \mathrm{MeV}={416\,}\mathrm{MeV}\)，故 \(\gamma_2=\gamma_1-416/105.66=10.22-3.94=6.28\)，\(v_2=c\sqrt{1-1/\gamma_2^2}=0.9872c\)。（用 PDG 空气能损表与标准大气数值积分得 \(\gamma_2=6.27\)。）
(b) 能损率均匀、密度均匀时，\(\gamma\) 随路程线性减小：\(\gamma(x)=\gamma_1-kx\)。由式 (3.5)，\(\dd\tau=\dd x/(c\sqrt{\gamma^2-1})\)（用了 \(\beta\gamma=\sqrt{\gamma^2-1}\)）。换元 \(\dd x=-\dd\gamma/k\)： \begin{align*} \tau&=\frac1{kc}\int_{\gamma_2}^{\gamma_1}\frac{\dd\gamma}{\sqrt{\gamma^2-1}}\\ &=\frac{\operatorname{arcosh}\gamma_1-\operatorname{arcosh}\gamma_2}{kc}. \end{align*} 实验室时间 \(\dd t=\dd x/v=\dd x\,\gamma/(c\sqrt{\gamma^2-1})\)，同样换元：\(t=\frac1{kc}\int_{\gamma_2}^{\gamma_1}\frac{\gamma\,\dd\gamma}{\sqrt{\gamma^2-1}}=\frac{\sqrt{\gamma_1^2-1}-\sqrt{\gamma_2^2-1}}{kc}\)。
(c) \begin{align*} kc&=3.94\times2.998\times10^{8}\ \mathrm{m/s}/1910\ \mathrm{m}\\&=6.18\times10^{5}\ \mathrm{s^{-1}}, \end{align*} 即 \(1/(kc)={1.617\,}\mathrm{\mu s}\)。\(\operatorname{arcosh}x=\ln(x+\sqrt{x^2-1})\)：\(\operatorname{arcosh}10.22=3.015\)，\(\operatorname{arcosh}6.28=2.524\)。于是 \(\tau=0.491\times{1.617\,}\mathrm{\mu s}={0.794\,}\mathrm{\mu s}\)，\(t=(10.17-6.20)\times{1.617\,}\mathrm{\mu s}={6.42\,}\mathrm{\mu s}\)。预期计数 \(563\,\ee^{-0.794/2.197}=392\) 个/小时；有效时间膨胀因子 \(t/\tau=8.1\)，介于 \(\gamma_1=10.2\) 与 \(\gamma_2=6.3\) 之间（Frisch 与 Smith 给出 \(8.4\pm2\)[6]）。按标准大气（密度随高度变化）数值积分得 \(\tau={0.791\,}\mathrm{\mu s}\)、393 个/小时：密度不均匀只改变 \(\gamma\) 随 \(x\) 变化的快慢分布，对固有时影响很小。 检验：\(\gamma_2\to\gamma_1\)（\(k\to0\)）时， \begin{align*} \tau&\to(\gamma_1-\gamma_2)/(kc\sqrt{\gamma_1^2-1})\\&=h/(c\beta_1\gamma_1)\\&=t_0/\gamma_1, \end{align*} 回到第 (2) 问。
(4) 由 \(N_2=N_1\ee^{-\tau/\tau_\mu}\)， \begin{align*} \tau&=\tau_\mu\ln(N_1/N_2)\\&=2.197\times\ln(563/412)\ \mu\mathrm{s}\\&=0.686\ \mu\mathrm{s}, \end{align*} 时间膨胀因子 \(t/\tau=6.42/0.686=9.4\)。两个计数独立，\(\ln(N_1/N_2)\) 的方差为 \(1/N_1+1/N_2\)（\(\delta\ln N=\delta N/N=1/\sqrt N\)），故 \(\sigma[\ln(N_1/N_2)]=\sqrt{1/563+1/412}=0.0648\)，相对误差 \(0.0648/0.312=21\%\)：\(9.4\pm1.9\)，与 (3) 的 8.1 相差 0.7\(\sigma\)。各计数 \(T_{\rm h}\) 小时，相对误差按 \(1/\sqrt{T_{\rm h}}\) 下降，要降到 10% 需 \(T_{\rm h}=(0.21/0.10)^2\approx4.3\) 小时。
(5) 不是。吸收体相同时，海平面选出的是“到达海平面时”速度为 \(v_1\) 的 \(\mu\) 子，它们在山顶时能量更高，是另一批粒子，其在山顶的计数率并不是 563。要选出同一批，海平面的吸收体要“少吃掉”空气已经吃掉的那部分能量，即 416 MeV：铁的面密度应减少 \(\Delta X_{\rm Fe}=416/1.54\ \mathrm{g/cm^2}={270\,}\mathrm{g/cm^2}\)，即约 34 cm 铁（用 PDG 铁的 CSDA 射程表计算两能量的射程差，得 270 g/cm2）。这一步是“比较同一批粒子的存活率”的关键，也说明为什么不能简单比较两地的全部 \(\mu\) 子计数：全部 \(\mu\) 子的能谱在两地不同。
(6) 取 \(\mu\) 子在某一时刻的瞬时静止系（速度 \(v(x)\) 处的惯性系）。在这个参考系中，大气以速度 \(v\) 迎面而来；地面系中长 \(\dd x\) 的一段空气，在 \(\mu\) 子系中被收缩为 \(\dd x/\gamma\)，它以速度 \(v\) 掠过 \(\mu\) 子所需的时间为 \((\dd x/\gamma)/v=\dd x/(\gamma v)\)——这正是式 (3.5) 中的 \(\dd\tau\)。对整段路程求和，\(\int\dd x/(\gamma v)=\int\dd x/(c\sqrt{\gamma^2-1})\)，与第 (3)(b) 问完全相同。\(\mu\) 子按自己的寿命 \(\tau_\mu\) 衰变，但它需要“经历”的时间只有 \(\tau\approx{0.79\,}\mathrm{\mu s}\)。两种参考系的描述给出同一个可观测量（到达海平面的数目）。
探测器位于胡夫金字塔中轴上、高出地面 21 m 处（王后墓室高度）；金字塔看作底宽 230 m、高 139 m 的正四棱锥，石灰岩 \(\rho={2.2\,}\mathrm{g/cm^3}\)。在探测器北上方有一个空洞，可以看作厚度为 \(t_{\rm v}\)（垂直于空洞自身平面量）、范围足够大的“板”，板面与水平面成 \(\alpha\) 角（\(\alpha>0\) 表示向北升高，\(\alpha<0\) 表示向南升高，大走廊对应 \(\alpha\approx-26^\circ\)）。只考虑南北竖直平面内、天顶角为 \(\theta\)（向北）的视线。在 \(\theta\approx20^\circ\)–\(35^\circ\) 内取 \(\eta'=1.96\)。
证明视线穿过空洞的长度为 \(\ell=t_{\rm v}/\cos(\theta+\alpha)\)，视线方向的岩石总长为 \(L(\theta)=(H_{\rm p}-z_0)/(\cos\theta+\frac{2H_{\rm p}}{B_{\rm p}}\sin\theta)\)（\(H_{\rm p}\)、\(B_{\rm p}\) 为金字塔高与底宽，\(z_0\) 为探测器高度）。
两个方向测得的超额为 \(\varepsilon(20^\circ)=(6.0\pm0.3)\%\)，\(\varepsilon(35^\circ)=(7.7\pm0.35)\%\)。“水平”（\(\alpha=0\)）与“平行于大走廊”（\(\alpha=-26^\circ\)）两种假设各预言比值 \(\varepsilon(35^\circ)/\varepsilon(20^\circ)\) 为多少？数据更支持哪一种，可信度有多高？由数据求 \(t_{\rm v}\)，并说明两个方向给出的 \(t_{\rm v}\) 是否自洽可以作为另一种检验。
某一分析窗口内“无空洞”的期望计数为每天 234 个，空洞造成的超额 \(\varepsilon=9.3\%\)。求只计统计误差时达到 5\(\sigma\) 的天数；若“无空洞”预言有 1% 或 2% 的不相关系统误差，结果如何？
空洞在 \(\tan\theta_x\) 方向的角宽为 \(w\)（顶帽形）。分析窗口宽 \(\Delta\) 取多大时显著性最大？若探测器角分辨不可忽略，最优窗口变宽还是变窄？
两台探测器在东西方向相距 \(b_{\rm d}={10.3\,}\mathrm{m}\)（东边的 A 位于大走廊轴线以东 5.8 m，西边的 B 位于以西 4.5 m）。设它们测得的超额中心分别位于 \(\tan\theta_x=-0.130\) 和 \(+0.100\)（向东为正），每个读数的不确定度为 \(\pm0.0125\)。求空洞中心高出探测器多少、在东西方向的位置，以及高度的误差。
(1) 取南北竖直平面内的坐标（北为正、向上为正）。视线方向 \(\uvec u=(\sin\theta,\cos\theta)\)。向北升高 \(\alpha\) 的板，其法线向南倾 \(\alpha\)：\(\uvec n=(-\sin\alpha,\cos\alpha)\)。穿过厚度 \(t_{\rm v}\) 的板，沿视线走过的长度 \(\ell\) 满足 \(\ell\,(\uvec u\cdot\uvec n)=t_{\rm v}\)，而 \(\uvec u\cdot\uvec n=\cos\theta\cos\alpha-\sin\theta\sin\alpha=\cos(\theta+\alpha)\)，所以 \(\ell=t_{\rm v}/\cos(\theta+\alpha)\)。检验：\(\alpha=0\) 时 \(\ell=t_{\rm v}/\cos\theta\)（斜穿水平板）；视线与板面平行（\(\theta+\alpha=90^\circ\)）时 \(\ell\to\infty\)。北坡面方程为 \(z=H_{\rm p}-(2H_{\rm p}/B_{\rm p})y\)；视线上的点 \((y,z)=(t\sin\theta,\ z_0+t\cos\theta)\)，代入解出 \(t=L(\theta)\) 即得所给公式。数值： \begin{align*} L(20^\circ)&=118/(0.940+1.209\times0.342)\\&=87.2\ \mathrm{m}, \end{align*} \(L(35^\circ)={78.0\,}\mathrm{m}\)。
(2) 由式 (8.4)，\(\varepsilon\approx\eta'\ell/L=\eta't_{\rm v}/[L\cos(\theta+\alpha)]\)，比值中 \(\eta'\) 与 \(t_{\rm v}\) 消去： \begin{align*} &\frac{\varepsilon(35^\circ)}{\varepsilon(20^\circ)}\\ &=\frac{L(20^\circ)}{L(35^\circ)}\cdot\frac{\cos(20^\circ+\alpha)}{\cos(35^\circ+\alpha)}\\ &=\begin{cases}1.118\times1.147=1.28, & \alpha=0,\\ 1.118\times1.007=1.13, & \alpha=-26^\circ.\end{cases} \end{align*} 实测比值 \(7.7/6.0=1.283\)，误差（两个独立相对误差平方相加）\(1.283\sqrt{(0.3/6.0)^2+(0.35/7.7)^2}=0.087\)。它与“水平”假设完全一致，与“平行于大走廊”相差 \((1.283-1.126)/0.087=1.8\sigma\)——倾向于水平，但还不足以排除倾斜。由 \(\varepsilon(20^\circ)\) 求 \(t_{\rm v}\)：水平假设 \begin{align*} t_{\rm v}&=\varepsilon L\cos\theta/\eta'\\&=0.060\times87.2\times0.940/1.96\\&=2.5\ \mathrm{m}, \end{align*} 由 \(35^\circ\) 方向同样得 2.5 m，自洽；倾斜假设由两个方向分别得 2.7 与 3.0 m，相差约 \(1.9\sigma\)（与比值检验基本是同一信息）。2017 年的数据不能区分这两种几何（图 1 中红、绿两种阴影），原因正在于此：两种假设预言的超额剖面差别只有十几个百分点，需要每个方向都有百分之几以下的统计精度。
(3) 由式 (8.5)，只计统计：\(\mu_0=(5/0.093)^2=2.9\times10^3\)，需 \(2.9\times10^3/234=12\) 天。由式 (8.7)，\(1/\mu_0=(\varepsilon/5)^2-\sigma_{\rm sys}^2\)。\(\sigma_{\rm sys}=1\%\)：\(1/\mu_0=3.46\times10^{-4}-1.0\times10^{-4}\)，\(\mu_0=4.1\times10^3\)，需 17 天；\(\sigma_{\rm sys}=2\%\)：\((\varepsilon/5)^2=3.46\times10^{-4}<4\times10^{-4}\)，天花板 \(\varepsilon/\sigma_{\rm sys}=4.6<5\)，再久也达不到 5\(\sigma\)，必须压低系统误差（例如与相邻方向比较、用无目标区域归一化）。
(4) 设单位角宽内本底计数为 \(r\)。窗口 \(\Delta\ge w\) 时，超额计数 \(\varepsilon rw\)，本底 \(r\Delta\)，\(n_\sigma=\varepsilon rw/\sqrt{r\Delta}\propto\Delta^{-1/2}\)；\(\Delta\le w\) 时，超额 \(\varepsilon r\Delta\)，\(n_\sigma=\varepsilon\sqrt{r\Delta}\propto\Delta^{1/2}\)。所以 \(\Delta=w\) 最优，\(n_\sigma^{\max}=\varepsilon\sqrt{rw}\)：窗口过宽，信号被稀释；过窄，丢掉信号。角分辨 \(\sigma_r\) 使空腔像被抹宽，信号溢出到窗口外，最优窗口变宽：\(\sigma_r=0.3w\) 时最优 \(\Delta\approx1.23w\)，最大显著性降为 0.77 倍（图 8b，数值计算）。
(5) 设空腔中心在 A、B 上方高度 \(z\) 处，东西坐标为 \(x_{\rm c}\)（以大走廊轴线为原点，向东为正）。A 在 \(x=+5.8\)，B 在 \(x=-4.5\)： \[\frac{x_{\rm c}-5.8}{z}=-0.130,\qquad \frac{x_{\rm c}+4.5}{z}=+0.100 .\] 两式相减：\(10.3/z=0.230\)，\(z={44.8\,}\mathrm{m}\)；代回得 \(x_{\rm c}=5.8-0.130\times44.8={-0.02\,}\mathrm{m}\)，即在大走廊轴线正上方。高度误差：\(z=b_{\rm d}/\Delta t\)，\(\Delta t=t_B-t_A\) 的误差 \(\sqrt2\times0.0125=0.0177\)， \begin{align*} \delta z&=z\,\delta(\Delta t)/\Delta t\\&=44.8\times0.0177/0.230\\&=3.4\ \mathrm{m}. \end{align*} 文献报告空腔中心距王后墓室地面 40–50 m[9]（本题读数为设定值，用来说明三角测量的原理）。基线越长、角分辨越好，定位越准。
集装箱上下各有一组径迹探测器，测每个 \(\mu\) 子的入射与出射方向（角分辨可忽略），不测动量。下方探测器之下还有一块 10 cm 厚的钢板，只有穿透它的 \(\mu\) 子才计入分析（用来剔除低能本底）。箱内是 1 m 厚的“水等效”货物，箱壁上下合计 0.6 cm 钢（钢均按铁计，密度 7.87 g/cm3）。怀疑货物中藏有一块边长 10 cm 的立方体金属块（铀或铅），只考虑竖直穿过块体中部的 \(\mu\) 子：它们穿过 90 cm 水、10 cm 金属块与上下箱壁。已知 \(X_0\)（cm）：水 36.1，铁 1.757，铅 0.561，铀 0.317；密度：铀 18.95、铅 11.35 g/cm3；能损（MeV cm\(^2\)/g，在相关能区取常数）：水 2.0，铁 1.5，铀 1.15，铅 1.18。入射动量谱用第 9.2 节的模型，表 3 可线性插值。
对 \(p={3\,}\mathrm{GeV/c}\) 的 \(\mu\) 子，分别求穿过“铀块路径”和“铅块路径”时的 \(\psi_0\)。
求两种假设下能被记录的 \(\mu\) 子的最小入射动量 \(p_{\rm c}\)，进而求实测投影角均方根之比 \(R\)。为什么 \(R\) 比第 (1) 问中单一动量的比值小？
同一个 \(\mu\) 子的 \(\psi_x\) 与 \(\psi_y\) 是否独立？求以 5\(\sigma\) 区分“铀”与“铅”所需的 \(\mu\) 子数。
只有从块体上表面进入、从下表面穿出的 \(\mu\) 子才走满 10 cm 厚度。取天顶角 \(\theta<45^\circ\)、积分强度 \(I(\theta)=J_{\rm v}\cos^2\theta\)，其中 \(J_{\rm v}={0.0083\,}\mathrm{cm^{-2}\cdot s^{-1}\cdot sr^{-1}}\) 是第 9.2 节谱模型中 \(p>{0.5\,}\mathrm{GeV/c}\) 的竖直积分强度，求这类 \(\mu\) 子的计数率与所需时间。（提示：借用实验题第 (1) 问的接收度公式。）
多厚的铅块与 10 cm 的铀块散射相同？怎样打破这种“厚度–材料简并”？
走私者用上下各 5 cm 厚的铅壳包住铀块以屏蔽其 \(\gamma\) 射线。这对“发现可疑重物”和“区分铀与铅”分别有利还是不利？
(1) 多层材料先把 \(x/X_0\) 相加：公共部分 \(90/36.1+0.6/1.757=2.493+0.341=2.836\)；铀块路径 \(2.836+31.58=34.42\)；铅块路径 \(2.836+17.82=20.65\)。\(\beta({3\,}\mathrm{GeV/c})=0.9994\)，前因子 \(13.6/(0.9994\times3000)={4.536\,}\mathrm{mrad}\)。由式 (6.5)： \begin{align*} \psi_0^{\rm U}&=4.536\sqrt{34.42}\,(1+0.038\ln34.42)\ \mathrm{mrad}\\ &=30.2\ \mathrm{mrad},\qquad \psi_0^{\rm Pb}=23.0\ \mathrm{mrad}, \end{align*} 单一动量下的比值 \((30.2/23.0)^2=1.73\)。
(2) 被记录的 \(\mu\) 子必须穿透全部材料与下方钢板。按题给能损，公共部分（90 g/cm2 水 + 4.7 g/cm2 箱壁）损失 \(180+7={187\,}\mathrm{MeV}\)；钢板 78.7 g/cm2 损失 118 MeV；铀块 \(189.5\times1.15={218\,}\mathrm{MeV}\)，铅块 \(113.5\times1.18={134\,}\mathrm{MeV}\)。所以入射动能至少为：铀 523 MeV，铅 439 MeV。换成动量 \(p_{\rm c}=\sqrt{(T+m_\mu c^2)^2-(m_\mu c^2)^2}/c\)：铀 0.62 GeV/c，铅 0.53 GeV/c。由表 3 线性插值：铀 \(p_{\rm eff}\approx{1.62\,}\mathrm{GeV/c}\)、\(\kappa\approx9.1\)；铅 \(p_{\rm eff}\approx{1.48\,}\mathrm{GeV/c}\)、\(\kappa\approx10.0\)。由式 (9.1)，均方根为 \begin{align*} &\text{铀：}30.2\times\frac{3.00}{1.62}=56\ \mathrm{mrad},\\ &\text{铅：}23.0\times\frac{3.00}{1.48}=46.5\ \mathrm{mrad},\\ &R=\Big(\frac{56}{46.5}\Big)^2=1.45 . \end{align*} \(R\) 比 1.73 小，是因为铅的能损小、截断低，“铅”假设下被记录的样本里含更多慢 \(\mu\) 子，它们散射大，部分抵消了铅辐射长度较长的效果。两个假设下记录到的不是同一批 \(\mu\) 子，这一点不能忽略。
(3) 不独立：\(\psi_x\)、\(\psi_y\) 共享同一个未知动量（第 9.3 节）。每个 \(\mu\) 子的 \(m=\tfrac12(\psi_x^2+\psi_y^2)\) 的相对方差为 \(2\kappa/3-1\)（式 (9.2)），按“铅”假设取 \(\kappa=10.0\)，得 5.7。由式 (9.3)： \[N\ge5.7\times\Big(\frac{5}{0.45}\Big)^2\approx700 .\] 若误把两投影当独立样本（因子 \((\kappa-1)/2=4.5\)），会低估为约 550。
(4) 实验题第 (1) 问的结果对任何两块相同的正方形平面都成立。这里两块“平面”就是块体的上下表面，取 \(D=d={10\,}\mathrm{cm}\)，\(\tan\theta\le1\) 正好满足 \(D\tan\theta\le d\)。计数率 \begin{align*} \dot N=J_{\rm v}d^2\int_0^{\pi/4}&\cos^2\theta\cdot\sin\theta\cos\theta\\ &\times\big(2\pi-8\tan\theta+2\tan^2\theta\big)\dd\theta . \end{align*} 三项积分分别为 \begin{align*} \int_0^{\pi/4}\cos^3\theta\sin\theta\,\dd\theta&=(1-\cos^445^\circ)/4\\&=0.1875, \end{align*} \(\int_0^{\pi/4}\sin^2\theta\cos^2\theta\,\dd\theta=\pi/32=0.0982\)，\(\int_0^{\pi/4}\sin^3\theta\cos\theta\,\dd\theta=\sin^445^\circ/4=0.0625\)，括号积分为 \begin{align*} &2\pi\times0.1875-8\times0.0982+2\times0.0625\\&=0.518. \end{align*} 于是 \begin{align*} \dot N&=0.0083\times100\times0.518\ \mathrm{s^{-1}}\\&=0.43\ \mathrm{s^{-1}}\\&=26\ \text{个/分钟}. \end{align*} 若把从上表面进入的 \(\mu\) 子全算上（只积第一项），会得到 59 个/分钟——其中 56% 从侧面穿出，没有走满厚度。所需时间 \(700/26\approx27\) 分钟。（按截断严格说铅假设应取 \(p>{0.53\,}\mathrm{GeV/c}\) 的强度，比 \(J_{\rm v}\) 小约 1%，可忽略。第 (1)–(3) 问都按竖直路径计算；走满厚度的 \(\mu\) 子按接收度加权，\(\avg{\sec\theta}=1.09\)，块体中的路程、\(x/X_0\) 与能损平均约大 9%，对 \(R\) 的影响很小，这里只作估计。）
(5) 散射只依赖 \(x/X_0\)：\(x_{\rm Pb}=10\times0.561/0.317={17.7\,}\mathrm{cm}\) 的铅与 10 cm 铀一样。打破简并的办法：(i) 用最近点法做三维成像，直接测出块体厚度（cm 级分辨即可区分 10 与 18 cm）；(ii) 比较不同入射角的 \(\mu\) 子——斜穿时路径变长，按已知几何拟合可同时得到尺寸与 \(X_0\)；(iii) 利用能损：铀块使截断更高，被“吃掉”的 \(\mu\) 子比例不同（第 (2) 问）。从 \(X_0\) 再推 \(Z\)：铀 \(X_0\)（cm）比铅小 44%，而 \(Z\) 只差 12%，因为以 cm 计的 \(X_0\) 还含密度（19.0 对 11.3 g/cm3）。
(6) 对“发现”有利：铅壳使穿过块体的 \(x/X_0\) 再增加 \(2\times5/0.561=17.8\)，散射更强，可疑重物更显眼。对“区分”不利：铅壳给两种假设加上同一个 \(x/X_0\)，单一动量下 \(R\) 从 1.73 降到 \(\big[\psi_0(34.42+17.8)/\psi_0(20.65+17.8)\big]^2=1.39\)，所需 \(\mu\) 子数增加约 \((0.73/0.39)^2\approx3.5\) 倍；而且“铅壳 + 铀”可被一块更大的纯铅冒充。所以实际检查先用散射找出“高 \(Z\) 重物”，再用其他手段（几何成像、诱发裂变的中子或 \(\gamma\) 测量等）确认成分。
上面的解答用入射动量计算 \(p_{\rm eff}\)，没有计入 \(\mu\) 子在块体与货物中的减速。脚本 ex3_numbers.py 用 PDG 射程表逐层积分减速，并用 Highland 公式的逐层形式 \(\psi_0^2\propto\int\dd x/[X_0(\beta pc)^2]\) 计算散射。块体放在货物中部（上下各 45 cm 水）时，得到 \(R=1.58\)，\(\kappa_{\rm Pb}=12.4\)，\(N\approx540\)，约 21 分钟，与上面的估算（27 分钟）同量级。若去掉下方钢板，最慢的 \(\mu\) 子也被记录，它们在路径末端散射极大，\(\kappa\) 升到约 20，所需时间增到约 47 分钟。慢 \(\mu\) 子的散射集中在路径末端，所以结果还与块体的竖直位置有关（脚本 ex3_position.py）：块体在货物顶部、中部、底部时，有钢板分别约为 19、21、25 分钟，无钢板约为 39、47、94 分钟。这说明剔除（或测出）慢 \(\mu\) 子对散射成像很重要。
探测器在火山脚下，视线仰角 \(10^\circ\)（天顶角 \(\theta=80^\circ\)），由地形图得这条视线在山体中的长度 \(L={500\,}\mathrm{m}\)。探测器对该方向的有效“面积 × 立体角”为 \(S_{\rm eff}\Delta\Omega={1.0\,}\mathrm{m^2}\times{4.0\times 10^{-3}\,}\mathrm{sr}\)，曝光 60 天，记录到 \(N=260\) 个穿山 \(\mu\) 子。在 \(\theta=80^\circ\)、\(E\approx{350\,}\mathrm{GeV}\) 附近，\(I(>E)=I_1(E/E_1)^{-\eta}\)，\(I_1={1.176\times 10^{-2}\,}\mathrm{m^{-2}\cdot s^{-1}\cdot sr^{-1}}\)，\(E_1={350\,}\mathrm{GeV}\)，\(\eta=1.99\)。能损参数同上。
求视线上岩石的平均密度。
求密度的统计误差。
怀疑“穿山 \(\mu\) 子”中混有与岩石厚度无关的本底（低能电子、被散射进来的 \(\mu\) 子、从背后入射被误判方向的粒子等），其强度 \(\Phi_{\rm bg}\) 未知。同一台探测器在同样的曝光中还记录了另一条视线：天顶角也约 \(80^\circ\)，山体中长 \(L_B={1000\,}\mathrm{m}\)，计数 \(N_B=85\)（\(S_{\rm eff}\Delta\Omega\) 相同）。在 \(E\approx{800\,}\mathrm{GeV}\) 附近 \(I(>E)=I_2(E/E_2)^{-\eta_2}\)，\(I_2={2.127\times 10^{-3}\,}\mathrm{m^{-2}\cdot s^{-1}\cdot sr^{-1}}\)，\(E_2={800\,}\mathrm{GeV}\)，\(\eta_2=2.14\)。假设两条视线上岩石平均密度相同、单位立体角的本底强度相同，求 \(\bar\rho\) 与 \(\Phi_{\rm bg}\)，以及 \(\bar\rho\) 的统计误差。与第 (1) 问比较，并推出“本底比例 \(f=\Phi_{\rm bg}/\Phi_{\rm sig}\) 造成的密度偏差”的一般公式。
若火山口下方这条视线上有一段 50 m 长的岩石（2.4 g/cm3）被平均密度 1.0 g/cm3 的含气泡岩浆取代，计数率变化多少？若想在一周内以 3\(\sigma\) 看出这一变化，\(S_{\rm eff}\Delta\Omega\) 至少要增大多少倍？
(1) 透射强度 \begin{align*} \Phi&=N/(S_{\rm eff}\Delta\Omega\,T)\\&=260/(1.0\times4.0\times10^{-3}\times5.184\times10^{6}\ \mathrm{s})\\&=1.254\times10^{-2}\ \mathrm{m^{-2}s^{-1}sr^{-1}}. \end{align*} 由 \(\Phi=I_1(E_{\min}/E_1)^{-\eta}\)： \begin{align*} E_{\min}&=350\times(1.254/1.176)^{-1/1.99}\ \mathrm{GeV}\\&=339\ \mathrm{GeV}. \end{align*} 由式 (4.8) 反解： \begin{align*} X&=\frac1b\ln(1+bE_{\min}/a)\\&=\frac{\ln(1+339/550)}{4.0\times10^{-6}\ \mathrm{cm^2/g}}\\&=\frac{0.480}{4.0\times10^{-6}}\ \mathrm{g/cm^2}\\&=1.200\times10^{5}\ \mathrm{g/cm^2}. \end{align*} 平均密度 \begin{align*} \bar\rho&=X/L=1.200\times10^{5}/(5.0\times10^{4})\ \mathrm{g/cm^3}\\&=2.40\ \mathrm{g/cm^3}. \end{align*} （用完整谱数值反演同样得 2.40。）
(2) \(bX=0.480\)，\(g=0.480\,\ee^{0.480}/(\ee^{0.480}-1)=1.257\)，\(\eta'=1.99\times1.257=2.50\)。由式 (8.3)，\(\delta\bar\rho/\bar\rho=1/(2.50\sqrt{260})=2.5\%\)，即 \(\bar\rho=(2.40\pm0.06)\,\mathrm{g/cm^3}\)。注意 \(\eta'\) 比金字塔大：辐射损失使 \(E_{\min}\) 随 \(X\) 增长得比线性快，密度测量反而更灵敏。
(3) 两个观测量（\(N_A=260\)、\(N_B=85\)）、两个未知量（\(\bar\rho\)、\(\Phi_{\rm bg}\)），可以解出。实测强度 \(\Phi_A={1.254\times 10^{-2}\,}\mathrm{}\)、 \begin{align*} \Phi_B&=85/(4.0\times10^{-3}\times5.184\times10^{6})\\&=4.10\times10^{-3}\ \mathrm{m^{-2}s^{-1}sr^{-1}}. \end{align*} 两式相减消去本底： \begin{align*} \Phi_A-\Phi_B&=\Phi_{\rm sig}(\bar\rho L_A)-\Phi_{\rm sig}(\bar\rho L_B)\\ &=8.44\times10^{-3}\ \mathrm{m^{-2}s^{-1}sr^{-1}}, \end{align*} 右边是 \(\bar\rho\) 的单调减函数，试几个值（\(E_{\min}\) 由式 (4.8)，强度由两段幂律）：
| \(\bar\rho\) | \(E_{\min,A}\)/GeV | \(E_{\min,B}\)/GeV | \(\Phi_{{\rm sig},A}\) | \(\Phi_{{\rm sig},B}\) | 差 |
|---|---|---|---|---|---|
| 2.40 | 339 | 886 | \(12.54\times10^{-3}\) | \(1.71\times10^{-3}\) | \(10.84\times10^{-3}\) |
| 2.60 | 375 | 1006 | \(10.25\times10^{-3}\) | \(1.30\times10^{-3}\) | \(8.95\times10^{-3}\) |
| 2.70 | 394 | 1070 | \(9.31\times10^{-3}\) | \(1.14\times10^{-3}\) | \(8.16\times10^{-3}\) |
线性插值得 \(\bar\rho=2.66\)（精确求根 2.663，用完整谱 2.663）， \begin{align*} \Phi_{\rm bg}&=\Phi_B-\Phi_{{\rm sig},B}\\&=4.10\times10^{-3}-1.20\times10^{-3}\\&=2.90\times10^{-3}\ \mathrm{m^{-2}s^{-1}sr^{-1}}, \end{align*} 占视线 A 真实透射强度的 30%。第 (1) 问忽略本底，密度偏低 0.26。误差：\(\Phi_A-\Phi_B\) 的统计误差 \(\sqrt{N_A+N_B}/(S_{\rm eff}\Delta\Omega\,T)={0.90\times 10^{-3}\,}\mathrm{}\)，除以上表中“差”对 \(\bar\rho\) 的斜率（约 7.7 × 10−3 每 g/cm\(^3\)），得 \(\delta\bar\rho\approx0.12\)，即 \(\bar\rho=(2.66\pm0.12)\,\mathrm{g/cm^3}\)。扣本底的代价是误差从 2.5% 增大到 4.4%：长视线 B 的计数主要用来测本底。这一解法依赖“两条视线密度相同”的假设；若不同，还需要第三个观测量（例如探测器背向山体、朝向天空下方的方向，那里只有本底）。
一般地，实测 \(\Phi_{\rm m}=\Phi_{\rm sig}(1+f)\)，\(\ln\Phi\) 偏高 \(\ln(1+f)\)，由 \(\dd\ln\Phi=-\eta'\dd\ln X\) 得 \[\frac{\delta\bar\rho}{\bar\rho}\approx-\frac{\ln(1+f)}{\eta'}=-\frac{\ln1.30}{2.50}=-10\%,\] 与上面 \(2.40\) 对 \(2.66\)（\(-9.8\%\)）一致。本底总是让密度偏低；视线越长（透射越弱），同样绝对强度的本底所占比例 \(f\) 越大，偏差越严重。
(4) 按第 (1) 问的密度 2.40 计算（忽略本底）。 \begin{align*} \Delta X&=-(2.4-1.0)\times5000\\&=-7.0\times10^{3}\ \mathrm{g/cm^2}, \end{align*} \(\Delta X/X=-5.8\%\)。一阶估计 \(\varepsilon\approx\eta'\times5.8\%=14.6\%\)；按幂律严格计算，\(E_{\min}\) 从 339 降到 314 GeV，\(\varepsilon=(314/339)^{-1.99}-1=16.1\%\)（\(\Delta X/X\) 已不算很小，一阶估计偏低 1.5 个百分点）。一周内 3\(\sigma\)：若“变化前”的计数率由长期测量精确知道，需 \(\mu_0=(3/0.161)^2=347\) 个/周，即 50 个/天；现在只有 \(260/60=4.3\) 个/天，\(S_{\rm eff}\Delta\Omega\) 需增大约 11 倍。若基线也只测了一周，两周的计数之差的方差是单周的 2 倍（第 5.2 节注 (i)），需增大约 25 倍（差值的方差为 \(\mu_0(2+\varepsilon)\)；按 \(2\mu_0\) 近似为 23 倍）；或者把时间分辨率放宽到几个月。这说明“动态 \(\mu\) 子成像”要求的探测器规模远大于“静态”。
11 实验与数据处理：\(\mu\) 子的天顶角分布
两块相同的正方形位置灵敏探测平面（边长 \(d={1.00\,}\mathrm{m}\)）水平放置，上下正对，间距 \(D={0.50\,}\mathrm{m}\)（图 11）。每个 \(\mu\) 子在两平面上的击中点决定其方向，按天顶角 \(\theta\) 每 \(5^\circ\) 分一组，计数 \(T={600\,}\mathrm{s}\)。假设 \(\mu\) 子强度 \(I(\theta)=I_0\cos^n\theta\)，与方位角无关。表 4 给出各组计数 \(N_i\) 与主办方按组中心 \(\theta_i\) 算出的接收度 \(G_i\)（第 8.1 节：该组所有方向的 \(\int S_{\rm eff}\,\dd\Omega\)，期望计数 \(\approx TG_iI(\theta_i)\)）。本题数据是按题述模型（\(I_0\)、\(n\) 取某组值）计算期望后，按泊松分布抽样生成的模拟数据，随机种子 20261017 在出题前固定。
一个方向为 \((\theta,\varphi)\) 的 \(\mu\) 子要被记录，必须同时穿过两块平面。证明对该方向的有效面积为 \begin{align*} S_{\rm eff}={}&(d-D\tan\theta\abs{\cos\varphi})\\&\times(d-D\tan\theta\abs{\sin\varphi})\cos\theta \end{align*} （\(D\tan\theta\le d\)），并由此证明第 \(i\) 组的期望计数为 \(\mu_i=T\int_{\text{组}i}I(\theta)\,g(\theta)\,\dd\theta\)， \begin{align*} g(\theta)={}&\sin\theta\cos\theta\\&\times(2\pi d^2-8dD\tan\theta\\&\quad+2D^2\tan^2\theta). \end{align*} 用组中心近似验证表中 \(30^\circ\)–\(35^\circ\) 一组的 \(G_i\)。
把模型线性化，说明为什么拟合时第 \(i\) 个点的权重应取 \(w_i=N_i\)，用加权最小二乘求 \(n\) 与 \(I_0\) 及其标准误差，并计算 \(\chi^2\)。
某同学只除以立体角与投影面积 \(2\pi d^2\sin\theta\cos\theta\,\Delta\theta\)，忘记了两平面重叠面积随 \(\theta\) 的缩小。他会得到多大的 \(n\)？偏差方向如何理解？
讨论系统误差（至少三项），说明它们使 \(n\) 偏大还是偏小，并提出改进方案。
| \(\theta\) 范围 | \(N_i\) | \(G_i\) | \(x_i\) | \(y_i\) |
|---|---|---|---|---|
| \(0^\circ\)–\(5^\circ\) | 925 | 23.23 | \(-0.0010\) | 4.195 |
| \(5^\circ\)–\(10^\circ\) | 2485 | 65.11 | \(-0.0086\) | 4.153 |
| \(10^\circ\)–\(15^\circ\) | 3665 | 99.96 | \(-0.0240\) | 4.113 |
| \(15^\circ\)–\(20^\circ\) | 4416 | 126.93 | \(-0.0474\) | 4.060 |
| \(20^\circ\)–\(25^\circ\) | 4738 | 145.38 | \(-0.0792\) | 3.995 |
| \(25^\circ\)–\(30^\circ\) | 4771 | 154.99 | \(-0.1199\) | 3.938 |
| \(30^\circ\)–\(35^\circ\) | 4105 | 155.72 | \(-0.1703\) | 3.783 |
| \(35^\circ\)–\(40^\circ\) | 3660 | 147.86 | \(-0.2315\) | 3.720 |
| \(40^\circ\)–\(45^\circ\) | 2771 | 132.04 | \(-0.3048\) | 3.555 |
| \(45^\circ\)–\(50^\circ\) | 1892 | 109.25 | \(-0.3922\) | 3.363 |
| \(50^\circ\)–\(55^\circ\) | 1110 | 80.90 | \(-0.4963\) | 3.130 |
| \(55^\circ\)–\(60^\circ\) | 507 | 48.89 | \(-0.6212\) | 2.850 |
(1) 方向 \((\theta,\varphi)\) 的径迹从上平面上的点 \((x,y)\) 落到下平面上的点 \((x-D\tan\theta\cos\varphi,\ y-D\tan\theta\sin\varphi)\)。被记录要求两个点都在各自的正方形内：在 \(x\) 方向，上平面的击中点要落在一个长 \(d-D\tan\theta\abs{\cos\varphi}\) 的区间内，\(y\) 方向同理。所以“能被记录的上平面击中点”的集合是一个矩形，面积 \((d-D\tan\theta\abs{\cos\varphi})(d-D\tan\theta\abs{\sin\varphi})\)（图 11b）。这是水平面积；垂直于径迹的面积要再乘 \(\cos\theta\)，得到 \(S_{\rm eff}\)。由强度定义，单位时间计数 \(\int I\,S_{\rm eff}\,\dd\Omega\)，\(\dd\Omega=\sin\theta\,\dd\theta\,\dd\varphi\)。对 \(\varphi\) 积分要用（附录 D） \begin{align*} &\int_0^{2\pi}\abs{\cos\varphi}\dd\varphi=\int_0^{2\pi}\abs{\sin\varphi}\dd\varphi=4,\\ &\int_0^{2\pi}\abs{\sin\varphi\cos\varphi}\dd\varphi=2, \end{align*} 得 \(\int_0^{2\pi}(\cdots)\dd\varphi=2\pi d^2-8dD\tan\theta+2D^2\tan^2\theta\)，即题给的 \(g(\theta)\)。组中心近似：\(G_i=g(\theta_i)\Delta\theta\)，\(\Delta\theta=5^\circ={0.08727\,}\mathrm{rad}\)。对 \(\theta_i=32.5^\circ\)：\(\tan\theta_i=0.6371\)，\(\sin\theta_i\cos\theta_i=0.4532\)，括号内 \(6.2832-2.5484+0.2029=3.9377\ \mathrm{m^2}\)， \begin{align*} G&=0.4532\times3.9377\times0.08727\\&=0.1557\ \mathrm{m^2 sr}, \end{align*} 与表值一致。（要求 \(D\tan\theta\le d\)，即 \(\theta\le63.4^\circ\)，本题的组都满足。）
(2) 期望计数 \(\mu_i\approx TG_iI_0\cos^n\theta_i\)，取对数： \begin{align*} y_i&\equiv\ln\frac{N_i}{TG_i}=\ln I_0+n\,x_i,\\ x_i&=\ln\cos\theta_i . \end{align*} 这是关于 \((\ln I_0,n)\) 的直线。\(y_i\) 的误差来自 \(N_i\) 的泊松涨落：\(\delta y_i=\delta N_i/N_i=1/\sqrt{N_i}\)，所以各点精度差别很大（\(N\) 从 507 到 4771），应以方差的倒数 \(w_i=1/\delta y_i^2=N_i\) 为权（附录 C 给出加权最小二乘的推导）。记 \([f]\equiv\sum_iw_if_i\)，由表中数据：\([1]=35045\)，\([x]=-5265.4\)，\([y]=135190.2\)，\([x^2]=1443.20\)，\([xy]=-18969.95\)，\(\Delta=[1][x^2]-[x]^2={2.2852\times 10^{7}\,}\mathrm{}\)。 \begin{align*} n&=\frac{[1][xy]-[x][y]}{\Delta}=2.058,\\ \sigma_n&=\sqrt{\frac{[1]}{\Delta}}=0.039,\\ \ln I_0&=\frac{[x^2][y]-[x][xy]}{\Delta}=4.1669,\\ \sigma_{\ln I_0}&=\sqrt{\frac{[x^2]}{\Delta}}=0.0079 . \end{align*} 即 \(n=2.06\pm0.04\)，\(I_0=(64.5\pm0.5)\ \mathrm{m^{-2}s^{-1}sr^{-1}}\)（误差只含统计）。\(\chi^2=\sum w_i(y_i-\ln I_0-nx_i)^2=12.7\)，自由度 \(12-2=10\)，\(\chi^2/\nu=1.27\)，拟合可以接受（\(\chi^2>12.7\) 的概率约 24%），不必放大误差。作为核对，用同一模型重复模拟 4000 次实验，拟合的 \(n\) 的标准差为 0.039，与公式一致。\(n\approx2\) 与 PDG 的 \(\cos^2\theta\) 相符[16]（图 12）。
(3) 他用 \(G_i'=2\pi d^2\sin\theta_i\cos\theta_i\Delta\theta\)，比真实接收度多了因子 \(2\pi d^2/(2\pi d^2-8dD\tan\theta+2D^2\tan^2\theta)\)，这个因子随 \(\theta\) 增大而增大，使 \(y_i'\) 在大角度处被压低，直线变陡。重新拟合得 \(n'=4.26\)，比正确值大 2.2——几何接收度在这里比统计误差重要 50 倍。若连 \(\cos\theta\) 投影也忘了（只除以立体角），还要再多出 1，得 \(n''=5.26\)（式 (7.6) 下的讨论）。重叠因子 \(1-\frac{4}{\pi}\frac{D}{d}\tan\theta+\frac{1}{\pi}\frac{D^2}{d^2}\tan^2\theta\) 在小角度处随 \(\theta\) 线性下降，而 \(\ln\cos\theta\approx-\theta^2/2\) 按平方变化，所以它不等价于任何固定幂次的 \(\cos\theta\)：错误版本的点不在一条直线上（图 12 中红点略呈弯曲），拟合斜率还依赖于所用的角度范围。
(4) 系统误差：
组中心近似：\(\cos^n\theta\) 与 \(g(\theta)\) 都在组内变化，用中心值代替组内积分有偏差。用无噪声的精确期望值做同样拟合，得 \(n=2.075\)，而模拟所用真值为 2.10：此近似使 \(n\) 偏小约 0.025（0.6\(\sigma\)）。改进：拟合时对每组做数值积分（或把组分细）。
角分辨与组间迁移：击中点分辨 1 cm 时，角分辨约 \(\sqrt2\times{1\,}\mathrm{cm}/{50\,}\mathrm{cm}\approx{0.03\,}\mathrm{rad}\approx1.6^\circ\)，与组宽 \(5^\circ\) 可比。真实分布在大角度处随 \(\theta\) 下降，事例从计数多的组“漏”到相邻组，平滑了分布，使 \(n\) 偏小（蒙特卡罗，脚本
exp_smear.py：每个击中点坐标加 1 cm 的高斯抹宽，按同样方法拟合，\(n\) 偏小约 0.025，与统计误差同量级）。改进：增大 \(D\)（代价是大角度接收度下降）或用模拟做反卷积。能量阈值：屋顶、探测器材料吸收低能 \(\mu\) 子；低能 \(\mu\) 子的角分布更陡，所以阈值越高，\(n\) 越小。结论“\(n\approx2\)”只对本装置的阈值成立。改进：在两平面间加已知厚度的铅板，测 \(n\) 随阈值的变化。
非 \(\mu\) 子本底与效率：低能电子（软成分）、偶然符合会增加计数；若它们在大角度处所占比例更大，大角度计数被抬高，\(n\) 偏小。斜穿时粒子在闪烁体中的路程长、信号大，探测效率随 \(\theta\) 增大，同样使 \(n\) 偏小。改进：加铅吸收电子，用三平面符合，测各角度的效率。
大气条件：气压升高，\(\mu\) 子穿过的空气更多，计数下降；斜向路径更长，大角度方向降得更多，使 \(n\) 偏大。测量期间应记录气压，必要时分时段修正。
12 小结
本讲的推导链如下（箭头表示“由……推出”）：
固有时 \(\dd\tau=\dd t/\gamma=m\,\dd x/p\)（式 (3.5)）+ 无记忆衰变 \(\Rightarrow\) 存活概率 \(\ee^{-\tau/\tau_\mu}\)（式 (3.6)）\(\Rightarrow\) \(\mu\) 子能从 15 km 高空到达地面；能量沿途变化时逐段积分（式 (7.10)，例题 1 的 \(\operatorname{arcosh}\) 公式）。
一次擦肩而过的冲量 \(2ze^2/(4\pi\varepsilon_0bv)\)（式 (4.2)）\(\Rightarrow\) 对碰撞参数积分，得阻止本领正比于 \(Z/A\)、反比于 \(\beta^2\)，再乘一个对数（式 (4.4)）\(\Rightarrow\) Bethe 公式与最小电离 \(\approx{2\,}\mathrm{MeV\cdot cm^2/g}\) \(\Rightarrow\) 加上辐射损失得 \(a+bE\) \(\Rightarrow\) \(E_{\min}(X)=(a/b)(\ee^{bX}-1)\)（式 (4.8)）。
初级谱 \(E^{-2.7}\) + \(\pi\) 的衰变/碰撞竞争（\(E_\pi^*\approx{113\,}\mathrm{GeV}\)，式 (7.4)）\(\Rightarrow\) \(\mu\) 子谱与角分布 \(\Rightarrow\) 积分谱斜率 \(\eta\)。
\(\Phi(X)=I(>E_{\min}(X))\) \(\Rightarrow\) 灵敏度 \(\eta'=\eta\,g(bX)\approx2\)（式 (8.2)）\(\Rightarrow\) 空腔超额 \(\varepsilon\approx\eta'\ell/L\)（式 (8.4)）与密度精度 \(1/(\eta'\sqrt N)\)（式 (8.3)）。
泊松方差 \(=\) 均值（式 (5.1)）\(\Rightarrow\) \(n_\sigma=\varepsilon\sqrt{\mu_0}\)（式 (8.5)）\(\Rightarrow\) 曝光时间（式 (8.6)）；系统误差给出天花板 \(\varepsilon/\sigma_{\rm sys}\)（式 (8.7)）。
卢瑟福小角散射 + 随机行走 \(\Rightarrow\) \(\avg{\psi^2}\propto Z(Z+1)\rho x/(A\,\beta^2p^2)\)（式 (6.2)）\(\Rightarrow\) 与轫致辐射同源，用 \(X_0\) 表示 \(\Rightarrow\) Highland 公式（式 (6.5)）\(\Rightarrow\) 动量弥散使峰度增大（式 (9.1)）\(\Rightarrow\) 所需 \(\mu\) 子数（式 (9.3)）。
贯穿全讲的两个“对数”值得记住：能量损失与多次散射都来自对碰撞参数的 \(\int\dd b/b\)；透射成像的灵敏度来自积分谱的对数斜率。决定一次 \(\mu\) 子成像成败的，往往不是物理原理，而是“每个像素能攒到多少个 \(\mu\) 子”和“对无目标模型的把握有多好”。
参考文献
附录 A 能量损失公式中的对数
为什么 \(b_{\min}\) 由量子力学决定
经典截断由 \(\Delta E(b_{\min})=W_{\max}\) 给出：令式 (4.3) 等于 \(2m_{\rm e}c^2\beta^2\gamma^2\)，得 \(b_{\min}^{\rm cl}=zr_{\rm e}/(\beta^2\gamma)\)。但电子是量子粒子：在重粒子参考系中电子以速度 \(v\)、动量约 \(\gamma m_{\rm e}v\) 迎面而来，由不确定关系，碰撞参数不可能被确定到优于 \(\hbar/\Delta p_{\max}\)，其中 \(\Delta p_{\max}=2\gamma m_{\rm e}v\) 是正碰时的最大动量转移，即量子截断 \(b_{\min}^{\rm q}\approx\hbar/(2\gamma m_{\rm e}v)\)。两者之比 \[\frac{b_{\min}^{\rm q}}{b_{\min}^{\rm cl}}=\frac{\hbar/(2\gamma m_{\rm e}v)}{zr_{\rm e}/(\beta^2\gamma)}=\frac{\hbar c\,\beta}{2z\,r_{\rm e}m_{\rm e}c^2}=\frac{\beta}{2z\alpha},\] 用到 \(r_{\rm e}m_{\rm e}c^2=\alpha\hbar c\)。对 \(\mu\) 子（\(z=1\)，\(\beta\approx1\)）比值约 69，所以真正的截断是较大的量子值。只有对慢速重离子（\(z\alpha/\beta\gtrsim1\)）玻尔的经典截断才适用。
\(b_{\max}\) 有多大
以标准岩石 \(\bar I=\hbar\bar\omega={136\,}\mathrm{eV}\) 为例，\(\bar\omega={2.07\times 10^{17}\,}\mathrm{s^{-1}}\)；\(\beta\gamma=3.5\) 时 \begin{align*} b_{\max}&=\gamma v/\bar\omega\\&=3.5\times3\times10^{8}\ \mathrm{m/s}/(2.07\times10^{17}\ \mathrm{s^{-1}})\\&\approx5\ \mathrm{nm}, \end{align*} 远大于原子间距（约 0.2 nm）。在这么大的范围内，粒子的电场会使沿途的大量原子极化，这些极化电荷屏蔽了远处的场，所以 \(\gamma\) 很大时 \(b_{\max}\) 不能无限增长——这就是密度效应修正 \(\delta\) 的物理来源，它在固体、液体中比在气体中重要（图 3 中水、岩石的曲线在高能处比“不含密度效应”的点线平坦）。
与 Bethe 公式对照
\begin{align*} \ln(b_{\max}/b_{\min}^{\rm q})&=\ln\big[(\gamma v/\bar\omega)(2\gamma m_{\rm e}v/\hbar)\big]\\&=\ln(2m_{\rm e}c^2\beta^2\gamma^2/\bar I). \end{align*} 式 (4.5) 方括号中 \(\tfrac12\ln(2m_{\rm e}c^2\beta^2\gamma^2W_{\max}/\bar I^2)\) 在 \(W_{\max}=2m_{\rm e}c^2\beta^2\gamma^2\) 时正好等于这个对数。\(-\beta^2\) 一项来自相对论运动学对大能量转移截面的修正；它们都只改变几个百分点。
附录 B 多次散射：从卢瑟福截面到 \(E_{\rm s}={21\,}\mathrm{MeV}\)
第一步：对数
第 6.2 节得到（\(z=1\)，包括电子贡献） \begin{align*} \frac{\dd\avg{\psi_{\rm sp}^2}}{\dd x}&=8\pi n_{\rm A}Z(Z+1)r_{\rm e}^2\Big(\frac{m_{\rm e}c^2}{\beta pc}\Big)^2\Lambda_{\rm s},\\ \Lambda_{\rm s}&=\ln\frac{b_{\max}}{b_{\min}} . \end{align*} \(b_{\max}\)：大于原子尺度时，核电荷被电子云屏蔽，取 Thomas–Fermi 屏蔽半径 \(a_{\rm TF}=0.885a_0Z^{-1/3}\)（\(a_0={0.529\times 10^{-10}\,}\mathrm{m}\) 为玻尔半径）。\(b_{\min}\)：小于原子核半径 \(R_{\rm N}\approx1.2A^{1/3}\,\mathrm{fm}\) 时，\(\mu\) 子穿进核内，所受力不再是点电荷的库仑力；对 GeV 量级的 \(\mu\) 子，其约化德布罗意波长 \(\hbar/p\approx{0.07\,}\mathrm{fm}\) 远小于 \(R_{\rm N}\)，量子截断不起作用。
第二步：与辐射长度比较
由式 (6.3)，\(1/X_0=4\alpha r_{\rm e}^2n_{\rm A}Z(Z+1)\Lambda_{\rm r}\)，\(\Lambda_{\rm r}=\ln(184Z^{-1/3})\)（这里 \(X_0\) 以长度计，\(n_{\rm A}\) 为原子数密度）。相除： \[\frac{\avg{\psi_{\rm sp}^2}}{x/X_0}=\frac{8\pi}{4\alpha}\Big(\frac{m_{\rm e}c^2}{\beta pc}\Big)^2\frac{\Lambda_{\rm s}}{\Lambda_{\rm r}}=\frac{2\pi}{\alpha}\Big(\frac{m_{\rm e}c^2}{\beta pc}\Big)^2\frac{\Lambda_{\rm s}}{\Lambda_{\rm r}} .\] 材料常数 \(n_{\rm A}Z(Z+1)\) 完全消去，只剩两个对数之比。
第三步：对数之比约为 2
数值（verify_numbers.py）：
| 材料 | \(a_{\rm TF}/R_{\rm N}\) | \(\Lambda_{\rm s}\) | \(\Lambda_{\rm r}\) | \(\Lambda_{\rm s}/\Lambda_{\rm r}\) |
|---|---|---|---|---|
| 岩石（\(Z=11\), \(A=22\)） | 6260 | 8.74 | 4.42 | 1.98 |
| 铁 | 3440 | 8.14 | 4.13 | 1.97 |
| 铅 | 1520 | 7.32 | 3.75 | 1.95 |
| 铀 | 1390 | 7.24 | 3.71 | 1.95 |
比值几乎与材料无关，约为 2。代入得 \(\avg{\psi_{\rm sp}^2}\approx(4\pi/\alpha)(m_{\rm e}c^2/\beta pc)^2\,x/X_0\)，即式 (6.4)， \begin{align*} E_{\rm s}&=m_{\rm e}c^2\sqrt{4\pi/\alpha}\\&=0.511\times41.5\ \mathrm{MeV}\\&=21.2\ \mathrm{MeV}. \end{align*} 这个“2”是一个数值事实，而不是结构上的必然：散射对数的下截断是核半径，轫致辐射对数 \(\ln(184Z^{-1/3})\) 的下截断则是电子的约化康普顿波长 \(\hbar/(m_{\rm e}c)={386\,}\mathrm{fm}\)；两个对数之比在 \(Z=11\)–92 内恰好都落在 1.95–1.98 之间，Rossi 与 Greisen 正是利用了这一点。
第四步：为什么 Highland 系数是 13.6 而不是 15.0
投影角的“全体均方根”包含了少数大角单次散射（卢瑟福截面 \(\propto1/\psi^4\) 的长尾）。用高斯拟合分布的中央 98%，得到的宽度比全体均方根小；Highland 与 Lynch–Dahl 对精确 Molière 分布拟合，得到系数 13.6 MeV 和随 \(x/X_0\) 缓慢增长的对数修正[21], [22]。
附录 C 统计工具的推导
泊松分布
沿用第 5.1 节的记号，\(N\) 是 \(n_{\rm s}\) 个独立伯努利变量之和，服从二项分布 \(\Pr(N)=\binom{n_{\rm s}}{N}q^N(1-q)^{n_{\rm s}-N}\)。令 \(n_{\rm s}\to\infty\)、\(n_{\rm s}q=\mu\) 固定、\(N\) 固定：\(\binom{n_{\rm s}}{N} q^N=\frac{n_{\rm s}(n_{\rm s}-1)\cdots(n_{\rm s}-N+1)}{N!}\frac{\mu^N}{n_{\rm s}^N}\to\frac{\mu^N}{N!}\)，\((1-q)^{n_{\rm s}-N}=(1-\mu/n_{\rm s})^{n_{\rm s}-N}\to\ee^{-\mu}\)。所以 \(\Pr(N)=\mu^N\ee^{-\mu}/N!\)。直接验证：\(\sum_NN\Pr(N)=\mu\)，\(\sum_NN(N-1)\Pr(N)=\mu^2\)，故 \(\operatorname{Var}N=\mu^2+\mu-\mu^2=\mu\)。
高斯分布的四阶矩
对 \(f(\psi)=\ee^{-\psi^2/2\sigma^2}/(\sqrt{2\pi}\sigma)\)，分部积分 \(\int\psi^4f\,\dd\psi=3\sigma^2\int\psi^2f\,\dd\psi=3\sigma^4\)，所以 \(\kappa=3\)。
混合分布的峰度
若 \(\psi\) 在给定 \(\sigma\) 时为高斯分布，而 \(\sigma\) 本身随机，则 \(\avg{\psi^2}=\avg{\sigma^2}\)，\(\avg{\psi^4}=3\avg{\sigma^4}\)，\(\kappa=3\avg{\sigma^4}/\avg{\sigma^2}^2\)。由 \(\operatorname{Var}(\sigma^2)=\avg{\sigma^4}-\avg{\sigma^2}^2\ge0\) 得 \(\kappa\ge3\)，等号仅在 \(\sigma\) 不随机时成立。\(\mu\) 子散射中 \(\sigma\propto1/(\beta p)\)，所以 \(\kappa=3\avg{(\beta p)^{-4}}/\avg{(\beta p)^{-2}}^2\)。
加权最小二乘
数据 \((x_i,y_i)\)，\(y_i\) 的标准差为 \(\sigma_i\)，模型 \(y=A+Bx\)。在高斯误差下，使 \(\chi^2=\sum w_i(y_i-A-Bx_i)^2\)（\(w_i=1/\sigma_i^2\)）最小的 \(A\)、\(B\) 是最可能的参数。令 \(\partial\chi^2/\partial A=\partial\chi^2/\partial B=0\)： \[A[1]+B[x]=[y],\qquad A[x]+B[x^2]=[xy],\] 其中 \([f]=\sum w_if_i\)。解得 \begin{align*} B&=\frac{[1][xy]-[x][y]}{\Delta}, \\ A&=\frac{[x^2][y]-[x][xy]}{\Delta}, \\ \Delta&=[1][x^2]-[x]^2 . \tag{C.1} \end{align*} \(B\) 是 \(y_i\) 的线性组合 \(B=\sum c_iy_i\)，\(c_i=w_i([1]x_i-[x])/\Delta\)；各 \(y_i\) 独立， \begin{align*} \sigma_B^2&=\sum c_i^2\sigma_i^2=\sum w_i([1]x_i-[x])^2/\Delta^2\\&=[1]\big([1][x^2]-[x]^2\big)/\Delta^2=[1]/\Delta. \end{align*} 同理 \(\sigma_A^2=[x^2]/\Delta\)。若模型正确且误差估计正确，\(\chi^2_{\min}\) 的期望值等于自由度 \(\nu\)（数据点数减参数个数）。
附录 D 几何：立体角与接收度
用投影倾角表示立体角
设方向单位矢量 \(\uvec n=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)\)，定义 \(t_x=\tan\theta_x=n_x/n_z\)，\(t_y=\tan\theta_y=n_y/n_z\)。射线与平面 \(z=1\) 交于点 \((t_x,t_y,1)\)，距原点 \(r=\sqrt{1+t_x^2+t_y^2}\)。平面上的面元 \(\dd t_x\dd t_y\) 对原点所张立体角等于面元在垂直于视线方向上的投影 \(\dd t_x\dd t_y\cos\chi\) 除以 \(r^2\)，其中 \(\chi\) 是平面法线（\(z\) 轴）与视线的夹角，\(\cos\chi=1/r\)。所以 \[\tag{D.1} \dd\Omega=\frac{\dd t_x\,\dd t_y}{(1+t_x^2+t_y^2)^{3/2}} .\] 检验：对整个平面积分， \begin{align*} &\int\!\!\int(1+t_x^2+t_y^2)^{-3/2}\dd t_x\dd t_y\\&=2\pi\int_0^\infty s(1+s^2)^{-3/2}\dd s=2\pi, \end{align*} 正好是上半球的立体角。
两个角积分
\(\int_0^{2\pi}\abs{\cos\varphi}\dd\varphi=4\int_0^{\pi/2}\cos\varphi\,\dd\varphi=4\)；\(\int_0^{2\pi}\abs{\sin\varphi\cos\varphi}\dd\varphi=4\int_0^{\pi/2}\tfrac12\sin2\varphi\,\dd\varphi=2\)。展开 \begin{align*} &(d-Dt\abs{\cos\varphi})(d-Dt\abs{\sin\varphi})\\&=d^2-dDt(\abs{\cos\varphi}+\abs{\sin\varphi})\\&\quad+D^2t^2\abs{\sin\varphi\cos\varphi}, \end{align*} 逐项积分得 \(2\pi d^2-8dDt+2D^2t^2\)。检验：\(D\to0\) 时回到 \(2\pi d^2\)（两平面重合，接收度就是平面面积乘立体角）。
附录 E 名词解释
| 术语 | 英文 | 一句话定义 |
|---|---|---|
| \(\mu\) 子 | muon | 带电轻子，质量 105.7 MeV/c2，平均寿命 2.197 μs |
| 初级/次级宇宙线 | primary/secondary cosmic rays | 来自太空的粒子/它们在大气中产生的粒子 |
| 空气簇射 | air shower | 一个高能初级粒子在大气中引发的粒子级联 |
| 固有时 | proper time | 随粒子运动的钟所走过的时间 |
| 衰变长度 | decay length | \(\beta\gamma c\tau\)，粒子衰变前的平均飞行距离 |
| 不透明度 | opacity | 沿路径的面密度 \(\int\rho\,\dd l\) |
| 米水当量 | meter water equivalent | 不透明度单位，1 m.w.e. \(={100\,}\mathrm{g/cm^2}\) |
| 阻止本领 | stopping power | 单位路程（或面密度）的平均能量损失 |
| 最小电离粒子 | minimum ionizing particle | 能损处于 Bethe 曲线极小值附近的粒子 |
| 平均激发能 | mean excitation energy | Bethe 公式中刻画原子电子束缚的能量 \(\bar I\) |
| 密度效应 | density effect | 介质极化屏蔽高速粒子电场、抑制能损增长 |
| 临界能量（\(\mu\) 子） | muon critical energy | 电离与辐射损失相等的能量 \(a/b\) |
| 连续减速近似 | CSDA | 把能损看作沿路径连续、确定的过程 |
| 辐射长度 | radiation length | 高能电子因轫致辐射能量降为 \(1/\ee\) 的面密度 |
| 多次库仑散射 | multiple Coulomb scattering | 大量小角库仑散射累积的随机偏折 |
| 强度 | intensity | 单位时间、单位垂直面积、单位立体角的粒子数 |
| 接收度 | acceptance | 探测器对某方向范围的“有效面积 × 立体角” |
| 显著性 | significance | 观测偏离无信号预言的标准差个数 |
| 最近点法 | point of closest approach (PoCA) | 用入射、出射径迹的最近点估计散射位置 |
| \(\mu\) 子成像 | muography | 用宇宙线 \(\mu\) 子的透射或散射为物体内部成像 |