引言¶
在前一讲中,我们用序参量区分有序与无序,用临界指数描述相变附近的变化,并通过标度律把这些指数联系起来。但还有一个问题没有解决。这些行为能否从一个具体的理论中算出来?
例如,知道序参量满足 \(m\sim(-t)^\beta\),只是知道它怎样被描述。要解释这样的幂律,我们需要知道系统为什么从偏好 \(m=0\) 转向选择 \(m\ne0\),以及为什么在接近转变时对扰动更加敏感。
本讲的两个核心概念,分别回答两个相接的问题。
Landau 理论用序参量的自由能描述不同状态的竞争,通过自由能极小值预测系统选择什么状态。 在最简单的情形下,自由能是一个四次多项式,求极小值便能得到连续相变及一组平均场临界指数。
Ginzburg 判据通过比较相关区域内的涨落与平均场预测的有序程度,检查“只围绕一个平均状态计算”是否可靠。 它不是另一套相变模型,而是对前一种近似的自洽性检验。
连接两者的关键,是把空间重新纳入理论。不同位置并不总是处于相同状态,它们之间还存在关联。本讲据此逐步推进,从确定变量一直走到检查近似。
确定序参量 → 写出自由能 → 求平均场状态 → 计入空间涨落 → 检查平均场
为了让这些概念对应到可见的现象,我们从一篇关于细菌涡旋晶格的实验论文出发。Wioland 等人在一篇发表于 Nature Physics 的论文中,让大量游动的细菌在彼此连通的微小圆形腔室内形成涡旋,研究这些局部旋转怎样组织成整体的有序排列[1]。他们发现,改变腔室之间的通道宽度,就能改变相邻涡旋对同向或反向旋转的偏好,并建立包含相互作用与涨落的模型解释这一现象,再通过平均场近似,将其简化为与本讲Landau 理论密切相关的四次势模型。
1. 细菌涡旋中的有序与涨落¶
这项细菌涡旋实验的意义,在于它把局部运动如何形成集体有序这一普遍问题变成了可以直接观察和调节的现象。单个腔室形成涡旋,并不能决定整个阵列怎样排列[1]。这说明,整体的有序不能只从孤立单元的性质判断,还取决于单元之间怎样相互影响。同时,有序区域中的旋转仍在不断波动,因此我们不仅要解释系统为什么偏好某种状态,还要判断用平均状态代替实际运动时忽略了什么,以及这些涨落是否足以影响预测的可靠性。实验本身并不是平衡系统,但它提出的这两个问题,将分别引出本讲平衡参照模型中的 Landau 理论与 Ginzburg 判据。
1.1 同样的细菌,为什么会形成不同的旋转排列?¶
局部旋转的形成与整个阵列的有序排列,是两个不同的问题。这篇论文的实验保持细菌种类和腔室基本结构不变,通过改变腔室之间的间隙,比较邻近涡旋的旋转关系。
Wioland 等人在微流控装置中排列了一组浅圆形腔室,并用狭窄通道连接相邻腔室。在适当的限制条件下,许多细菌的运动形成一个整体旋转的涡旋;腔室之间的流动则使这些涡旋相互影响[1]。实验首先提出了一个很具体的问题。一个腔室中的涡旋,会偏好与相邻涡旋同向旋转,还是反向旋转?
腔室间隙为 \(6\,\mu\mathrm m\)。视频前半部分左侧为明场影像,右侧叠加旋转变量的伪彩色,其中绿色与紫色对应相反旋转方向,颜色不是细菌种类或密度的标记;后半部分展示局部流场。应注意局部的交替排列,而不是寻找覆盖整个装置的完美棋盘。视频源自: Wioland, H., Woodhouse, F. G., Dunkel, J., & Goldstein, R. E. (2016). Ferromagnetic and antiferromagnetic order in bacterial vortex lattices. Nature physics, 12(4), 341-345.
现在把间隙加宽,排列方式发生变化,局部相邻腔室更倾向出现同向旋转。
间隙放宽,出现的是偏好同向旋转的区域,不是所有涡旋永远保持同一方向。视频源自: Wioland, H., Woodhouse, F. G., Dunkel, J., & Goldstein, R. E. (2016). Ferromagnetic and antiferromagnetic order in bacterial vortex lattices. Nature physics, 12(4), 341-345.
论文把同向旋转排列称为“铁磁序”,把交替旋转排列称为“反铁磁序”。这些名称描述的是排列关系,并不意味着细菌具有这里所讨论的磁性。 两种排列都可能有序;真正变化的是邻居之间偏好的关系。
有序并不等于静止,也不等于没有涨落。 局部旋转强度会变化,有序区域会分裂、合并或翻转。我们既要描述稳定的统计倾向,也要描述偏离这种倾向的变化。后者正是平均场近似需要接受检验的原因。
1.2 从一幅流场到一个序参量¶
视频让我们看见了不同的排列,但理论还需要能够比较和计算的变量。由于当前关心的是腔室的整体旋转,而不是每个细菌的运动细节,可以先把一个腔室内的复杂流场压缩为一个有正负的实数,再检查怎样的整体平均能够识别不同的有序。
用 \(V_i\) 表示第 \(i\) 个腔室的旋转变量,正负号区分两种旋转方向,绝对值表示旋转强弱。
原论文通过粒子图像测速得到流场,再根据腔室内的角动量式空间平均定义 \(V_i\);流速还按同次实验的均方根速度归一化[1]。因此,\(V_i\) 不是某一个细菌的速度,也不是简单的像素颜色平均。这里首先要把握的物理含义是,一个数描述一个腔室的集体旋转。
如果研究同向旋转,可以再对 \(N\) 个腔室求平均,定义
若许多腔室偏好同一方向,它们的贡献能够累积,\(m\) 就可能明显偏离零;若正负方向相互抵消,\(m\) 则接近零。于是,一个庞大的运动系统开始拥有了一个可计算的整体描述。
这里的平均只是定义一个整体观测量,并没有假定各个腔室的状态相同。后面用一个均匀值代表不同位置的状态来求解,才引入本讲所说的平均场近似。
但这里必须停下来检查变量是否合适。一个完美的交替旋转排列,也可能满足 \(m=0\),却显然不是无序。对正方晶格,给两种交替位置分别加上正负权重,可以定义
其中 \(i_x,i_y\) 是腔室的整数格点坐标。原先相邻涡旋的正负交替被这个权重抵消,它们便能在 \(m_{\mathrm{alt}}\) 中同号累积。这就是针对交替排列选择的序参量。
这个例子具体说明了第 4 讲中“序参量必须对应所研究的有序”的含义。不能先随意选一个平均值,再把它等于零解释为没有结构。
同向排列与交替排列需要不同的序参量,但一旦选定所研究的有序,就可以提出一个共同的问题——这种有序怎样从无序背景中建立起来? 接下来,我们先固定一种有序模式,用它对应的序参量 \(m\) 描述有序程度,研究系统为什么从偏好 \(m=0\) 转向偏好 \(m\ne0\)。理解这个基本过程之后,再回到论文讨论不同排列之间的竞争。
1.3 从描述有序到解释有序¶
序参量告诉我们系统处于怎样的状态,却没有解释它为什么选择这种状态。要从描述走向预测,还需要比较不同有序程度的稳定性。对于平衡系统,自由能正是完成这一比较的工具。
细菌依靠持续消耗能量维持运动,因此我们先在平衡参照模型中建立这套方法,再回到实验。实验启发的有正负的集体变量、邻近区域之间的相互作用,以及围绕有序状态的涨落,将作为贯穿两者的线索。
这个参照模型考虑有限温度下、短程相互作用引起的普通连续相变,序参量只有一个实数分量,而且在没有外场时,\(m\) 与 \(-m\) 等价。我们将看到,对称性与稳定性如何引出一个关于序参量的四次多项式,也就是最简单的标量四次理论,常写作 \(\phi^4\) 理论。
接下来第 2—5 节依次完成建立自由能、求出平均场状态、检验涨落。第 6 节再回到细菌涡旋,理解这套描述怎样进入论文的有效模型。
1.4 平均值相同,排列也相同吗?¶
序参量的选择会决定哪些结构被保留下来。把同向、交替和无序的旋转排列放在一起,可以直观看到,普通平均值能够识别整体旋转,却可能把另外两种截然不同的排列混在一起。理解这种差别,才知道后面自由能中的 \(m\) 究竟在描述什么。
可以先对图中的排列做一次具体计算。设每个腔室的旋转强度均为 \(v>0\),蓝色对应 \(+v\),橙色对应 \(-v\)。左图的 \(16\) 个腔室都是蓝色,因而 \(m=v\)。中图与右图则各有 \(8\) 个蓝色腔室和 \(8\) 个橙色腔室,在这个等强旋转的理想示例中,两者都给出
相同的结果来自平均运算本身。它累加每个腔室的旋转,却不记录这些腔室位于哪里。即使把中图的蓝色与橙色重新排列,破坏原有的交替规则,只要两种颜色的数量不变,\(m\) 就不会改变。普通平均值保留了整体旋转的偏向,却丢掉了旋转方向怎样分布的位置信息。 因此,\(m=0\) 只能说明正负贡献抵消,不能独自回答排列是否有序。
交替序参量补入的正是这种位置信息。把棋盘的两类位置分别赋予 \(+1\) 和 \(-1\) 的权重,并让正权重对应中图的蓝色位置,原来的 \(+v\) 保持不变,原来的 \(-v\) 则乘以 \(-1\),也变成 \(+v\)。于是,原本在普通平均中抵消的贡献,在交替平均中全部累积,得到 \(m_{\mathrm{alt}}=v\)。这个计算说明,序参量不只是求平均,还规定了我们要识别哪一种排列规则。
选定有序模式后,还要区分序参量的符号与大小。如果把左图所有涡旋同时反转,\(m\) 会从 \(v\) 变成 \(-v\),但所有腔室仍然同向旋转,并没有变得更无序。对交替排列整体反转,也只是交换棋盘上的两种颜色,使 \(m_{\mathrm{alt}}\) 改变符号。正负号区分同一种有序的两个方向,绝对值才衡量这一选定模式的有序程度。 这就是后面比较 \(m\) 与 \(-m\) 时所指的两种状态。
由此也能明确 Landau 理论接下来要解决的问题。它不是用一个数囊括所有可能的排列,而是先选定一种有序及其序参量,再研究这种有序能否建立、建立后有多强。在前面约定的平衡参照模型中,若没有偏好某个方向的外场,正负两个方向应具有相同的自由能;但系统究竟偏好零序参量,还是某个非零值,仍需要计算。下一节的自由能 \(f(m)\),就是为比较这些不同的有序程度而引入的。
2. Landau 理论如何描述状态的竞争¶
序参量能够区分状态,却还不能解释系统为什么选择其中某一种。Landau 理论用自由能补上这一环节,把不同有序程度之间的竞争转化为一个可以比较和求解的问题。我们先理解这种有效描述的思路,再用对称性与稳定性约束自由能,最后说明寻找极小值时采用了怎样的近似。
2.1 为什么不从每一个粒子的运动算起?¶
直接求解全部微观运动,未必是理解相变最有效的起点。如果不同系统的有序可以由相同类型的变量描述,就有理由先寻找它们共有的结构,再把材料差异放入有效参数。Landau 的相变理论正是沿着这一思路建立的,它使建模不再完全依赖于微观问题是否能够精确求解。
列夫·达维多维奇·朗道(Lev Davidovich Landau,1908—1968)是苏联理论物理学家。1937 年,他在相变理论工作中把序参量和对称性置于核心位置[2]。这一方法的重要性,不在于用多项式替代了一条复杂曲线,而在于改变了提问的顺序,先识别相变前后改变的有序,再为这个有序建立最简的有效描述。
Landau 后来与 Lifshitz 合著的《统计物理》系统介绍了这种方法[3]。1950 年,Ginzburg 与 Landau 又把能够随空间变化的序参量用于超导理论[4]。本讲使用的实标量模型并不是那套超导方程的全部内容,但共享同一种研究策略。不必先解出全部微观细节,也可以用序参量、对称性和稳定性约束理论。
以刚才的旋转方向为直观参照,微观问题可能包含每个细菌的位置、朝向和周围流动;宏观描述则可以先考察系统是否偏好某种整体旋转。当我们只研究这个问题时,许多微观细节不必逐一出现,而是进入有效描述的系数中。
这种压缩并不保证答案正确。它只给出了一个能够计算的起点。Landau 理论的力量与局限,必须在同一条推理链中理解。
2.2 多项式的每一项分别意味着什么?¶
第 3 讲说明,在给定温度等热力学条件下,平衡状态由合适的自由能决定。现在需要把这一原则落实到序参量上,用一个函数比较不同有序程度的代价。对称性决定这个函数允许出现哪些项,稳定性则决定哪些项不能省略;四次多项式将从这些条件中产生。
设系统被约束在某个均匀序参量 \(m\),用 \(f(m)\) 表示这个约束状态的自由能密度。它回答的不是“系统此刻在哪里”,而是“维持不同的有序程度,分别具有怎样的自由能代价”。
对于连续相变,序参量在转变附近从零连续出现。这时,如果粗粒化后的局域描述可以作解析展开,便可以按 \(m\) 的幂次组织它。零外场下两种方向等价,要求
因此,\(m\)、\(m^3\) 等奇次项不能单独出现。保留能够描述普通连续相变的最低阶项,得到
这就是本讲的 Landau 自由能,其中取 \(u>0\)。各项在状态竞争中承担着不同的作用。
- \(f_0\) 与 \(m\) 无关,是规则背景;它不改变极小值的位置。
- \(r\) 决定原点附近的曲率。\(r>0\) 时,较小的非零 \(m\) 增加自由能;\(r<0\) 时,离开零点反而降低自由能。
- \(u>0\) 提供最低阶稳定作用。当 \(r<0\) 鼓励非零有序时,四次项阻止 \(m\) 无限制增大。
- \(h\) 是与 \(m\) 共轭的外场。它通过 \(-hm\) 偏好一种方向,使原本对称的两个状态不再等价。
这里的“外场”不必是磁场,而是指自由能中与序参量线性耦合的控制量。对旋转变量而言,可以把它理解为选择一种旋转手性的偏置;这不表示论文已经施加或标定了这样的实验控制。
读图时最重要的不是记住三条曲线,而是看清二次项与四次项之间的竞争。二次项决定是否倾向建立有序,四次项决定这种增长在哪里停止。 当二次项改变符号,最低点便可以从原点连续移向两侧。
这里展开的是保留序参量后的有效局域自由能,不是已经完成所有长尺度涨落统计的精确热力学自由能。因此,即使以一个光滑的局域多项式为起点,最终的平衡自由能仍然可以在相变处出现非解析行为。
2.3 求极小值,就是怎样的近似?¶
写出自由能之后,还需要说明它怎样给出平衡状态。寻找极小值并不是另加的一条规则,而是来自不同状态之间的统计权重比较。关键在于,我们先只比较空间均匀的构型;这种限制使问题容易求解,也留下了之后必须检查的空间涨落。
在平衡参照模型中,一个均匀状态的统计权重与
成正比。这里 \(\mathcal V\) 是系统体积,\(T\) 是绝对温度,\(k_B\) 是 Boltzmann 常数;不用 \(V\) 表示体积,是为了与腔室的旋转变量 \(V_i\) 区分。当体积很大时,很小的自由能密度差也会产生很大的权重差。于是,在只允许均匀状态的计算中,最低点主导平衡答案。
先令 \(h=0\)。驻点条件为
它给出 \(m=0\),以及在 \(r<0\) 时才存在的两个非零解。再通过二阶导数判断这些解的稳定性,有
当 \(r>0\) 时,\(m=0\) 的曲率为正,是稳定极小值。当 \(r<0\) 时,原点曲率为负,变成不稳定的极大值;新的稳定解为
两种符号代表两个等价的有序分支。自由能没有偏好哪一个方向,但一个宏观有序相可以选择其中之一,这就是自发对称性破缺。在有限系统中,若观察时间足够长,两种方向之间的翻转可能使总平均重新为零;讨论非零自发序参量时,指的是先选定一个有序分支,或者按热力学极限后的零场极限定义它。
这里的“平均场”近似已经可以明确辨认。我们用一个均匀值 \(m_0\) 描述平衡背景,并用它的极小值方程求状态,尚未计算空间涨落对这个方程的反馈。 系统实际是否处处接近这个背景,仍然没有得到检验。
因此,“Landau 理论”与“平均场近似”相关,却不完全同义。前者提供有效自由能的形式;后者是求解它的一种近似。即使保留同样的四次理论,认真处理涨落后,也可能得到不同于均匀极小值的临界行为。
3. 从平均场预测到临界行为¶
有了稳定极小值,便可以把自由能的几何结构转化为可观测的变化。最低点的位置决定序参量,附近的曲率决定对外场的响应,而沿平衡分支变化的自由能决定熵和比热。把这些结果与前一讲的定义对照,就能看见临界指数如何从同一个理论中被算出来。
3.1 序参量为什么按平方根出现?¶
最低点从原点移向两侧,说明有序开始建立,但还没有说明它出现得有多快。要得到临界指数,需要把自由能系数的变化与临界距离联系起来,再考察稳定解的领先幂次。这一步将把图中的连续分岔转化为序参量的平方根规律。
前一讲已经引入约化温度,用来表示与临界点的相对距离。现在先在均匀平均场理论内部,让 \(r\) 在温度 \(T_0\) 处变号,写成
\(T_0\) 是这一近似预测的临界温度;保留这个下标,是因为涨落可能把真实临界温度移到另一个位置。\(a\) 表示二次系数在临界附近变化的斜率,\(u\) 的平滑变化则暂取为临界附近的正值。
代入刚才的极小值,在选定的正向有序分支上有
与定义 \(m_0\sim(-t_0)^\beta\) 比较,得到
这个指数并不是另行假设的。它来自极小值方程中 \(|r|m\) 与 \(um^3\) 的平衡。两者相当,意味着 \(m^2\sim |r|/u\)。临界指数第一次从自由能的结构中被算了出来。
这里计算的是理想化平衡模型的集体序参量,不是说单个腔室出现稳定旋转,就已经发生了一个热力学连续相变。单个变量的双稳态和无限多自由度的集体相变,需要始终分开。
3.2 势阱变平,为什么响应会变大?¶
只知道最低点的位置还不够。同样稳定的两个状态,对扰动的敏感程度也可能不同。为描述这种差异,我们考察一个微小外场怎样移动最低点,并把位移的大小与势阱曲率联系起来。响应的增强将由此成为可以计算的结果。
恢复一个很小的外场 \(h\),极小值条件变成
这就是把外部偏置与有序程度联系起来的状态方程。定义线性响应率 \(\chi=\partial m/\partial h\),在固定温度下对状态方程求导,得到
在零场稳定分支上,令
便得到
\(\kappa\) 是最低点附近的曲率,也就是偏离这个状态时的恢复强度。曲率大,偏离最低点需要较大代价;曲率小,同样的外场便能把平衡值推移得更远。响应发散在自由能图像中的含义,就是稳定点附近越来越平。
无序侧 \(m_0=0\),所以 \(\kappa_+=r\);有序侧 \(m_0^2=-r/u\),所以 \(\kappa_-=2|r|\)。因此
两侧都给出 \(\gamma_{\mathrm{MF}}=1\),但振幅不同。到 \(r=0\) 时,状态方程只剩 \(h=um^3\),于是又得到临界等温线指数 \(\delta_{\mathrm{MF}}=3\)。三个指数来自同一个极小值方程,而不是三套互不相关的经验规律。
3.3 没有潜热,为什么仍然是相变?¶
序参量与响应揭示了有序状态的变化,热力学上的相变还需要从平衡自由能中辨认。本讲的普通连续相变没有有限的熵跳变,但自由能的更高阶导数仍可能不连续。把稳定解代回原来的多项式,可以看清没有潜热时,比热异常怎样出现。
把稳定解代回自由能,扣除规则背景后得到
无序侧的最低点仍在原点;有序侧选择非零 \(m_0\),获得一个有限的自由能降低。两条分支在 \(t_0=0\) 处具有相同的函数值和一阶温度导数,但二阶导数不同。
这里可以用第 3 讲的热力学关系理解第 4 讲的相变分类。熵密度为 \(s=-\partial f_{\mathrm{eq}}/\partial T\),比热密度为 \(C=-T\partial^2f_{\mathrm{eq}}/\partial T^2\)。这里熵没有有限跳变,因此没有潜热;比热的有序相关部分却出现有限跳变,对应 \(\alpha_{\mathrm{MF}}=0\)。这不是说全部比热为零,也不是说 \(\alpha=0\) 一定代表对数发散;具体形式必须由计算决定。
潜热可以在这里作一个简短对照。恒压可逆相变中,同一份物质的两相在共存点满足 \(\Delta G=0\),因而有
其中 \(G\) 是 Gibbs 自由能,\(T_{\mathrm{tr}}\) 是相变温度,\(L\) 是吸收的潜热。缓慢加热标准大气压下的纯冰水混合物,冰未融尽时,输入的热量可以改变两相比例而不使温度明显上升。加热不必等于升温,因为能量还可以用于改变物质的组织方式。 本节的连续相变没有这种有限的熵跳变,奇异性却仍可出现在更高阶导数中。
判断相变不能只看有没有吸热平台。一阶相变要求适当自由能的某个一阶导数发生跳变,但这个导数不一定是熵,因此也不一定伴随潜热。这里的关键是,本讲算出的自由能与一阶温度导数都连续,异常从更高阶导数中显现。
至此,平均场给出的答案已经相当完整。但仔细检查推导,会发现它几乎没有用到空间维度。这不是维度不重要,而是我们还没有允许不同位置表现得不一样。 接下来必须补上这一点。
4. 空间涨落与相关长度¶
前面的计算把整个系统压缩为一个均匀序参量,因此能够描述整体偏好,却看不见不同区域之间的差异。现在需要恢复空间依赖,让局部偏离及其相互影响进入自由能。邻近区域的协调代价与局域恢复强度共同决定相关长度,也为估计涨落的大小提供依据。
4.1 邻近区域之间,还缺少一个协调代价¶
回到视频,整个阵列并不总处于一种统一旋转状态。一片区域可能偏向一个方向,另一片区域则可能不同。要保留这些结构,既需要让序参量随位置变化,也需要描述邻近区域之间的相互作用;仅仅把同一个局域势写在每个位置,并不能解释空间协调。
因此,我们把均匀序参量推广为 \(m(\mathbf x)\),其中 \(\mathbf x\) 表示位置,\(m(\mathbf x)\) 表示附近一小块区域的有序程度。这正是第 1 讲中粗粒化的基本想法——减少微观细节,但保留仍然重要的空间结构。粗粒化并不等于把整个系统替换成一个常数。
仅仅给每个位置写一个相同的局域势,还不够。如果邻近区域可以毫无代价地独立选择正负方向,就没有理由形成延伸的同向有序区域。对于偏好相邻变量相近的相互作用,需要给空间变化增加代价。
在离散格点上,这一点可以直接从一条耦合键看出。取 \(J>0\),则
最后两项可以并入局域二次系数;第一项则明确惩罚相邻位置的差异。当变化足够缓慢,差分可以换成空间导数,得到 Landau–Ginzburg 自由能泛函
“泛函”意味着输入不再是一个数,而是整个空间分布,输出仍是一个自由能。积分把各处的贡献加起来,\(f(m)\) 就是第 2 节的局域四次自由能。\(c>0\) 描述空间变化的代价,\(d\) 是空间维度,\(\nabla m\) 描述场在空间中改变得有多快。与 \(m\) 无关的背景项不影响下面的变分和涨落,可以省略。
这两个部分的分工由此变得清楚。局域势决定每一处偏好怎样的取值,梯度项决定相邻位置怎样协调。 对交替排列,适合平滑描述的是先消去正负交替的序参量场,而不是强行把原始交替变量当作缓变场。
4.2 相关长度来自两种代价的竞争¶
梯度项使局部变化能够影响邻近区域,但这种影响不会在所有距离上同样强。局域势倾向于把每一处拉回稳定值,空间耦合则倾向于让相邻位置共同变化。比较这两种作用,可以得到影响衰减的特征长度,并解释它为什么在临界点附近增长。
选择一个稳定背景 \(m_0\),把实际场写成
\(\delta m\) 表示偏离背景的变化。由于 \(m_0\) 已经满足极小值方程,线性项消失;只保留二次项,得到
这里出现的仍是上一节的曲率 \(\kappa\)。第一项反对场离开局部稳定值,第二项反对这种偏离在空间中变化得太快。只保留二次项称为高斯近似,因为相应的涨落权重是高斯形式。
设某个偏离在长度 \(L\) 上显著改变。它的梯度量级约为 \(\delta m/L\),于是两项的相对重要性由 \(\kappa\) 与 \(c/L^2\) 决定。它们相当的长度为
这就是当前高斯理论中的相关长度。它不是人为画出的团块直径,而是局域恢复与空间协调竞争产生的特征尺度。这个长度也会出现在小扰动的平衡方程中。远离扰动源的一维方向上,有
衰减解为 \(\delta m\propto e^{-x/\xi}\)。在这个平衡高斯模型中,涨落—响应关系把静态响应与自发涨落的连通关联联系起来,二者具有相同的衰减长度。因此,这个响应长度也就是相关长度。更高维度的关联函数还带有几何因子,但离开临界点时,长距离衰减仍由这个 \(\xi\) 控制。
因此,势阱变平既使一个很小的外场能够造成很大响应,也使局部偏离能够影响更远的位置。在当前近似下,二者由同一个曲率联系起来,满足
这些是本模型的高斯结果,不是对所有系统的无条件恒等式。由于平均场曲率满足 \(\kappa\propto|t_0|\),所以 \(\xi\propto|t_0|^{-1/2}\),得到 \(\nu_{\mathrm{MF}}=1/2\)。前一讲中独立定义的响应与相关长度,现在有了共同的物理来源。
4.3 为什么有了高斯涨落,还需要继续检查?¶
相关长度说明涨落能够延伸多远,却还没有告诉我们涨落的幅度是否足够小。高斯计算保留二次项,预先要求更高阶项只带来小修正。为检查这个要求,需要估计一个相关区域的平均值会偏离背景多少,再把这种偏离与有序背景作比较。
考虑一个边长约为 \(\xi\) 的区域,并把其中的场作空间平均,记作 \(\overline m_\xi\)。在选定的有序相中,它的统计均值为 \(m_0\),但不同构型下的区域平均值仍会变化。用
表示这种变化的方差。尖括号表示对许多可能构型取统计平均,不是再次对同一张图作空间平均。
为什么选 \(\xi\),而不是整个样品的大小?很大的样品可以包含许多相关区域,即使每个区域内部涨落明显,总平均仍可能稳定。宏观平均很稳定,不等于区域内部可以忽略关联。 检查临界理论,需要关注决定临界行为的相关尺度。
在尺度 \(\xi\) 上,局域二次项与梯度项量级相同。因此,一个相关区域的平均值偏离背景 \(\delta\overline m_\xi\) 时,其自由能代价在量级上为
其中 \(\xi^d\) 是相关体积。把这一二次代价放入平衡权重 \(e^{-\Delta F_\xi/(k_BT)}\),再与标准高斯形式 \(e^{-x^2/(2\sigma^2)}\) 对照,便可看出二次项系数越大,分布越窄,方差越小。这里的二次项系数量级为 \(\kappa\xi^d/(k_BT)\),因此
这里的 \(\sim\) 表示省略了依赖区域形状和平均方式的有限常数。结果适用于相关长度远大于粗粒化单元大小的情形;我们关心的是它如何随 \(\kappa\) 与 \(\xi\) 改变,而不是给某一种形状求精确系数。
这个估计也没有把区域内部当成彼此独立的粒子。若把同一区域分成 \(n\) 个等体积的粗粒化单元,记各单元相对背景的偏离为 \(\delta m_i\),则平均值的方差按定义为
除了每个单元自身的方差,这个和还保留了不同单元之间的交叉关联。只有单元彼此独立时,才能直接丢掉 \(i\ne j\) 的项;而我们选择相关区域,恰恰是因为其中的变化不能视为完全独立。前面的自由能估计通过梯度项保留了这种协调,再求整个区域平均值的涨落,避免了把微观粒子数简单代入独立抽样公式。
现在需要比较的背景 \(m_0\) 与典型偏离 \(\sigma_\xi\) 都已经出现。下一步不再增加一种新模型,而是检验刚才的近似能否成立。
5. Ginzburg 判据与平均场的适用范围¶
现在已有两类答案,一类是平均场预测的有序背景,另一类是围绕它计算出的相关区域涨落。Ginzburg 判据把二者放在同一尺度上比较,用被忽略部分的大小检查原近似是否自洽。沿着这个比较,可以进一步理解维度为何重要,以及平均场失效的临界范围有多宽。
5.1 相对涨落决定近似是否自洽¶
有限温度下存在涨落,并不足以否定平均场近似。真正需要判断的是,典型涨落相对于被当作零阶答案的背景是否很小。这个比较既有直接的物理意义,也能够从自由能展开中高阶项的相对大小得到解释。
如果 \(\sigma_\xi\ll |m_0|\),一个相关区域的平均值通常只在背景附近轻微变化。此时,以这个背景作为零阶答案是合理的起点。
如果 \(\sigma_\xi\) 已经与 \(|m_0|\) 同量级,变化就不能继续被称为“小偏离”。均匀极小值仍然可以写出来,却失去了原先的近似依据。用方差与背景平方定义无量纲比值,要求
这就是 Ginzburg 判据的核心[5]。分子衡量被忽略的相关区域涨落,分母衡量被保留的有序背景;二者的比较,决定以平均场为起点的展开是否自洽。
为什么用背景作比较,也可以从多项式本身看出。在有序极小值处,\(\kappa=2u m_0^2\)。围绕它展开时,三次项相对于二次项的量级是 \(\delta m/m_0\),四次项则由这个比值的平方控制。因此,当典型偏离已经接近背景,原先忽略的非线性项就不再是小修正。
下面从有序侧估计这个比值,因为该侧有非零 \(m_0\) 可作比较。无序侧并不是没有对应的检验,而是不能直接除以零背景,需要从响应等量的涨落修正来判断。
5.2 四维为什么会自然出现?¶
判据给出了比较的方法,但还需要知道这个比值怎样随临界距离改变。背景、恢复强度与相关体积都在变化,其中相关体积的增长直接依赖空间维度。把这些变化合在一起,四维就会作为不同标度行为的分界出现,而不再只是一个需要记住的数字。
计算之前还要区分临界点位置的改变,因为涨落可能移动临界温度。以下把规则的临界点平移吸收到参数定义中,改用
表示相对于实际临界点的距离。平均场自洽区内的领先估计仍写成 \(m_0^2\sim a|t|/u\)、\(\kappa\sim a|t|\);这里 \(a\) 已允许包含有限的参数重定义。临界点移动与临界指数改变是两件不同的事。 我们现在检验的是后者所涉及的奇异涨落,而不是把临界温度的偏移误当作失效。
把第 4 节的方差估计代入判据,得到
再使用 \(\kappa\sim a|t|\) 和 \(\xi\sim\sqrt{c/(a|t|)}\),便得到
这个式子可以一步一步读。背景平方按 \(|t|\) 下降,恢复强度也按 \(|t|\) 下降;相关体积却按 \(|t|^{-d/2}\) 增大。三种变化合在一起,才产生指数 \((d-4)/2\)。维度不是后来附加的标签,而是通过“参与共同涨落的区域有多大”进入计算。
将有限系数记为 \(\overline g_0\),可以把结果简写为
\(\overline g_0\) 是无量纲的涨落强度前因子,与微观参数、归一化及所选窗口有关,不是所有材料共有的常数。例如,取参考长度 \(\xi_0=\sqrt{c/a}\),便有 \(\overline g_0\sim k_BT_c/(a^2\xi_0^d/u)\)。其中 \(a^2/u\) 具有自由能密度的量纲,乘以体积 \(\xi_0^d\) 才是能量,因此这个比值确实无量纲。
按维度比较,可以区分三种不同的临界趋势。
- \(d<4\) 时,在存在有限温度连续相变的条件下,越接近临界点,长波长相对涨落越大,平均场最终失去自洽性。若系统根本没有这样的相变,就不能把结论理解成只需修正几个指数。
- \(d>4\) 时,长波长奇异修正相对减弱,普通连续相变的领先临界指数可以取平均场值;这不意味着所有短尺度修正都消失。
- \(d=4\) 时,简单幂次计数给不出增大或减小,必须进一步检查对数修正,不能直接宣布涨落无关。
因此,当前这类短程标量四次理论的上临界维度为 \(d_c=4\)。这里的“四维”是理论的维度分界,不是说细菌实验实际位于四维,也不是说每一种相变都有同一个上临界维度。
5.3 Ginzburg 区有多宽?¶
对于四维以下、存在有限温度连续相变的系统,相对涨落在接近临界点的过程中增大,但它在多近的位置才变得不可忽略,还取决于模型参数。Ginzburg 数把这个问题转化为临界温区的估计,也帮助我们理解平均场为什么可以在一段范围内有效,却在更接近临界点时失去可靠性。
在 \(d<4\) 时,令 \(\mathcal R_G\) 达到量级 \(1\),得到 Ginzburg 数
它给出平均场开始失去自洽性的临界温区量级。如果 \(\mathrm{Gi}\ll1\),就可能存在一个窗口
左边保证涨落尚小,右边保证足够接近临界点,低阶 Landau 展开仍然合适。“平均场有用”与“临界点附近平均场失效”因此并不矛盾。实验可能在相当宽的区间看到平均场式行为,只有进入更窄的临界区,才明显看到偏离。
以三维为例,同样取 \(\overline g_0=0.01\)。在 \(|t|=10^{-2}\) 时,\(\mathcal R_G\simeq0.1\);接近到 \(|t|=10^{-4}\),它便达到 \(1\)。这里比较的是同一个近似内部的大小,而不是从一张实验照片估出了 Ginzburg 数。\(\mathrm{Gi}\) 也是交叉尺度,不是一条新的热力学相边界。
还需要注意,相对涨落增大并不要求每一种绝对涨落都发散。三维的这个高斯估计给出
随着 \(\xi\) 增大,我们平均的区域也变大,所以区域平均值的绝对方差可以变小;但背景平方下降得更快,比值仍然增大。真正失去可靠性的,是“背景远大于其自身起伏”这个条件。
5.4 判据给出了边界,重整化群要回答边界之后的问题¶
判定平均场失效,不等于已经算出了正确的临界行为。Landau 理论提供可计算的起点,Ginzburg 判据说明这个起点的适用条件;当条件不再满足时,还需要一种重新组织涨落计算的方法。这一尚未解决的问题,将把我们引向重整化群。
但判据没有给出新的 \(\beta\) 或 \(\nu\)。当 \(\mathcal R_G\) 不再很小时,只补上一次很大的修正也不可靠,因为后续修正可能同样重要。这时需要重新组织计算,逐步处理不同尺度的涨落,同时更新剩余变量之间的相互作用。
这也回应了第 2 讲关于有效参数为什么需要随着观察尺度改变的问题。不是因为我们想给同一个系统换一套符号,而是因为被消去的涨落会改变剩余部分的有效描述。Ginzburg 判据指出平均场的边界;重整化群则提供跨越这个边界的方法。
6. 用四次模型理解细菌涡旋¶
有了自由能、相互作用与涨落的区分,再读开篇的论文,就不必同时理解所有实验和计算细节。我们先辨认通道几何怎样改变涡旋之间的作用,再检查作者的有效模型与本讲平衡理论的联系,最后用一个小规模计算展示其中的基本机制。实验事实、模型近似和教学演示将在这个过程中分别说明。
6.1 通道几何怎样改变有效相互作用¶
开篇两段视频最值得解释的差别,是相邻涡旋的排列偏好随间隙改变。论文把这一变化追溯到边界附近的运动路径,再通过腔室与支柱两类变量构造有效相互作用。理解这条机制链,就能看清为什么只改变通道几何,也可能改变最终的有序排列。
实验发现,正方晶格在间隙约 \(8\,\mu\mathrm m\) 两侧表现出不同的平均近邻旋转关联。在更大的间隙、接近约 \(20\,\mu\mathrm m\) 以后,单个腔室对稳定涡旋的约束又逐渐失效。这两种变化不能混为一谈。前者涉及相邻旋转偏好的改变,后者涉及局部涡旋本身还能否维持[1]。读原文时还需注意,它用 \(\chi\) 表示归一化近邻关联,而本讲第 3 节的 \(\chi\) 是外场响应率;两者不是同一个观测量。
机制的关键在腔室边界附近。细菌不仅在腔室内部形成整体流动,还在边界附近形成沿边运动。窄间隙下,这些边界运动主要留在各自腔室周围,相邻腔室之间的作用偏好相反的整体旋转。间隙加宽后,边界运动的路径改变,可以沿腔室之间的星形支柱周围绕行,进而使围绕同一支柱的腔室偏好同向旋转[1]。
因此,作者的完整模型不只包含腔室旋转 \(V_i\),还包含支柱附近的环流 \(P_j\)。这不是为了让方程更复杂,而是因为能够改变有效相互作用方向的那一部分运动,不能一开始就被省略。 两套变量共同构成模型,参数从实验流场的时间序列中估计,再检验其能否再现实验的近邻关联趋势。
论文随后对支柱环流作近似消去,得到只保留腔室旋转的有效四次模型
符号 \(\langle i,j\rangle\) 表示每对相邻腔室计一次。第一项决定邻居的相对偏好。\(J>0\) 时同号取值降低模型能量,\(J<0\) 时异号取值降低模型能量。后面的局域势决定单个旋转变量的偏好强度及稳定性,其中 \(b>0\)。这里沿用论文的 \(a,b\),其中 \(a\) 是局域二次系数,不同于第 3—5 节中 \(r=a t\) 的温度斜率 \(a\)。
于是,方程中的两部分便与前面的推导对应起来。局域四次势与邻近耦合,分别对应第 2 节的状态竞争和第 4 节的空间协调。 论文中复杂的运动,并没有被解释成一个孤立多项式;多项式是整个相互作用模型的局域部分。
6.2 论文与 Landau 理论的联系及其边界¶
论文的约化模型已经包含熟悉的四次势,但形式相似并不意味着可以照搬全部热力学结论。需要分别辨认局域旋转偏好、集体有序和随机涨落,也要区分有效模型中的近似与真实实验的非平衡性质。这些区别决定了本讲哪些结论能够帮助解释论文,哪些还需要独立检验。
先看局域势
当 \(a<0\) 时,它偏好两个相反方向的非零旋转;当 \(a>0\) 时,它偏好接近零的旋转。这个结构与第 2.2 节的单阱、双阱示意相同。但单个腔室偏好非零旋转,并不自动意味着整个晶格建立了同向或交替有序。 后者还取决于耦合、噪声与空间关联。
即使作最简单的均匀近似,这个区别也能看见。设每个格点有 \(q\) 个邻居,并暂取 \(J>0\)、\(V_i=m\),则每格点的模型能量为
因子 \(q/2\) 来自每条键只能算一次。集体均匀状态的二次系数已经是 \(a-qJ\),不再只是局域的 \(a\)。不过,这仍只是限制到均匀构型的模型能量,不是完成有限噪声统计后得到的精确自由能或临界条件。
再看涨落。原论文的完整随机动力学中,腔室与支柱具有不同的有效噪声强度,数据推断也不支持把二者视为同一个热平衡温度。消去支柱后的简化描述具有平衡式概率权重,但它是一个有效近似,不能据此宣称活细菌回到了普通热平衡[1]。
这里涉及的两次“平均场”也需要分清。论文对支柱变量作平均场式消去后,仍然保留了许多相互作用的 \(V_i\) 及其随机变化;第 2 节的均匀极小值近似则进一步把场的统计问题压缩为一个背景。前一种简化并没有自动完成后一种求解。
因此,从复杂流动中选出标量变量、为其写出对称的四次势,再通过相互作用和涨落检验模型,构成了这篇论文与本讲最直接的联系。论文中的有序区域和翻转也提醒我们,平均状态之外还有需要处理的结构。但仅凭这些画面,不能测出 \(\mathrm{Gi}\),不能证明上临界维度为四,也不能把间隙处的关联换号直接当作本讲推导的热平衡连续相变。
6.3 两个相邻涡旋的复现¶
为了把局域偏好、相互作用与涨落联系起来,我们从论文的有效模型中只保留两个相邻腔室。程序输出一张动态图和一张统计图。动态图展示涡旋怎样旋转,统计图则借鉴论文图 1j、k 的思路,分别量化相邻涡旋的方向关联与旋转强度[1]。从观察运动到计算统计量,我们便能检查图像中的直觉是否得到支持。这里复现的是基本机制,不是原文实验曲线或完整晶格模拟。
两个变量的模型能量为
在这个简化的平衡式模型中,指定同一个有效噪声尺度 \(\Theta>0\),概率密度为
\(Z_2\) 是使全部概率加起来等于 \(1\) 的归一化因子。\(\Theta\) 用来控制模型中的随机分散程度,不是培养液的实际温度。这个概率写法与第 3 讲的配分函数结构相同;在这里,它的适用条件是我们明确选择了这个简化模型,而不是直接从活细菌的非平衡性推出了 Boltzmann 分布。
程序固定 \(a=-1\)、\(b=1\)、\(\Theta=0.3\),动画比较 \(J=-0.5\) 与 \(J=0.5\),统计图则在 \(-0.8\le J\le0.8\) 内扫描耦合。横轴直接使用模型参数 \(J\),不换算为通道宽度,因为这里没有重新拟合实验的几何与参数关系。
要量化“同向还是反向”,先看两个变量乘积的符号。乘积为正表示同向,为负表示反向。对这里只有一对邻居的情形,论文的归一化近邻关联简化为
其中 \(\operatorname{sgn}\) 对正数取 \(+1\),对负数取 \(-1\),在零点取 \(0\)。因此,\(C\) 是同向出现的概率与反向出现的概率之差,越接近 \(+1\) 表示越偏好同向,越接近 \(-1\) 表示越偏好反向。这里用 \(C\) 而不沿用论文的 \(\chi\),以免与前面表示外场响应率的符号混淆。
方向关联还不能告诉我们旋转有多强。为此,另取两个腔室的均方根旋转强度
平方使正反旋转都贡献正值,不会因为整体翻转而抵消。\(C\) 描述邻居之间的关系,\(V_{\mathrm{rms}}\) 描述局部旋转的强度,两者不能互相替代。 统计图的实线通过联合概率密度求这两个平均量,而不是从一段短动画的访问次数估计它们。
要让同一模型运动起来,还需要规定变量怎样随时间改变。最简单的选择是让它一方面沿着模型能量降低的方向变化,另一方面不断受到随机扰动。在一小段模型时间 \(\Delta\tau\) 内,取
第一行来自 \(-\partial H_2/\partial V_i\)。其中 \(-aV_i-bV_i^3\) 把局部旋转拉向四次势偏好的取值,\(JV_j\) 则让相邻涡旋参与决定变化方向。第二行是随机扰动,\(\zeta_1,\zeta_2\) 在每一步从标准正态分布中独立抽取,均值为零、方差为一。它让变量不会永远停在一个最低点上。
随机增量的方差为 \(2\Theta\Delta\tau\),所以振幅与 \(\sqrt{\Delta\tau}\) 成正比,而不是与 \(\Delta\tau\) 成正比。这样,把许多独立的小步合起来时,累积方差随经历的时间增长,不会因为把计算步长缩小就人为抹去噪声。
这里把阻尼系数吸收到时间单位中,\(\tau\) 不是实验秒数。这是一个最简单的随机松弛模型,也称为过阻尼 Langevin 动力学。噪声振幅中的 \(\Theta\) 与概率密度里的 \(\Theta\) 相同,使连续时间模型具有上面的稳态权重;代码采用足够小的时间步近似它。动画与统计图不是两个拼接起来的例子,而是同一能量函数的时间图像与统计图像。
统计图还用圆点给出随机动力学的计算结果,作为对实线的另一种检验。每个耦合值模拟 \(2048\) 组独立的涡旋对,先经历 \(60\) 个模型时间单位的初始演化,再各取一个末态计算统计量;积分步长取 \(0.0025\)。误差条表示这些独立样本给出的标准误,其中均方根的误差通过平方根函数传播。这样,样本量来自独立模拟,而不是把动画中彼此关联的相邻帧当成独立测量。
为了把变量的变化看得更清楚,动态图采用与论文补充视频相似的并列视图。左列让固定的灰度示踪纹理随旋转场连续运动,右列在同一画面上显示模型速度箭头。示意场的角速度与 \(V_i\) 成正比,因此两种视图由同一条模型轨迹驱动;纹理不是每一帧重新生成的噪点,也不是实验中的细菌影像。这里只改变变量的呈现方式,不另外引入一个细菌运动模型。
完整代码见 5.vortex_tutorial.py。安装 numpy、matplotlib 和 Pillow 后,运行 python 5.vortex_tutorial.py,图片输出在脚本所在目录。
# 第5讲 两个相邻涡旋的模型与可视化
# 同一个能量函数用于计算联合概率分布,并驱动含噪声的旋转轨迹。
import json
import os
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
ROOT = Path(__file__).resolve().parent
OUTPUT = ROOT / "vortex_reproduction"
# ============================================================
# 第一部分 局域四次势与两个涡旋的联合概率
# ============================================================
# v 表示有正负的旋转强度。a 决定原点附近的曲率,b > 0 保证大振幅时稳定。
# 默认 a = -1、b = 1 时,孤立涡旋的局域势在 v = -1 和 v = 1 处有最低点。
def local_potential(v, a, b=1.0):
return 0.5 * a * v**2 + 0.25 * b * v**4
def pair_density(j, a=-1.0, b=1.0, noise=0.3, extent=2.6, points=401):
# j 对应耦合 J,noise 对应正文的有效噪声尺度 Theta,不是培养液温度。
# 在 [-extent, extent] 内取 points 个网格点,枚举两个旋转变量的组合。
if b <= 0 or noise <= 0 or extent <= 0 or points < 3:
raise ValueError("Require b, noise, extent > 0 and points >= 3")
v = np.linspace(-extent, extent, points)
v1, v2 = np.meshgrid(v, v)
# H_2 = U(V_1) + U(V_2) - J*V_1*V_2,耦合符号决定偏好同号还是异号。
energy = local_potential(v1, a, b) + local_potential(v2, a, b) - j*v1*v2
# 减去最低能量只改变所有权重的公共倍数,不改变归一化后的概率。
# 这样指数不大于零,可避免大正指数造成的数值溢出。
weight = np.exp(-(energy - energy.min()) / noise)
# 用二维梯形积分近似归一化因子,每个方向的端点权重取一半。
# measure 包含网格面积,因此应让 density * measure 的总和等于 1,
# 而不是直接让所有网格点的概率密度相加等于 1。
quadrature = np.ones(points)
quadrature[[0, -1]] = 0.5
measure = np.outer(quadrature, quadrature) * (v[1] - v[0])**2
density = weight / np.sum(weight * measure)
return v, density, measure
# ============================================================
# 第二部分 用随机动力学生成旋转轨迹
# ============================================================
# Euler-Maruyama 方法把一次变化拆为能量梯度驱动的漂移与随机增量。
# dt 是模型时间步长;每经过 substeps 步记录一帧,先丢弃 burn_steps 步,
# 以减小从零状态出发的影响。这不保证一段有限轨迹已访问全部稳态区域。
def simulate_pair(j, frames=240, substeps=20, dt=0.005, burn_steps=10000,
seed=2026, replicas=1, a=-1.0, b=1.0, noise=0.3):
if min(frames, substeps, replicas) < 1 or dt <= 0 or burn_steps < 0:
raise ValueError("Require positive sizes and dt, and nonnegative burn-in")
if b <= 0 or noise < 0:
raise ValueError("Require b > 0 and noise >= 0")
# 固定种子使演示可重复;replicas 表示并行模拟的独立涡旋对数量。
rng = np.random.default_rng(seed)
state = np.zeros((replicas, 2))
# path 的三个维度依次是记录帧、涡旋对编号、两个腔室的旋转变量。
path = np.empty((frames, replicas, 2))
for step in range(burn_steps + frames * substeps):
# state[:, ::-1] 交换两列,为每个涡旋取出邻居的旧状态。
# 两个变量都从同一时刻的状态计算漂移,避免先后更新引入人为差异。
drift = -a*state - b*state**3 + j*state[:, ::-1]
# 随机增量的方差为 2*noise*dt,因此振幅必须随 sqrt(dt) 缩放。
state += dt*drift + np.sqrt(2*noise*dt)*rng.normal(size=state.shape)
if step >= burn_steps and (step-burn_steps+1) % substeps == 0:
path[(step-burn_steps)//substeps] = state
if not np.isfinite(path).all():
raise RuntimeError("Unstable trajectory; reduce the integration step")
return path
# ============================================================
# 第三部分 计算并绘制方向关联与旋转强度
# ============================================================
def audit_layout(fig, name, exclude_axes=()):
fig.canvas.draw()
fig.set_layout_engine("none")
# 可选的排版检查默认关闭,普通运行不需要额外的检查模块。
if os.environ.get("VORTEX_FIGURE_AUDIT") == "1":
from audit_panel_alignment import require_matplotlib_panel_alignment
require_matplotlib_panel_alignment(
fig, json_out=str(OUTPUT / f"{name}.alignment.json"),
exclude_axes=list(exclude_axes), strict=True,
)
def save_figure(fig, name, exclude_axes=()):
audit_layout(fig, name, exclude_axes)
fig.savefig(ROOT / f"{name}.png", dpi=300)
fig.savefig(OUTPUT / f"{name}.pdf")
fig.savefig(OUTPUT / f"{name}.svg")
plt.close(fig)
print(f"Saved image: {name}.png")
def equilibrium_statistics(j, **kwargs):
v, density, measure = pair_density(j, **kwargs)
v1, v2 = np.meshgrid(v, v)
# 把概率密度乘以积分面积,得到各个网格点的统计权重。
weight = density * measure
# 对单独一对邻居,论文的归一化关联简化为 sign(V1*V2) 的概率平均。
# 平方后再求均方根,则只测旋转强弱,不会因方向相反而抵消。
alignment = np.sum(np.sign(v1*v2) * weight)
rms = np.sqrt(np.sum(0.5*(v1**2+v2**2) * weight))
return np.array([alignment, rms])
def sample_statistics(j, replicas=2048, seed=2026, dt=0.0025, burn_time=60):
if replicas < 2 or dt <= 0 or burn_time < 0:
raise ValueError("Require replicas >= 2, dt > 0 and burn_time >= 0")
# 每组独立涡旋对经历初始演化后只取一个末态,不把相邻动画帧当成独立样本。
state = simulate_pair(j, frames=1, substeps=1, dt=dt,
burn_steps=round(burn_time/dt),
replicas=replicas, seed=seed)[0]
alignment = np.sign(state[:, 0]*state[:, 1])
square = np.mean(state**2, axis=1)
rms = np.sqrt(square.mean())
means = np.array([alignment.mean(), rms])
# square 是每组涡旋对的平均平方强度,样本单位是涡旋对而不是单个腔室。
# 平方根的导数为 1/(2*rms),由此把 square 均值的标准误传播到均方根。
errors = np.array([alignment.std(ddof=1), square.std(ddof=1)/(2*rms)])
return means, errors/np.sqrt(replicas)
def plot_statistics():
# 实线由概率积分计算;七个耦合值另外运行随机动力学,得到点和误差条。
# 两种计算没有互相拟合,也没有使用实验数据。
couplings = np.linspace(-0.8, 0.8, 161)
theory = np.array([equilibrium_statistics(j) for j in couplings])
sampled_j = np.array([-0.8, -0.5, -0.25, 0, 0.25, 0.5, 0.8])
samples = [sample_statistics(j, seed=3100+i) for i, j in enumerate(sampled_j)]
means, errors = (np.array(items) for items in zip(*samples))
fig, axes = plt.subplots(1, 2, figsize=(7.08, 3.8), sharex=True)
fig.subplots_adjust(left=0.11, right=0.97, bottom=0.19, top=0.78, wspace=0.42)
for k, ax in enumerate(axes):
curve, = ax.plot(couplings, theory[:, k], color="#356B91", lw=2,
label="Steady-state integral")
points = ax.errorbar(sampled_j, means[:, k], yerr=errors[:, k],
fmt="o", color="#BA553D", ms=4, capsize=2.5,
elinewidth=1, label="Simulation (mean +/- SE)")
ax.set(xlabel=r"Coupling, $J$", xlim=(-0.87, 0.87), xticks=[-0.8, 0, 0.8])
ax.annotate("ab"[k], xy=(0, 1), xycoords="axes fraction", xytext=(-14, 12),
textcoords="offset points", weight="bold")
axes[0].axhline(0, color="0.75", lw=0.8, zorder=0)
axes[0].set(ylabel=r"Rotation alignment, $C$", ylim=(-1.12, 1.12),
yticks=[-1, -0.5, 0, 0.5, 1])
# 无噪声时,能量最低点给出 |V1| = |V2| = sqrt((|J|-a)/b)。
# 灰色虚线只保留这个最优幅度,蓝色实线则保留有限噪声下的概率分散。
minimum, = axes[1].plot(couplings, np.sqrt(1+np.abs(couplings)), "--",
color="0.45", lw=1.3, label="Zero-noise minimum")
axes[1].set(ylabel=r"R.m.s. rotation, $V_{\mathrm{rms}}$", ylim=(0.85, 1.42),
yticks=[0.9, 1.05, 1.2, 1.35])
fig.legend(handles=[curve, points, minimum], loc="upper center",
bbox_to_anchor=(0.53, 0.99), ncol=2, fontsize=9,
columnspacing=1.6, handlelength=2.5)
save_figure(fig, "vortex_pair_statistics")
# 返回数值用于可选的复核记录;默认运行仍只输出图片,不输出数据表格。
return {"couplings": couplings.tolist(), "theory": theory.tolist(),
"sampled_couplings": sampled_j.tolist(), "means": means.tolist(),
"standard_errors": errors.tolist(), "replicas_per_point": 2048,
"dt": 0.0025, "burn_time": 60, "seeds": list(range(3100, 3107)),
"parameters": {"a": -1.0, "b": 1.0, "noise": 0.3},
"data_origin": "Teaching-model simulation; not experimental data"}
# ============================================================
# 第四部分 把旋转变量转换为连续运动的涡旋画面
# ============================================================
# 这一部分只负责可视化,不求解细菌运动或流体方程。
# 灰度纹理与速度箭头由同一条 V_i 轨迹驱动,不额外引入随机动力学。
def texture_geometry():
x, y = np.meshgrid(np.linspace(-2.04, 2.04, 448),
np.linspace(-1.04, 1.04, 228))
# site 区分左右腔室,radius 与 angle 是相对各自腔室中心的极坐标。
site = (x >= 0).astype(int)
local_x = x - np.where(site == 0, -0.98, 0.98)
radius = np.hypot(local_x, y)
angle = np.arctan2(y, local_x)
# 构造腔室和通道的灰度轮廓;mask 让运动纹理在腔室边缘平滑淡出。
wall_distance = np.minimum(radius-0.96,
np.maximum(np.abs(x)-0.14, np.abs(y)-0.18))
wall = (0.72 + 0.16*np.exp(-((wall_distance-0.026)/0.018)**2)
- 0.16*np.exp(-((wall_distance+0.006)/0.012)**2))
mask = np.clip((0.954-radius)/0.012, 0, 1)
return radius, angle, site, wall, mask
def tracer_textures(seed=71):
# 为两种耦合情形、每种情形的两个腔室各生成一张固定的合成纹理。
# 傅里叶滤波压低高频噪点,使示踪纹理平滑;不是逐帧生成新噪声。
rng = np.random.default_rng(seed)
wave = 2*np.pi*np.fft.fftfreq(256)
k2 = wave[:, None]**2 + wave[None, :]**2
noise = rng.normal(size=(2, 2, 256, 256))
texture = np.fft.ifft2(np.fft.fft2(noise)*np.exp(-0.65*k2)).real
texture /= texture.std(axis=(-2, -1), keepdims=True)
return texture
def texture_frame(phases, texture, geometry):
radius, angle, site, wall, mask = geometry
# phases 是旋转变量随模型时间的累积量。对当前像素反向追踪旋转前的位置,
# 再从固定纹理取值,便能显示连续运动。径向因子控制示意场的角速度分布。
theta = angle - (1.2-0.08*radius**2)*phases[site]
px = np.clip(127.5 + 122*radius*np.cos(theta), 0, 254.999)
py = np.clip(127.5 + 122*radius*np.sin(theta), 0, 254.999)
ix, iy = px.astype(int), py.astype(int)
fx, fy = px-ix, py-iy
# 采样位置通常不落在整数像素上,用周围四个像素作双线性插值,避免跳动。
value = ((1-fx)*(1-fy)*texture[site, iy, ix]
+ fx*(1-fy)*texture[site, iy, ix+1]
+ (1-fx)*fy*texture[site, iy+1, ix]
+ fx*fy*texture[site, iy+1, ix+1])
return np.clip(wall*(1-mask) + (0.52+0.11*value)*mask, 0, 1)
def animate_vortices():
frames, substeps, dt = 240, 20, 0.005
spacing = substeps * dt
paths = np.stack([simulate_pair(j, frames=frames, substeps=substeps,
dt=dt, seed=2026+k)[:, 0]
for k, j in enumerate([-0.5, 0.5])])
# 对旋转强度作时间累加,近似积分得到纹理运动需要的旋转相位。
phases = np.cumsum(paths, axis=1)*spacing
geometry, textures = texture_geometry(), tracer_textures()
fig, axes = plt.subplots(2, 2, figsize=(9.6, 5.6), facecolor="white")
fig.subplots_adjust(left=0.025, right=0.985, bottom=0.065, top=0.885,
wspace=0.035, hspace=0.22)
fig.text(0.26, 0.965, "Synthetic tracers", ha="center", fontsize=12)
fig.text(0.75, 0.965, "Model velocity", ha="center", fontsize=12)
fig.text(0.035, 0.018, "Simulation", fontsize=10, color="0.3")
clock = fig.text(0.965, 0.018, "", ha="right", fontsize=10, color="0.3")
qx, qy = np.meshgrid(np.arange(-0.7, 0.71, 0.28), np.arange(-0.7, 0.71, 0.28))
keep = qx**2+qy**2 < 0.76**2
qx, qy = qx[keep], qy[keep]
images, arrows = [], []
# 上下两排比较负耦合与正耦合;左右两列显示纹理和叠加速度箭头的画面。
for case, j in enumerate([-0.5, 0.5]):
pair_images, pair_arrows = [], []
for column, ax in enumerate(axes[case]):
letter = "abcd"[2*case+column]
ax.set_title(letter + (fr" $J={j:+.1f}$" if column == 0 else ""),
loc="left", fontsize=11, pad=5)
ax.set_axis_off()
pair_images.append(ax.imshow(np.zeros_like(geometry[0]), origin="lower",
cmap="gray", vmin=0, vmax=1, extent=(-2.04, 2.04, -1.04, 1.04)))
for center in [-0.98, 0.98]:
pair_arrows.append(axes[case, 1].quiver(center+qx, qy, 0*qx, 0*qy,
color="#A32A28", angles="xy", scale_units="xy", scale=6,
width=0.003, headwidth=3.2, headlength=4.2, pivot="mid"))
images.append(pair_images)
arrows.append(pair_arrows)
def update(frame):
for case in range(2):
gray = texture_frame(phases[case, frame], textures[case], geometry)
images[case][0].set_data(gray)
images[case][1].set_data(0.78+0.32*(gray-0.6))
for site, arrow in enumerate(arrows[case]):
# 圆周运动的速度为 (-omega*y, omega*x),与纹理使用相同的角速度。
omega = (1.2-0.08*(qx**2+qy**2))*paths[case, frame, site]
arrow.set_UVC(-omega*qy, omega*qx)
clock.set_text(fr"$\tau={(frame+1)*spacing:04.1f}$")
update(0)
audit_layout(fig, "vortex_pair_dynamics")
width_pixels = 648
fig.set_dpi(width_pixels/fig.get_figwidth())
shades = np.linspace(0, 255, 16)
blend = np.linspace(0, 1, 8)[:, None]
colors = np.vstack([np.repeat(shades[:, None], 3, axis=1),
(1-blend)*[163, 42, 40]+blend*240])
palette = Image.new("P", (1, 1))
palette.putpalette(colors.astype("uint8").ravel().tolist() + [0]*(768-colors.size))
movie = []
for frame in range(frames):
update(frame)
fig.canvas.draw()
rgb = Image.fromarray(np.asarray(fig.canvas.buffer_rgba())[:, :, :3])
movie.append(rgb.quantize(palette=palette, dither=Image.Dither.NONE))
movie[0].save(ROOT / "vortex_pair_dynamics.gif", save_all=True,
append_images=movie[1:], duration=50, loop=0, optimize=True)
if (ROOT / "vortex_pair_dynamics.gif").stat().st_size >= 10_000_000:
raise RuntimeError("GIF exceeds the 10 MB limit; reduce width_pixels")
update(frames//2)
fig.savefig(OUTPUT / "vortex_pair_dynamics.pdf")
fig.savefig(OUTPUT / "vortex_pair_dynamics.svg")
plt.close(fig)
print("Saved animation: vortex_pair_dynamics.gif")
# ============================================================
# 第五部分 统一绘图风格并输出结果
# ============================================================
def main():
OUTPUT.mkdir(exist_ok=True)
plt.rcParams.update({
"font.family": "sans-serif", "font.sans-serif": ["DejaVu Sans"],
"font.size": 11, "svg.fonttype": "none", "pdf.fonttype": 42,
"axes.spines.top": False, "axes.spines.right": False,
"axes.linewidth": 0.8, "legend.frameon": False,
"figure.facecolor": "white",
})
statistics = plot_statistics()
if os.environ.get("VORTEX_FIGURE_AUDIT") == "1":
(OUTPUT / "statistics_data.json").write_text(json.dumps(statistics, indent=2))
animate_vortices()
if __name__ == "__main__":
main()
观看时,先沿左列比较上下两排,再用右列的箭头确认每个腔室的转向。负耦合使相邻涡旋偏好相反的方向,正耦合使它们偏好相同的方向;但这些偏好并不是把两者锁死的机械齿轮。旋转强弱始终变化,短时也可以偏离更常见的相对排列。相互作用决定什么更有利,噪声使实际状态围绕这种偏好不断变化。
一段短动画只能展示一次有限时间的运动,不能保证充分经历所有可能状态。下面把观察转成统计问题,用概率积分得到的曲线与独立随机模拟得到的点,比较不同耦合下的方向关联和旋转强度。
先看左图的关联换号。 当 \(J<0\) 时,反向旋转降低耦合能,\(C<0\);当 \(J>0\) 时,同向旋转更有利,\(C>0\)。在 \(J=0\) 处,两个腔室各自运动,同向与反向没有统计偏好,故 \(C=0\)。这把动画中的排列倾向变成了一条可计算的曲线,但有限的两个变量只产生平滑变化,不能把曲线穿过零点叫作热力学相变。
再看右图为什么左右对称。 把 \(J\) 换成 \(-J\),同时把第二个涡旋反转,模型能量保持不变。因此,同向与反向的偏好互换,而 \(V_1^2\)、\(V_2^2\) 不变,均方根旋转强度也不变。只测旋转强弱,我们就无法区分这两种排列。这与第 1 节的序参量问题相呼应,也解释了为什么论文需要分别统计局部旋转与邻居之间的关联。
右图还保留了一条只看能量最低点的灰色虚线。在有利的相对方向上,令两个涡旋的幅度均为 \(v\),则 \(H_2=(a-|J|)v^2+bv^4/2\)。求极小值得到
这是忽略噪声时的最优幅度,不是对整个概率分布求得的均方根。在当前参数范围内,蓝色实线低于这条虚线,因为有限噪声使变量也访问最低点之外的取值。找到能量最低点,与算出有涨落时的统计平均,并不是同一件事。 这个差别直接连回本讲的主线,但还不是对空间涨落的 Ginzburg 检验。
两个腔室不能产生热力学极限中的真正相变,也不能用来提取本讲的临界指数或 Ginzburg 数。这个小计算的任务,是把局域偏好、相互作用与概率涨落三件事分别展示出来。只有先分清它们,才有必要继续处理原论文的双晶格随机动力学与参数推断。
6.4 从选择变量到检验近似¶
细菌涡旋算例把本讲的几个核心问题放在了同一个模型中。局域四次势描述单个旋转变量偏好什么取值,相互作用决定邻居之间怎样配合,随机运动则显示系统不会始终停在最有利的状态上。理解这三种作用之后,还需要分清三个层次——如何建立模型、如何近似求解,以及如何判断近似是否可靠。它们分别对应 Landau 理论、本讲的平均场求解与 Ginzburg 判据。
Landau 理论首先提供建立有效描述的方法。 不必追踪每个细菌,而是选择能够辨认所研究有序的变量,再用对称性与稳定性约束势函数。在涡旋模型中,正反旋转的等价性对应偶次项,正的四次项限制旋转幅度,邻近耦合则使局部偏好形成相互关联的排列。这里与 Landau 理论相通的是这套建模思路,而不是“凡是出现四次多项式,就已经完成了平均场计算”。
平均场近似进一步规定怎样求解模型。 第 2 节用一个均匀背景 \(m_0\) 描述系统,通过自由能极小值求出它,却没有计算空间涨落对这个答案的反馈。相比之下,刚才的两涡旋程序直接保留两个变量的联合分布,并不因为计算了统计平均就变成平均场。统计图中的灰色虚线只取两涡旋的无噪声最优幅度,蓝色实线则保留有限噪声下的概率分散;两者的差别帮助我们看见“只保留最优状态”遗漏了什么,但灰色虚线本身不是完整晶格的平均场解。
Ginzburg 判据要检验的,则是这些被忽略的涨落能否继续作为小修正。 动画中存在波动,并不足以说明平均场失效;关键是波动相对于平均场背景有多大。回到具有空间延展的平衡参照模型,还必须考虑哪些位置一起涨落,并在相关长度 \(\xi\) 决定的区域内比较 \(\sigma_\xi\) 与 \(|m_0|\)。只有这样,才从“看见涨落”走到了第 5 节的自洽性检验。
这里也要分清统计图的误差条与 Ginzburg 判据中的涨落。误差条表示我们把统计平均估计得有多准确,\(\sigma_\xi\) 则表示相关区域内的平均序参量本身起伏得有多大。 增加独立模拟的次数可以缩短误差条,却不会因此消除系统自身的涨落。所以,模拟点接近积分曲线,检查的是两种求解方式对同一模型的相容性,并不是已经证明平均场可靠。
由此再回看论文,四次势、相互作用与随机动力学就不再是彼此分离的内容。它们依次让我们描述状态的竞争、理解有序如何建立,并看见最优状态之外仍有统计结构。两涡旋算例把这种区别具体展示出来;要进一步使用 Ginzburg 判据,则需要恢复空间关联,检查这些涨落是否足以改变以均匀背景为起点的预测。
总结¶
回到引言的问题,临界行为能否从一个具体理论中算出来?本讲给出的答案包含两个相接的部分。我们能够通过一个简单的有效自由能算出相变与临界指数,但还必须检查求解时忽略的涨落是否足够小。得到一个答案,与知道这个答案何时可信,是同一条推理链上的两个任务。
Landau 理论提供描述状态竞争的自由能。 先选定所研究的有序,再由对称性决定允许出现哪些项,由稳定性决定哪些非线性项不能省略。最简单的四次理论由此建立起来。它的意义不只是把自由能画成一条曲线,而是让系统为什么偏好零序参量或非零序参量成为可以计算的问题。
本讲的平均场近似把这个问题简化为均匀极小值的求解。 最低点的位置给出有序程度,附近的曲率决定响应强弱,平衡分支的自由能给出热力学变化。平均场临界指数就是沿着这条链算出来的。因此,Landau 理论与平均场近似不能混为一谈,前者规定有效描述的结构,后者是我们最先采用的求解方式。
恢复空间依赖,才有办法检查这个求解方式。 Landau–Ginzburg 自由能同时保留局域恢复代价与邻近区域之间的协调代价,二者共同决定相关长度。接近临界点时,势阱变平、响应增强、相关区域扩大,原先被均匀描述省略的空间涨落便可能变得不可忽略。
Ginzburg 判据把“能否忽略”变成一个相对大小的比较。 在有序侧,以相关区域内平均序参量的方差除以背景平方,要求
它问的不是有没有涨落,而是涨落是否小到足以维持以平均场背景为起点的展开。对本讲的普通短程标量四次理论,四维以下足够接近临界点时,这个条件一般不再成立。平均场的失败,并不意味着序参量和有效自由能失去了意义,而是说明必须改进处理涨落的方式。
细菌涡旋为这条主线提供了可见的例子。四次势与耦合帮助我们解释旋转偏好,动画与统计图让我们区分最优状态和有涨落时的平均;真正的 Ginzburg 检验还要把这种检查推进到相关区域与空间尺度。至此,本讲的逻辑可以归结为
用 Landau 理论写出自由能 → 用平均场近似求出背景 → 恢复空间涨落 → 用 Ginzburg 判据检验近似。
当涨落不再是小修正,接下来的任务不是抛弃有效描述,而是研究它怎样随着观察尺度改变。这正是重整化群要继续回答的问题。下一讲进入伊辛模型的世界——从 1D 精确解到 2D 临界点,通过一个明确的微观模型比较平均场与更精确的答案。随后,块自旋与粗粒化——Kadanoff 的直观 RG 图像将把这里的尺度判断变成实际的变换。
参考文献¶
[1] WIOLAND H, WOODHOUSE F G, DUNKEL J, GOLDSTEIN R E. Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[J]. Nature Physics, 2016, 12: 341–345. DOI: https://doi.org/10.1038/nphys3607
[2] LANDAU L D. On the theory of phase transitions. I[J]. Zhurnal Eksperimental'noi i Teoreticheskoi Fiziki, 1937, 7: 19–32. https://cds.cern.ch/record/480039
[3] LANDAU L D, LIFSHITZ E M. Statistical Physics: Part 1[M]. 3rd ed. Oxford: Pergamon Press, 1980. https://www.sciencedirect.com/book/monograph/9780080230399/statistical-physics
[4] GINZBURG V L. Nobel Lecture: On superconductivity and superfluidity (what I have and have not managed to do) as well as on the “physical minimum” at the beginning of the XXI century[J]. Reviews of Modern Physics, 2004, 76(3): 981–998. DOI: https://doi.org/10.1103/RevModPhys.76.981
[5] TONG D. Statistical Field Theory[EB/OL]. University of Cambridge [2026-09-13]. https://davidtong.org/teaching/statistical-field-theory/











