量子密钥分发（QKD）听起来离日常的数据分析工作很远，但它的数学骨架其实非常朴素：不过是一些复向量、酉矩阵和概率运算，而这些恰恰是R语言最擅长的事情。以E91协议为代表的纠缠型量子密钥分发，要求通信双方Alice和Bob各自对纠缠粒子做随机测量，再通过经典信道比对部分结果来检测窃听者。整个过程完全可以在R里做成一个蒙特卡洛仿真，直观地看到窃听行为如何抬高误码率、隐私放大如何压缩密钥。本文就把这套流程完整地实现一遍。
量子态的R语言表示方法
单个量子比特的状态可以用一个二维复向量描述。在R中，我们直接用c()构造列向量，用%*%做矩阵乘法。对于纠缠态，最重要的是Bell态之一，也就是通常写作Phi+的那个态。它是一个四维向量，含义是两个粒子要么同时处于0态，要么同时处于1态，且概率各占一半。测量结果的关联性正是密钥协商安全性的来源。
测量算符用投影矩阵表示。E91协议里双方使用三组测量基，彼此夹角为45度和22.5度。由于R原生支持复数运算，构造含复分量的测量基毫无障碍。下面这段代码定义了基本工具函数：
library(MASS)
# 单比特基矢
ket0 <- matrix(c(1, 0), ncol = 1)
ket1 <- matrix(c(0, 1), ncol = 1)
# 构造双粒子直积态
tensor <- function(a, b) {
kronecker(a, b)
}
# Bell态 |Phi+> = (|00> + |11>) / sqrt(2)
bell_phi_plus <- (tensor(ket0, ket0) + tensor(ket1, ket1)) / sqrt(2)
# 以角度theta构造测量基（二维实向量的旋转即可覆盖E91需求）
basis <- function(theta) {
v1 <- matrix(c(cos(theta), sin(theta)), ncol = 1)
v2 <- matrix(c(-sin(theta), cos(theta)), ncol = 1)
list(v1, v2)
}
有了这些积木，任意两体测量结果的联合概率都可以通过投影算符的期望值算出来：Re(t(Conj(v)) %*% P %*% bell_phi_plus)的形式，其中P是对应投影矩阵的张量积。这套写法的好处是完全向量化，不需要写任何显式循环就能扩展到更多协议变体。
E91协议的完整仿真流程
E91的核心流程是：纠缠源产生大量Bell对，分别发送给Alice和Bob；双方各自独立地随机选择一组测量基进行测量；随后通过公开信道公布基矢选择，保留双方使用匹配角度的回合作为原始密钥，再随机抽样部分回合估计误码率，判断是否存在窃听。
窃听者Eve的行为可以这样模拟：她以概率p对粒子进行截获并用自己的基测量，再把粒子转发出去。一旦她测量的基与合法方不同，纠缠关联就被破坏，误码率随之上升。下面的函数实现了单轮仿真：
# 计算在给定双方基角(theta_a, theta_b)下
# Alice与Bob测量结果一致的概率
prob_match <- function(theta_a, theta_b) {
ba <- basis(theta_a); bb <- basis(theta_b)
p <- 0
for (i in 1:2) for (j in 1:2) {
proj <- kronecker(ba[[i]], bb[[j]]) %*%
t(kronecker(ba[[i]], bb[[j]]))
amp <- t(Conj(kronecker(ba[[i]], bb[[j]]))) %*% bell_phi_plus
# 同向偏振记为结果一致：1-1 或 2-2
if (i == j) p <- p + Mod(amp[1, 1])^2
}
p
}
# 单轮E91仿真：n对粒子，窃听概率p_eve，基噪sigma
simulate_e91 <- function(n = 1000, p_eve = 0, sigma = 0) {
angles <- c(0, pi/8, pi/4) # 双方共用的三组基角
results <- data.frame(alice = integer(n), bob = integer(n),
match = logical(n))
for (k in 1:n) {
ia <- sample(1:3, 1); ib <- sample(1:3, 1)
pm <- prob_match(angles[ia], angles[ib])
# 窃听破坏关联：窃听回合一致概率退化为0.5附近
if (runif(1) < p_eve) pm <- pm * 0.5 + 0.25
# 信道噪声轻微扰动
pm <- pmin(1, pmax(0, pm + rnorm(1, 0, sigma)))
same <- rbinom(1, 1, pm)
results$alice[k] <- sample(0:1, 1)
results$bob[k] <- ifelse(same == 1, results$alice[k],
1 - results$alice[k])
results$match[k] <- (ia == ib)
}
results
}
运行仿真后，只保留match为真的回合，这些回合里Alice和Bob的比特在无窃听、无噪声时应当完全一致。抽样比对一段子集，如果误码率超过协议阈值（通常取11%左右），双方直接放弃本轮密钥；否则进入密钥提炼阶段。这个判定逻辑完全对应真实E91设备中的后处理流程。
窃听强度与密钥质量的数值分析
仿真的价值在于可以做参数扫描。我们把窃听概率从0逐步拉到1，每个点重复多轮取平均，观察保留回合的误码率变化。理论上，全功率窃听会把一致概率推向0.75，对应约25%的误码率，远超安全阈值，因此会被协议自动检出。
sweep <- sapply(seq(0, 1, by = 0.1), function(p) {
n <- 4000
keep <- simulate_e91(n, p_eve = p, sigma = 0.02)
keep <- keep[keep$match, ]
mean(keep$alice != keep$bob) # 误码率QBER
})
plot(seq(0, 1, by = 0.1), sweep, type = "b", pch = 19,
xlab = "窃听概率", ylab = "误码率 QBER",
main = "E91协议窃听检测灵敏度")
abline(h = 0.11, lty = 2) # 安全阈值
从曲线可以看出两个关键现象。第一，误码率随窃听概率近似线性上升，这解释了为什么QKD的安全性检测是定量的而非定性的——窃听者哪怕只碰十分之一的粒子，也会留下可统计的痕迹。第二，信道本身的噪声会抬高基线误码率，因此工程上必须在窃听阈值与设备损耗之间留出余量，这也是现实量子网络要部署诱骗态和纠错编码的原因。
最后一步是隐私放大。在R里最简单的做法是把保留比特按块做异或或哈希压缩，压缩比例根据观测到的误码率由安全码率公式决定：r约等于1减去两倍的QBER。写一个循环按块异或，几十行代码就能完成，得到的最终密钥在信息论意义上对Eve不可知。把这套仿真跑通之后，你不仅得到了一份可复现的R代码，更重要的是建立了对纠缠、测量关联和安全证明之间的直觉，这种直觉读多少论文都不如亲手算一遍来得扎实。