摘要:节律行为, 即系统行为呈现随时间的周期变化, 在我们的周围随处可见. 不同节律之间可以通过相互影响、相互作用产生自组织, 其中同步是最典型、最直接的有序行为, 它也是非线性波、斑图、集群行为等的物理内在机制. 不同的节律可以用具有不同频率的振子(极限环)来刻画, 它们之间的同步可以用耦合极限环系统的动力学来加以研究. 微观动力学表明, 随着耦合强度增强, 振子同步伴随着动力学状态空间降维到一个低维子空间, 该空间由序参量来描述. 序参量的涌现及其所描述的宏观动力学行为可借助于协同学与流形理论等降维思想来进行. 本文从统计物理学的角度讨论了耦合振子系统序参量涌现的几种降维方案, 并对它们进行了对比分析. 序参量理论可有效应用于耦合振子系统的同步自组织与相变现象的分析, 通过进一步研究序参量的动力学及其分岔行为, 可以对复杂系统的涌现动力学有更为深刻的理解.
关键词: 同步/
序参量/
涌现/
分岔

English Abstract

全文HTML

--> --> -->

节律, 即系统呈现的随时间周期变化, 广泛存在于我们生活的世界, 从钟摆、星体运动到生物钟, 都是不同系统典型的动力学行为. 复杂与非线性系统会表现出多种多样的自组织与群体性行为^[1-3], 而不同的节律之间也可以通过相互作用涌现出各种集体的自组织现象, 其中同步是最典型、最直接的有序行为, 它也是诸如非线性波、时空斑图乃至各种生物集群行为(例如鸟群、鱼群、蜂群、蚁群、人类社会等)的内在物理机制^[4,5].
作为一种最基本的协同现象, 同步研究可以追溯到1673年Huygens发现的关于两个相邻钟摆同步的讨论. 在后来的很多不同物理问题中, 人们都发现了同步现象. 电子与无线电工程的发展很大程度上促进了关于同步的研究. 相比于现象观察及实验研究, 同步的理论研究相对要滞后得多, 直到20世纪初关于极限环的研究才得到了突破. 极限环是非线性耗散系统的一类典型时间振荡解和吸引子, 不同时间振荡之间可以通过相互作用产生新的协同. 由于大量看似不同的系统都具有可以用极限环描述的时间振荡, 人们逐渐认识到, 在各种不同同步行为的背后应该具有共同的物理机理^[5]. Wiener^[6]在其专著中关于“脑电波与自组织系统”一章指出, 脑电波的出现是来自于不同频率的锁定(pulling together of frequencies), 而且这种行为与其他诸如萤火虫同步闪光、蟋蟀同声鸣唱等现象有着共同的机制. 而如何在极限环基础上来研究驱动或相互作用的振子间的同步在理论上就成为重要课题.
Winfree^[7]于1967年提出了全局耦合振子模型. 他认识到, 在同步问题上, 极限环之间相互作用起作用的关键自由度是其相位, 因此只需要考虑耦合的相位振子模型即可揭示同步的动力学实质. Winfree模型的动力学方程可写为

${\dot \theta _i} = {\omega _i} + \left( {\sum\limits_{j = 1}^N {X({\theta _j})} } \right)Z({\theta _i}), $

这里用i = 1, 2, ···, N来标记不同振子, {θ_i}代表振子的相位, {ω_i}为振子的自然频率, 它们各不相同, 设符合某个统计分布函数g(ω). Winfree发现, 自然频率分布较窄时振子之间会相互同步^[7]. 这样, 在众多不同物理背景下出现于不同体系中的很多现象就都可以利用耦合相振子系统的同步动力学来得到非常好的解释. 大量相互作用振子出现的整体同步行为是一种典型的集体涌现, 如何解析刻画同步涌现是一个挑战性课题.
Kuramoto^[8-10]考虑振子数目N → ∞ (热力学极限), 自然频率$\{ {\omega _i}\} $为单峰分布g(ω), 振子间相互作用是耦合强度为K的全局性平均场形式, 并取最简单的相位差的正弦函数sin?(?θ)来描述相互作用. 平均场耦合相振子模型可写为

${\dot \theta _i} = {\omega _i} + \frac{K}{N}\sum\limits_{j = 1}^N {\sin ({\theta _j} - {\theta _i})} . $

该模型被后来研究者称为Kuramoto模型^[11]. 我们还引入振子的集体自然频率$\bar \omega = \displaystyle\frac{1}{N}\mathop \sum \nolimits_{i = 1}^N {\omega _i}$. 一般情况下自然频率分布函数设为关于$\bar \omega$的对称分布. 为讨论方便, 通常设$\bar \omega=0$, 分布对称性满足$g\left( \omega \right) = g\left( { - \omega } \right)$.
尽管振子的自然频率各不相同, 有相互作用时(K ≠ 0)各振子的实际振动频率都会相应地从K = 0时的自然频率{ω_i}处发生偏离, 并随耦合强度变化而变化. 可定义振子的平均频率Ω_i为

${\Omega _i} = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_0^T {{{\dot \theta }_i}(t){\rm{d}}t}. $

为描述振子的整体同步情况(相干程度), Kuramoto引入如下的序参量, 将其定义为所有振子相位复函数的平均场:

$z = R\exp ({\rm{i}}\varTheta ) = \frac{1}{N}\sum\limits_{j = 1}^N {\exp ({\rm{i}}{\theta _j})},$

其中复序参量的模R = |z|描述振子的相干性强弱, $\varTheta $为一任意相位.
在耦合强度很弱时, 只有自然频率极其靠近的振子才会同步, 但它们所占的比例几乎可以忽略, 大量振子的Ω_i都不等, 它们在任一时刻的相位都均匀分布于0— 2π之间, 如图1(a)所示, 此时R = 0. 随着耦合强度的增加, 越来越多的振子会同步, 平均频率Ω_i相等, 这些同步的振子相位之间会靠近并保持固定相位关系, 振子不再均匀分布, 如图1(b)所示. 当所有Ω_i都相等时, R就不为零, 表明此时振子之间可保持固定的相位关系. 存在一个临界的耦合强度K_c, 当K ≤ K_c时R = 0; 当K ≥ K_c时, R ≠ 0. 在很强的耦合下, 振子相位会靠得很近, 形成整体的同步大集团, 如图1(c)所示. 大量耦合振子可以通过相互作用克服自然频率不同带来的无序而涌现出同步态实质上就是一种典型的非平衡相变. 在临界K_c处发生的相变在理论上可以处理. Kuramoto利用统计物理学方法和自洽方程成功求解^[9,10], 理论上得到了临界耦合强度K_c, 并得到在K_c附近序参量R的临界行为.
图 1 耦合振子同步示意图　(a) 在耦合强度很弱时, 大量振子不同步, 任一时刻相位均匀分布于0?2π之间; (b)随着耦合强度的增加, 越来越多的振子会同步, 振子不再均匀分布; (c) 在很强的耦合下, 振子相位会靠得很近, 形成整体的同步大集团; (d) 序参量随耦合强度的变化
Figure1. A schematic diagram of synchronization of coupled oscillators: (a) Most of the oscillators are asynchronous and evenly distributed along the circle; (b) with increasing the coupling, more and more oscillators are synchronized and are no longer evenly distributed; (c) under a strong coupling, oscillators form a single synchronous cluster, and the phases of oscillators are close to each other; (d) dependence of the order parameter on the coupling strength

结合具体的物理背景, Kuramoto模型的应用范围和领域大大拓展 ^[11-17]. 人们研究了不同耦合函数^[18,19]、不同耦合网络拓扑结构^{[14,15,20-25]}、脉冲耦合^[26]等情况下的同步, 探讨了同步涌现的玻璃态^[27]、非对称耦合及其行波和驻波波态^[28,29]、惯性效应^[30,31]、阻挫效应^[32-35]、时间延迟效应^[36]、外场驱动效应^[37]、时变耦合^[38]等与实际物理背景密切相关的问题. 以Kuramoto模型为代表的耦合振子系统集体动力学问题时至今日一直是重要的研究热点.
对于大量振子同步的研究有几种不同层次和方法. 利用非线性动力学^[39-41]和同步分岔树方法^[42-44]可以对耦合振子的微观动力学进行研究. 当振子数目$N \gg 1 $时, 微观动力学分析会变得很繁杂, 可以利用统计物理和宏观方法开展研究. 早期Kuramoto基于统计物理学方法和自洽方程理论成功地对同步相变进行了解析研究^[9]. 但自洽方程方法建立在同步态为定态的前提下, 近年来人们发现了大量非定态集体动力学行为, 因此有必要建立一套反映集体动态行为的统计与宏观理论^[45]. 20世纪90年代初, Watanabe与Strogatz发现^[46], 一类具有对称性的耦合振子系统的高维动力学可以通过引入M?bius变换(后称为Watanabe-Strogatz变换, 简称WS变换)精确降维到三维空间, 这意味着系统的部分可积性. Ott与Antonsen^[47,48]提出了序参量拟设理论, 将高维微观动力学研究降维到二维序参量空间加以研究. 人们随后意识到, OA拟设下的二维动力学实际上是WS变换下的三维动力学的进一步降维^[49,50].
耦合振子系统同步的序参量动力学研究大致包括如下几个重要和基础问题. 第一, 序参量是表征整体行为的特征量, 那么对于大自由度的耦合振子系统, 序参量是如何涌现的？第二, 序参量及其涌现的物理意义是什么？第三, 如何用序参量刻画不同的同步或有序状态及其转变？第四, 如果有不止一个序参量共存, 这些序参量之间如何竞争？第五, 对于实际系统, 如何引入恰当的表征集体行为的序参量来刻画有序行为？如何从复杂系统的动态大数据中重构序参量动力学？本文将集中于上述的前三个基本问题, 对大量振子组成的系统产生同步作为一种典型的涌现行为开展宏观和统计层面的研究, 总结、回顾并对比几种不同的序参量理论框架^[51].

首先来看Kuramoto如何通过引入序参量并利用求解自洽方程来解析处理同步问题. 自洽方法的基础是通过假设系统存在一个不随时间改变的定态, 在此状态下序参量为一个待定的定值, 通过序参量的定义和系统定态的运动方程, 得到待定的序参量的值, 并在分析的过程中得到此序参量与相应定态的存在条件. 自洽方法从其方法本身便限定了其适用范围, 虽然只能用于对定态的分析, 却可以不受具体动力学的限制, 是振子系统分析中广泛使用的方法之一.
首先, 假定振子数N足够多情况下, 序参量(4)式与N无关且不随时间变化. 考虑到Kuramoto模型(2)式相互作用的平均场形式, 很容易可以将方程(2)重新写为

${\rm{d}}{\theta _i}/{\rm{d}}t = {\omega _i} + KR\sin (\varTheta - {\theta _i}). $

如果能定出R, 则方程(5)完全可以求解, 但这不是一件简单的事. 一个可行的办法是带着未知R继续讨论, 建立一个关于R的方程来将其求解. 该方程即为自洽方程(self-consistent equation)^[8,51]. 下面的讨论即是围绕这一主题展开.
显然R = 0对应于无相互作用的情形(均匀分布解), 这个非相干态总是系统的一个解, 但不总是稳定. 当$K > {K_{\rm{c}}}$时均匀分布解失稳. 另外一个解是Ω_i都相等时的解, 此时振子之间可以保持固定的相位, R ≠ 0, 所有振子都以集体频率$\bar \omega $转动, $\varTheta = \bar \omega t$. 此解在$K \leqslant {K_{\rm{c}}}$时不稳定. 只要有Ω_i不相等, ${\theta _i}$就总是均匀分布于0—2π之间.
引入变量${\phi _i} = {\theta _i} - \bar \omega t$, (5)式可写为

${\rm{d}}{\phi _i}/{\rm{d}}t = {\omega _i} - \bar \omega - KR\sin {\phi _i}, $

此方程正是过阻尼情形下的单摆方程. 它有两个解. 1) 同步解. 当方程(6)描述的第i个振子满足$\left| {{\omega _i} - \bar \omega } \right| \leqslant KR$时, 该振子的相位?就会保持定值${\phi _i} = {\sin ^{ - 1}}[({\omega _i} - \bar \omega )/(KR)]$, 这意味着所有满足该条件的振子都会以$\bar \omega$频率运动, 它们处于同步状态. 2) 非同步解. 当方程(6)描述的振子i满足$\left| {{\omega _i} - \bar \omega } \right| > KR$时, 则该振子的相位?就会随时间变化, 且凡是自然频率满足该条件的振子都处于非同步状态.
振子数N→∞时, ${\phi _i}$在0—2π之间会形成分布. 设分布函数为P(?, ω, t), 它不仅依赖于?, 还依赖于振子自然频率ω. 平均场由分布函数可表为

$R{{\rm{e}}^{{\rm{i}}\varTheta }} = {{\rm{e}}^{{\rm{i}}\bar \omega t}}\int\nolimits_0^{2{\text{π}}} {\int\nolimits_{ - \infty }^\infty {{{\rm{e}}^{{\rm{i}}\phi }}P(\phi,\omega,t)g(\omega ){\rm{d}}} } \omega {\rm{d}}\phi, $

单振子分布函数满足连续性方程:

$\frac{{\partial P(\phi,\omega,t)}}{{\partial t}} \!=\! - \frac{\partial }{{\partial \phi }}\big\{ [\omega - \bar \omega - KR\sin \phi ]P(\phi,\omega,t) \big\}. $

考虑到振子具有自然频率分布, 如果只考察?的统计分布, 则需要进一步对频率做积分$P(\phi, t) = $$ \displaystyle\int {P(\omega, \phi, t)g(\omega ){\rm{d}}\omega } $. 这些振子相位的分布直接决定了相关的平均量, 例如序参量, 因而很重要的一点就是如何确定分布函数. 根据上面讨论的两类定态解, 可把分布P(?)分解为同步与非同步两部分:

$P(\phi ) = {P_{\rm{s}}}(\phi ) + {P_{{\rm{as}}}}(\phi ). $

同步振子的相位?趋于不动点, 与时间无关, 因此${P_{\rm{s}}}(\phi )$可由自然频率分布得到:

$\begin{split} &{P_{\rm{s}}}(\phi ) = g(\omega )\left| {\frac{{{\rm{d}}\omega }}{{{\rm{d}}\phi }}} \right| = KRg(\bar \omega + KR\sin \phi )\cos \phi,\\ & \qquad \qquad \phi \in \left[ { - \frac{{\text{π}}}{2},\frac{{\text{π}}}{2}} \right]; \end{split}$

非同步振子的相位?_i则随时间变化, 因为?_i随时间变化是非均匀的, 单位时间内探测到?在$\phi \to $? + d?之间的概率反比于相速度$|\dot \phi |:\; p(\phi, \omega ) \propto |{\dot \phi } |^{-1}$. 把运动方程(6)式代入并归一化可得到

$p(\phi,\omega ) = \frac{{\sqrt {{{(\omega - \bar \omega )}^2} - {{(KR)}^2}} }}{{2{\text{π}}\left| {\omega - \bar \omega - KR\sin \phi } \right|}}, $

因此

${P_{{\rm{as}}}}(\phi) = \int\nolimits_{| {\omega - \bar \omega } | > KR} {g(\omega )p(\phi,\omega ){\rm{d}}\omega } .$

令$x = \omega - \bar \omega $, 并考虑到$g(\bar \omega + x) = g(\bar \omega - x)$, 将(11)式代入(12)式积分可写为

${P_{{\rm{as}}}}(\phi ) = \int_{KR}^\infty {\frac{{g(\bar \omega + x)x\sqrt {{x^2} - {{(KR)}^2}} }}{{{\text{π}}[{x^2} - {{(KR\sin \phi )}^2}]}}{\rm{d}}x}, $

将前面引入的平均场用分布写出来, 并利用(9)式可得

$\begin{split} R{{\rm{e}}^{{\rm{i}}\varTheta }}\; &= \int_{ - {\text{π}}}^{\text{π}} {{{\rm{e}}^{{\rm{i}}(\phi + \bar \omega t)}}P(\phi ){\rm{d}}\phi } \\ &= \int_{ - {\text{π}}}^{\text{π}} {{{\rm{e}}^{{\rm{i}}\phi + {\rm{i}}\bar \omega t}}[{P_{\rm{s}}}(\phi ) + {P_{{\rm{as}}}}(\phi )]{\rm{d}}\phi } .\end{split}$

可以看出(14)式中${P_{{\rm{as}}}}(\phi )$部分积分为零, 只有${P_{\rm{s}}}(\phi )$对积分有贡献. 进一步将积分实部和虚部分离, 可得到(注意R为实):

$R = KR\int_{ - \frac{{\text{π}}}{2}}^{\frac{{\text{π}}}{2}} {{{\cos }^2}\phi g(\bar \omega + KR\sin \phi ){\rm{d}}\phi }, $

$0 = KR\int_{ - \frac{{\text{π}}}{2}}^{\frac{{\text{π}}}{2}} {\cos \phi \sin \phi g(\bar \omega + KR\sin \phi ){\rm{d}}\phi }. $

由(16)式可以确定$\bar \omega $ (实际上在原方程中可通过对称性确定). (15)式是一个自洽方程, 由此可以确定R和同步临界耦合强度${K_{\rm{c}}}$. 当$K \geqslant {K_{\rm{c}}}$时, R由零连续变为一个小量, 因此可把$g(\omega )$在$\bar \omega $附近展开:

$g(\bar \omega + KR\sin \phi ) \approx g(\bar \omega ) + \frac{{g''(\bar \omega )}}{2}{(KR)^2}{\sin ^2}\phi + O({R^4}),$

其中$g''(\bar \omega ) = {\left. {{{\rm{d}}^2}g(\omega )/{\rm{d}}{\omega ^2}} \right|_{\omega = \bar \omega }}$. 代入(17)式可以得到

$1 = \frac{{{\text{π}}K}}{2}g(\bar \omega ) - \frac{1}{{16}}{\text{π}}{K^3}{R^3}g''(\omega ) + O({R^4}).$

当$K \to {K_{\rm{c}}}$时, 由于R²项为高阶小量, R→0, 由此得到同步临界耦合强度为

${K_{\rm{c}}} = 2/[{\text{π}}g(\bar \omega )].$

将(19)式代回(18)式可以确定在临界点K_c附近R的行为:

$R \approx \sqrt {\frac{{8g(\bar \omega )(K - {K_{\rm{c}}})}}{{g''(\bar \omega ){K^3}}}} . $

例如对典型的自然频率符合洛伦兹分布$g(\omega ) = {\{ {\text{π}}[{(\omega - \bar \omega )^2} + {\gamma ^2}]\} ^{ - 1}}\gamma $的情况, 临界耦合为K_c = 2γ, 序参量$R = \sqrt {(1 - 2\gamma /K} )$. 可以看出R在K_c附近为典型的连续相变$R \propto {(K - {K_{\rm{c}}})^{1/2}}$, 如图1(d)所示. 需要指出的是, 这是一种非平衡相变, 是系统有序(耦合)压倒无序(自然频率随机分布)的结果.

耦合振子同步微观动力学的研究表明, 大量耦合振子随着耦合强度的增加会经历由部分同步到整体同步的过渡, 相空间维数随同步进程而逐渐降低, 整体同步时, 系统在相空间的复杂运动会落到一个极低维的空间^[42-44]. 降维意味着系统发生同步的时候, 只需要少数变量即可刻画耦合振子系统的同步. 该结果为多振子体系同步的宏观及序参量描述提供了事实基础.
统计物理学建立了从微观到宏观之间的联系, 进一步通过统计定律来计算微观量的统计平均来得到宏观热力学量^[27]. 热力学和统计物理学的思想和方法为处理耦合振子同步转变问题提供了思路. 由大量相互作用振子组成的系统通常是非平衡系统, 同步是大量振子整体动力学从无序向有序的非平衡相变. 非平衡行为研究自20世纪中期发展起来的以耗散结构理论^[52,53]、协同学^[54,55]以及多学科分支形成的自组织理论等为描述众多的非平衡相变现象的共同本质提供了重要依据.
耗散结构理论认为, 处于非平衡状态的系统会在一定范围内维持原有热力学状态, 一直到远离平衡到一定临界值, 热力学分支失稳使系统进入到有序的结构分支. 非平衡相变现象在自然界各个领域有形形色色的表现, 物理表现十分不同, 而内含的数学实质是有明显的规律, 由此形成了以反映扩散方程等为核心的刻画宏观自组织过程的耗散结构理论^[53].
Haken的协同学则关注于由大量自由度构成的复杂系统在外参量的驱动下和在子系统之间的相互作用下如何以自组织的方式在宏观尺度上形成空间、时间或功能有序结构的条件、特点及其演化规律^[54]. 协同学的基本原理是支配原理(slaving principle), 它认为, 协同系统的状态由一组状态参量来描述, 这些状态参量弛豫时间尺度是不相同的, 慢变的线性不稳定模称为慢变量(slow variable), 而快变的线性稳定模称为快变量(fast variable). 当系统接近于发生显著质变的临界点时, 慢模数目会减少为只有一个或少数几个, 这些慢变量可以完全确定系统的宏观行为并表征系统的有序化程度, 故而称为序参量. 而为数众多的快模则由慢模/序参量所支配, 并可将其绝热消去, 由此可以建立少自由度的协同学基本方程. 在序参量方程基础上, 我们就可来研究协同系统的各种非平衡定态/非定态、稳定性及其非平衡相变^[55].
支配原理在数学上即绝热消去, 即可以对快速变化的变量进行平均加以消去, 只保留变化慢(绝热)的变量. 该方法在物理上有着深刻含义. 下面以如下的n维动力学系统为例来简单介绍一下支配原理:

$\dot {{x}} = {{A}}{{x}} + {{B}}({{x}}), $

其中${{x}}\left( t \right) = {({x_1}\left( t \right), {x_2}\left( t \right), \cdots, {x_n}\left( t \right))^{\rm{T}}}$为n维状态矢量, A为$n \times n$常数矩阵, B(x)包含x的二次以上代数式矢量. 设x = 0为方程的解. 设矩阵A的本征值为$\left\{ {{\lambda _i}, i = 1, 2, \cdots, n} \right\}$, 其中本征值按照其实部由大到小排列. 如果所有本征值实部均为负, 即$\left\{ {{\rm{Re}}{\lambda _i} < 0, i = 1, 2, \cdots, n} \right\}$, 那么x = 0为稳定解. 改变参数使A的m个模失稳, 而其他模仍然保持稳定, 即假设该模对应的本征值$\{0 < {\rm{Re}}\lambda _k^u \ll 1$,$ {\rm{Re}}\lambda _i^s < 0, k = 1, 2, \cdots, m, i = 1, 2, \cdots, n - m\}$, 且${\rm{Re}}\lambda _k^u \ll \left| {{\rm{Re}}\lambda _i^s} \right|$. 通过引入一个T矩阵进行如下的线性变换

$\tilde {{A}} = {{{T}}^{ - 1}}{{AT}}, $

将A对角化

$\tilde {{A}} = \left( {\begin{array}{*{20}{c}} {{{{\lambda}} ^u}}&{\bf{0}} \\ {\bf{0}}&{{{{\lambda}} ^s}} \end{array}} \right), $

这里λ^{u, s}分别为m × m与(n – m) × (n – m)对角子矩阵. 则(21)式可改写为

$\dot {{u}} = {\lambda ^u}{{u}} + {{{B}}^u}({{u}},s), \tag{24a}$

$\dot {{s}} = {{{\lambda}} ^s}{{s}} + {{{B}}^s}({{u}},{{s}}), \tag{24b}$

这里x(t) = (u(t), s(t))^T, 其中u(t) = (u₁(t), u₂(t), ···, u_m(t))^T, s(t) = (s₁(t), s₂(t),··· , s_n–m(t))^T, ${{{B}}_u} \!= \!{(B_1^u, B_2^u, \cdots, B_m^u)^{\rm{T}}}, \;{{{B}}_s}\! = \!{(B_1^s, B_2^s, \cdots, B_{n - m}^s)^{\rm{T}}}.$ s(t)为满足(24b)式的快变量, 而u(t)为满足(24a)式的慢变量. 按照支配原理绝热消去, 令(24b)式左边$\dot {{s}} = 0$, 由此n – 1个方程可将n – m维矢量s作为u的函数解出, 并代入(24a)式可得到m维的非线性方程

$\dot {{u}} = {\lambda ^u}{{u}} + {{{B}}^u}({{u}},{{s}}({{u}})). $

一般实际问题中首先失稳的模往往很少, 常常是一两个, 从(24)式到(25)式, 在理论上可以看作是从多变量方程到少数序参量方程的约化, 在动力学上也产生了极大简化, 这是支配原则在讨论相变点附近行为时给出的有益结果.
需要注意的是, 变量的快慢不能理解为其他诸如振荡的快慢, 而是弛豫的快慢, 因此严格来说应该称为快/慢弛豫变量. 慢变模式有时可能本身就是高频振荡, 该模对于扰动响应的弛豫时间长短才是判断该模为快或慢的标准.
序参量有着几何上的意义, 支配原理也密切联系着中心流形定理. 中心流形定理是一种常用的几何降维方法, 它利用流形与对应子空间相切的特性, 求出系统在中心流形上的约化方程. 对于高维动力系统来说, 通过传统的分岔行为很难直接研究其动力系统. 为了更好地抓住所要研究问题的本质, 一般采取中心流形定理等降维措施将其化为低维方程再进行研究.

下面将从耦合振子微观动力学方程出发, 利用统计力学方法建立分布函数方程, 并通过统计平均引入各阶序参量, 建立相应的序参量方程, 并进一步讨论对方程的降维.
2

-->

4.1.广义序参量

从一般形式的全局耦合振子系统运动方程出发:

${\dot \theta _j}(t) = F({{\alpha}},{\theta _j},{{\beta}},{\gamma _j}) + {\xi _j}(t),$

其中j = 1, 2, ···, N, ${\xi _j}(t)$为作用于第j个振子上的随机噪声, 设为振子间无关联的高斯白噪声, 满足$\left\langle {{\xi _i}(t)} \right\rangle = 0$, $\left\langle {{\xi _i}(t){\xi _j}(t')} \right\rangle = 2 D{\delta _{ij}}\delta (t - t')$. ${{\beta}} = \{{\beta _1}, $${\beta _2}, \cdots, {\beta _m}\} $为一组均匀控制参量, 这些参量对所有振子均相同, 如通常考虑振子间相互作用强度相同. ${{\gamma}} = \left\{ {{\gamma _1}, {\gamma _2}, \cdots, {\gamma _N}} \right\}$为一组非均匀控制参量, 它们对不同振子不一样. 例如振子的自然频率通常各不相同, 此时{γ_i = ω_i}.
定义一组广义序参量α = {α_n}, 其分量为如下的n阶序参量:

${\alpha _n} = \frac{1}{N}\mathop \sum \limits_{j = 1}^N \exp ({\rm{i}}n{\theta _j}).$

当n = 1时, α₁即是Kuramoto引入的相干因子. n > 1时α_n为高阶序参量^[56]. 在热力学极限N→∞下, 可以引入振子相位的密度分布函数ρ(γ, θ, t), ρ(γ, θ, t)dθ为一个振子在时刻t相位处于θ→θ + dθ的概率或时刻t相位处于θ→θ + dθ内的振子数密度. 它满足(12)式对应的Fokker-Planck方程:

$ \dfrac{\partial \rho}{\partial t} = - \dfrac{\partial (\rho v)}{\partial \theta} + D \dfrac{{\partial ^2}\rho }{\partial {\theta ^2}}, $

其中速度场v = F(α, θ, β, γ). 在没有外加噪声(D = 0)时, 分布函数方程退化为如下的连续性方程

$ \dfrac{\partial \rho }{\partial t} + \dfrac{\partial (\rho v)}{\partial \theta} = 0. $

在一般情况下振子是非全同的, 即分布函数ρ(γ, θ, t)与γ有关, 系统总的分布函数需要对所有非均匀参量求和, $\rho \left( {\theta, t} \right) = {\displaystyle\sum \nolimits _i}\rho \left( {{\gamma _i}, \theta, t} \right)$. 例如, 如果$\left\{ {{\gamma _i} = {\omega _i}} \right\}$, 则分布函数为

$\rho \left( {\theta,t} \right) = {\sum\limits _i}\rho \left( {{\omega _i},\theta,t} \right) = \int\limits \rho \left( {\omega,\theta,t} \right){\rm{g}}\left( \omega \right){\rm{d}}\omega,$

其中g(ω)为自然频率分布. ρ(θ, t)包含了振子系统集体行为的所有信息. 在同步问题宏观层面, 我们最为关注序参量α = {α_n}. 对于均匀系统即没有非均匀参量γ, 利用分布函数, 序参量可表为如下积分:

${\alpha _n} = \int {\exp ({\rm{i}}n\theta )\rho (\theta,t)} {\rm{d}}\theta . $

如果系统有非均匀性, 则需要同时对非均匀参量γ求和. 例如, 如果非均匀参量是振子自然频率, 则有

$\begin{split}{\alpha _n} ={}&\int {{\alpha _n}(\omega,t)g(\omega ){\rm{d}}\omega } \\={} & \int {\exp ({\rm{i}}n\theta )\rho (\omega,\theta,t)g(\omega ){\rm{d}}\omega } {\rm{d}}\theta. \end{split}$

可以看到, n阶序参量α_n实际上就是exp(inθ)的统计平均, 或称为exp(iθ)的n阶矩. 由于宏观量的各阶矩描述与分布函数描述等价, 可由一方信息推知另外一方的信息. 还可以看到, 广义序参量实际上就是分布函数ρ(θ, t)的傅里叶变换系数.
2 -->

4.2.Ott-Antonsen拟设

下面先考虑全同振子系统的序参量运动方程. 利用序参量的定义(27)式, 对其进行时间求导, 并利用振子运动方程(26)式可以得到

${\dot \alpha _n} = \frac{{{\rm{i}}n}}{N}\mathop \sum \limits_{j = 1}^N \exp ({\rm{i}}n{\theta _j})F({{\alpha}},{\theta _j},{{\beta}} ). $

由于函数F(α, θ, β)是循环变量相位θ 的2π周期函数, 因此可做傅里叶展开:

$F({{\alpha}},{\theta _j},{{\beta}} ) = \mathop \sum \limits_{k = - \infty }^\infty {f_k}({{\alpha}},{{\beta}} )\exp ({\rm{i}}k{\theta _j}),$

将展开式代入运动方程(33)中可得序参量运动方程为

${\dot \alpha _n} = {\rm{i}}n\mathop \sum \limits_{k = - \infty }^\infty {f_k}(\alpha,\beta ){\alpha _{k + n}},~~n=1,2,\cdots$

考虑最简单的耦合形式, 即(34)式的傅里叶分解中只需包含最低阶项:

$\begin{split} F({{\alpha }},\theta ) ={}& {f_1}({{\alpha}} )\exp ({\rm{i}}\theta ) + {f_{ - 1}}({{\alpha}} )\exp ( - {\rm{i}}\theta ) \\& + {f_0}({{\alpha}} ), \end{split}$

则序参量运动方程(35)简化为

${\dot \alpha _n} = {\rm{i}}n[{f_1}({{\alpha}} ){\alpha _{n + 1}} + {\bar f_1}({{\alpha}} ){\alpha _{n - 1}} + {f_0}({{\alpha}} ){\alpha _n}],$

其中n ≥ 0, ${\alpha _{ - n}} = {\bar \alpha _n}$.
广义序参量可看成是相位振子的集体坐标变量, 这意味着相位运动方程(26)通过上述变换可以化为运动方程(35)和(37), 它们都是一组耦合的序参量方程, 处理该序参量方程组的难度等价于相位运动方程(26), 对其的简化需要新的条件.
方程(37)显然存在一个平庸的非相干解α_n ≡ 0, 对应于耦合振子的非同步态. 随着耦合强度的增加, 耦合振子会产生同步, 整体运动在相空间也会塌缩到一个低维空间中. 在广义序参量空间来看, 系统也必然会在低维空间运动. 根据协同学原理, 这些广义序参量中可能只有少数为慢变量, 其余为快变量, 其中慢变的广义序参量会成为系统状态的真正序参量. 由于快变量的变化都依赖于慢变量, 因此在发生同步转变的区域附近, 各阶序参量之间应存在一定的关系, 它们均依赖于慢变的序参量. 一种最简单的可能情形是所有的各阶序参量α_n都依赖于α₁, 不妨设为α_n = G(α₁, n). 将其代入运动方程(37)中并比较每一阶函数的傅里叶展开系数可以得到^[56]

${\alpha _n} = G({\alpha _1},n) = \alpha _1^n. $

则(37)式的多个运动方程可简并为单一方程,

${\dot \alpha _1} = {\rm{i[}}{f_1}({\alpha _1})\alpha _1^2 + {\bar f_1}({\alpha _1}) + {f_0}({\alpha _1}){\alpha _1}]. $

(38)式正是Ott-Antonsen (OA)拟设(ansatz)^[47,48].
2 -->

4.3.Poisson和不变子流形

OA拟设有什么物理意义呢？下面通过讨论分布函数来进行分析.
由于分布函数是2π循环变量相位的函数, 可将其进行傅里叶展开, 而展开系数即为广义序参量α_n:

$\rho (\theta,t) = \frac{1}{{{\rm{2{\text{π}} }}}}\left[ {1 + \sum\limits_{n = 1}^\infty {\left( {{{\bar \alpha }_n}(t){{\rm{e}}^{{\rm{i}}n\theta }} + {\alpha _n}(t){{\rm{e}}^{ - {\rm{i}}n\theta }}} \right)} } \right]. $

一般来说, 知道了各阶傅里叶系数即广义序参量α_n, 就可以利用上述求和来得到分布函数, 通常这需要无穷阶的傅里叶系数.
OA拟设认为相位分布函数ρ的傅里叶展开系数即广义序参量α_n相互之间并不独立, 且满足幂函数关系:

${\alpha _n}(t) = \alpha _1^n(t),\;\;{\bar \alpha _n}(t) = \bar \alpha _1^n(t). $

将其代入(40)式可以得到

$\begin{split} \rho (\theta,t) \!=\!& \frac{1}{{2{\text{π}}}}\bigg\{1 \!+\! \sum\limits_{n = 1}^\infty \big[\bar \alpha _1^n(t)\exp ({\rm{i}}n\theta) \\ & \!+\! \alpha _1^n(t)\exp ( - {\rm{i}}n\theta )\big] \bigg\}, \end{split}$

(42)式右边的求和为幂级数, 可以得到振子分布为泊松和形式

$\rho (\theta,t) = \frac{1}{{2{\text{π}}}}\frac{{1 - {r^2}}}{{1 - 2r\cos (\theta - \varTheta ) + {r^2}}}, $

其中r为序参量α₁(t)的幅度, Θ为集体相位. 因而分布ρ(θ, t)完全由α₁(t)决定.
由于具有关系(38)式的解满足简并的运动方程(39), 则在动力学演化过程中形式(38)一直得到满足. 振子的分布虽然随着系统的演化而改变, 但将始终具有泊松和分布的形式. 如果系统的初始相密度分布满足泊松分布, 那么不管在任何时刻系统的相密度分布将始终保持着这一性质. 序参量关系式(38)式及简并运动方程(39)被称为动力系统的Poisson和不变子流形. 这个不变子流形一个很重要的特点是虽然α₁(t)可以不含时, 也可以含时, 但分布形式随时间演化保持不变. 这一结果将Kuramoto自洽理论仅仅讨论定态的结果拓展至一般情形.
2 -->

4.4.非全同振子系统的OA拟设

对于振子非全同的情形, 若振子自然频率各不相同, 即{γ_i = ω_i}, 设它们满足分布g(ω). 那么当N → ∞时, 引入密度分布函数ρ(ω, θ, t), 广义序参量相应写为(32)式, 其中α_n(ω, t)为ρ(ω, θ, t)的n阶傅里叶展开式, 也可以理解为自然频率在ω → ω + dω之间的局域序参量. 利用连续性方程可以得到如下递归方程:

${\dot \alpha _n} = {\rm{i}}n\mathop \sum \limits_{j = - \infty }^\infty {f_j}\left( {{{\alpha}},{{\beta}},\omega } \right){\alpha _{j + n}}(\omega,t). $

同样如果耦合函数只包含一阶傅里叶系数, 则可以引入如下局域序参量的Ott-Antonsen拟设:

${\alpha _n}(\omega,t) = \alpha _1^n(\omega,t)$

是方程(44)的一组特解. 一阶序参量为

${\alpha _1}(t) = \int {{\alpha _1}(\omega,t)g(\omega ){\rm{d}}\omega }. $

当分布函数g(ω)的形式为ω的有理分式时, 可将ω从实轴延拓到复ω平面中去, 在不引起发散的情形下(|α1(t)| ≤ 1)直接得到α₁(t)的演化方程.

耦合相振子系统在振子数$N \gg 1 $时的微观动力学研究会变得很困难. OA方法给出了一种有效的将高维耦合振子动力学降维的方案, 但高维系统降维到二维的序参量空间来研究是有前提和成立条件的.
OA方法之所以成功, 其关键在于Kuramoto系统本身的微观动力学具有部分可积性. 这要从1994年Watanabe和Strogatz的工作谈起. 他们研究了全局耦合约瑟夫森结方程(与Kuramoto模型系统同类). 研究发现, N维方程中的每一条轨迹只能局限在一个三维子空间中, 这意味着原始的高维微观态可以通过一定方法降维至低维的宏观态^[46], 这就需要他们提出的后来被称为Watanabe-Strogatz(WS)变换的方法来实现. WS变换的提出在一定程度上为人们寻求高维动力系统的低维解提供了方向, 但早期没有引起人们足够重视, 且数学与物理意义均不明确. 2009年, 受Ott与Antonsen工作的启发, Marvel等^[49,50]成功将这类问题的解推广到一般形式, 给出了WS变换的数学意义, 并清晰地给出OA拟设的数学依据.
2

-->

5.1.WS变换

WS变换实际上来自于复数M?bius变换^[57]. 定义复分数的M?bius变换$F:{\mathbb{C}} \to {\mathbb{C}} $为

$F\left( z \right) = \frac{{az + b}}{{cz + d}},~a,b,c,d \in {\mathbb{C}},~ad - bc \ne 0$

M?bius变换有很多很好的性质, 其中之一就是存在且保持从直线到直线、从圆到圆的映射. 由所有满足上述条件的变换函数可以构成一个群, 称为M?bius变换群,

$\begin{split} M\left( {{{\mathbb{C}}_\infty }} \right) ={}& \Big\{ F:{{\mathbb{C}}_\infty } \to {{\mathbb{C}}_\infty }\Big|F(z) \!=\! \frac{{az + b}}{{cz + d}},\\&~~~ad - bc\! \ne \!0 \Big\}. \end{split}$

考虑一个子群, 该子群包含了那些将单位开圆盘上的复数一对一地映射到自身的所有分数线性变换.
将耦合函数(36)式对应的动力学方程重新写为如下形式:

${\dot \theta _j} = f\exp ({\rm{i}}{\theta _j}) + g + \bar f\exp ( - {\rm{i}}{\theta _j}),$

其中j = 1, 2, ···, N, f = f₁(α)为光滑复函数, $\bar f$是其复共轭, g = f₀(α)为实函数, 它们均不依赖于指标j. 方程(49)定义了一个 N 维动力系统. Marvel和Strogatz^[49]指出, 方程(49)的解满足如下M?bius群在复空间内单位圆上的含时M?bius变换M_t:

$\exp \left[ {{\rm{i}}{\theta _j}(t)} \right] = {M_t}\left[ {\exp ({\rm{i}}{\varphi _j})} \right], $

其中含时M?bius变换表达为

${M_t}(w) = \frac{{\exp \left[ {{\rm{i}}\psi (t)} \right]w + \alpha (t)}}{{1 + \bar \alpha (t)\exp \left[ {{\rm{i}}\psi (t)} \right]w}}, $

这里{φ_j, j = 1, 2, ···, N}是系统的一组运动常数, ψ(t)是实参量函数, α(t)是复函数, $\left| {\alpha \left( t \right)} \right| \leqslant 1$, $\bar \alpha \left( t \right)$是α(t)的复共轭. (50)式和(51)式的变换就被称为WS变换. 对任意时刻t都对应一个变换M_t, 这些变换构成的集合满足群的性质, 全体变换的集合{M_t}构成一个代数系统即M?bius变换群. 系统(49)式的演化完全被变换(50)式和(51)式所支配.
2 -->

5.2.三维动力学

方程(49)满足变换(50)式和(51)式, 说明了系统任一时刻的状态{θ_j(t)}可由一组运动不变量{φ_j}在WS变换{α(t), ψ(t)}的作用下完全确定. 而运动常数或不变积分{φ_j}的存在则反映了WS变换下系统的可积性. M?bius群本质上是一个三参量(ψ(t), Reα(t), Imα(t))的李群, 相轨迹{θ_j(t)}被{ψ(t), α(t)}以及{φ_j}唯一地确定.
下面从代数方程出发去推导ψ(t), α(t)的运动方程. 根据方程(50)有

${\theta _j}(t) = - {\rm{i}}\ln {M_t}\left( {\exp \left[ {{\rm{i}}{\varphi _j}(t)} \right]} \right), $

则对时间的导数为

$ \begin{split} {{\dot \theta }_j}(t)\; & = - {\rm{i}}\frac{1}{{{M_t}\left[ {\exp ({\rm{i}}{\varphi _j})} \right]}}\frac{{{\rm{d}}{M_t}\left[ {\exp ({\rm{i}}{\varphi _j})} \right]}}{{{\rm{d}}t}} \\ & = R\exp ({\rm{i}}{\varphi _j}) + \frac{{\dot \psi + {\rm{i}}\bar \alpha \dot \alpha - \alpha ({\rm{i}}{\dot{ \bar \alpha}} - \bar \alpha \dot \psi )}}{{1 - |\alpha {|^2}}} \\ &\quad + \bar R{{\exp (-{\rm{i}}{\varphi _j})}}, \end{split} $

其中

$R \equiv \frac{{{\rm{i}}{\dot {\bar \alpha }}- \bar \alpha \dot \psi }}{{1 - |\alpha {|^2}}}, $

与(49)式对比可得

$f = R = \frac{{{\rm{i}}{\dot {\bar {\alpha}}} - { {\bar\alpha {\dot \psi}}} }}{{1 - |\alpha {|^2}}},\tag{55a}$

$\;\;g = \frac{{\dot \psi + {\rm{i}}\bar \alpha \dot \alpha - \alpha ({\rm{i}}{\dot {\bar {\alpha }}}- \bar \alpha \dot \psi )}}{{1 - |\alpha {|^2}}}. \tag{55b}$

解出$\dot \alpha $和$\dot \psi $, 则有

$\begin{split}& \dot \alpha = {\rm{i}}(f{\alpha ^2} + g\alpha + \bar f), \\ & \dot \psi = f\alpha + g + \bar f\bar \alpha,\end{split} $

因此通过WS变换, 可以将原始的N维相振子方程约化为三维的闭合方程, 而振子数目N可以是有限的, 也可以是无限大.
2 -->

5.3.从WS变换到OA拟设

WS变换是一个从N维相空间向三维相空间的严格动力学变换. 不同于统计方法例如基于序参量的OA拟设, WS方法并不依赖于任何近似条件或者特定状态, 而是将整个系统的动力学完整而严格地投影到低维系统中^[56]. 当振子数N → ∞时, OA拟设下单一复序参量α₁满足的方程为(39)式, 这是一个二维的闭合实方程. 可以看到, WS变换下的三维方程(56)式中第一式即复序参量方程(39), 因此一个自然的问题是, 这两者之间是否存在某种联系？在什么情况下WS变换可以退化到OA拟设？下面来分析一下WS变换和OA拟设这两种降维方案之间的区别和联系.
对于N振子系统(49)式, 通过WS变换可以将序参量α₁(t)重新表示为

$\begin{split}{\alpha _1}(t)\; & = \frac{1}{N}\mathop \sum \limits_{j = 1}^N \exp \left[ {{\rm{i}}{\theta _j}(t)} \right]\\ & = \frac{1}{N}\mathop \sum \limits_{j = 1}^N \frac{{\exp \left[ {{\rm{i}}\psi (t)} \right]\exp ({\rm{i}}{\varphi _j}) + \alpha (t)}}{{1 + \bar \alpha (t)\exp \left[ {{\rm{i}}\psi (t)} \right]\exp ({\rm{i}}{\varphi _j})}}, \end{split}$

从一组不变量{φ_j}出发, 若初始状态已知, 则系统的动力学情况可以由三维方程(56)描述. 取特定均匀分布{φ_j = 2π(j – 1)/N}, 则序参量(57)式可简化为

${\alpha _1}(t) = \alpha (t)(1 + I), $

方程中

$I = \frac{{1 - |\alpha (t){|^{ - 2}}}}{{1 \pm {{\left( {\bar \alpha (t)\exp \left[ {{\rm{i}}\psi (t)} \right]} \right)}^{ - N}}}}, $

其中“–”对应于N为偶数的情况, “+”对应于N为奇数的情况. 当$N \gg 1 $, $I \ll 1 $, 序参量就可以取近似${\alpha _1}(t) \approx \alpha (t)$. 因而, 当运动常数{φ_j}取均匀测度时, α₁(t)与α(t)的方程完全一致, 方程(56)中α(t)和ψ(t)的演化解耦, WS变换的三维方程流形退化到二维OA流形(39)式.
需要指出的是, 当运动常数{φ_j}不取均匀测度时, α(t)和ψ(t)的时间演化就会相互耦合而无法简单解耦, α(t)的动力学行为将会受到ψ(t)的影响, 此时系统的动力学在三维空间中进行, α(t)和ψ(t)的运动会变得非常复杂^[50,58]. 采用直角坐标, 设α = x+iy, f = Ref+iImf, 则三维方程(40)化为

$\begin{split}& \dot x = - uy + {\mathop{\rm Im}\nolimits} f(1 - {x^2} - {y^2}),\\& \dot y = ux + {\mathop{\rm Re}\nolimits} f(1 - {x^2} - {y^2}),\\& \dot \psi = u,\end{split}$

这里$u = 2 x\operatorname{Re} f + g - 2 y\operatorname{Im} f$.
图2给出了方程(60)的运动轨道在α平面的Poincare截面落点分布. 可以很清楚看到在|α|比较小的区域均为环面, α的演化为定态或周期振荡. 但在|α|较大的范围, 运动轨道是不规则、混沌的^[58]. 在环面区域, 类似于哈密顿系统, 环面对应于低维的运动, OA拟设可以成立. 在混沌区域, α(t)和ψ(t)的演化相互耦合, OA拟设不适用.
图 2 序参量α相空间的Poincare截面落点分布, 可以看到闭合环面和混沌散点
Figure2. The Poincare section of the order parameter α in phase space, where one can find the closed tori and chaotic scattered points.

2

-->

6.1.经典Kuramoto模型的同步相变

对于经典Kuramoto模型, 即平均场耦合振子系统, Kuramoto已成功地通过自洽方程理论得到了临界点及其临界行为. 下面利用OA拟设来研究该同步转变问题. 由于自然频率非均匀, 因此需用非均匀参数情况下的拟设. 利用序参量

$z\left( t \right) = {N^{ - 1}}\sum\limits_{j = 1}^N \exp \left( {{\rm{i}}{\theta _j}} \right) = \displaystyle\int\nolimits {\alpha _1}\left( {\omega, t} \right)g\left( \omega \right){\rm{d}}\omega $

定义式可知, 耦合函数中的分量分别为

${f_1} \equiv - Kz/(2{\rm{i)}},~~{\rm{ }}{\bar f_1} \equiv Kz/(2{\rm{i)}},~~{\rm{ }}{f_0} \equiv \omega,$

即Kuramoto模型耦合函数只包含到一阶傅里叶分量. 由(39)式可以得到

${\dot \alpha _1}(\omega,t) = [2{\rm{i}}\omega {\alpha _1}(\omega,t) + Kz(t) - K\bar z(t)\alpha _1^2(\omega,t)]/2, $

如果g(ω)为 Lorentz 分布

$g(\omega ) = 1/[{\text{π}}({\omega ^2} + 1)], $

由(63)式有

$z(t) = \int_{ - \infty }^\infty {\frac{{{\alpha _1}(\omega,t){\rm{d}}\omega }}{{{\text{π}}({\omega ^2} + 1)}}} .$

(64)式的积分可以将ω延拓到复平面的上半平面进行, 利用留数定理可以得到$z(t) = {\alpha _1}(\omega = i, t)$, 将其代入到(63)式和 (64)式中有

$\dot z = z(K - 2 - K|z{|^2})/2, $

此即序参量满足的动力学方程. 这个复方程可以写成振幅和幅角两个实方程. 可以看到, 该动力系统存在临界点K_c = 2, 当K ≤ K_c时, 系统有唯一不动点z ≡ 0, 此即非相干态; 当K ≥ K_c时, z ≡ 0解失稳, 系统分岔到新解

$|z| = \sqrt {(K - {K_{\rm{c}}})/K} . $

该解即同步解. 对比Kuramoto自洽方程得到的结果可见, 二者完全一致.
2 -->

6.2.奇异态研究

奇异态(chimera state)是近年来发现的一类对称性破缺导致的时空斑图涌现行为, 它描述了结构全同的单元(如相振子, 其振子的自然频率以及耦合方式都相同)在非局域耦合下会产生相干(coherence)和非相干(incoherence)共存的态^[59,60]. 这种对称性自发破缺的现象在很多生物系统中可以看到, 例如海豚与其他海洋哺乳动物、迁徙的候鸟等都具有一类有趣的半脑睡眠现象^[61], 即它们可以在左(右)脑休息的时候让右(左)脑保持清醒的状态. 人们在Kuramoto模型以及其他很多振子系统特别是神经系统中均发现了奇异态^[62-67], 并在光学混沌实验 ^[68]和化学Belousov-Zhabotinsky (BZ)反应实验^[69]实现.
结合相振子系统, 我们考虑如下的空间一维Ginzburg-Landau方程

$\frac{{\partial \phi (x,t)}}{{\partial t}} = \omega - \int_{{\text{π}}}^{\text{π}} {G(\Delta x)\sin (\Delta \phi + \alpha ){\rm{d}}x'}, $

这里?(x, t)为t时刻位于空间x处振子的相位, ω为振子自然频率, 设为空间均匀且ω = 0. 积分代表空间不同位置振子间的耦合, $\Delta x = x - x'$, $\Delta \phi = \phi \left( {x, t} \right) - \phi \left( {x', t} \right)$. α为振子间的相移(阻挫), 设0 < α < π/2. 非负耦合核函数G(x) ≥ 0反映空间不同位置振子之间的非局域相互作用, 通常设为归一化且随|x|增加而衰减的偶函数:

$\displaystyle\int\nolimits_0^\infty G\left( x \right){\rm{d}}x = 1, G\left( { - x} \right) = G\left( x \right).$

下面用OA拟设方法来研究系统(67)式的奇异态动力学^[62,70,71]. 设振子分布函数为f(x, ?, t), 它满足连续性方程

$\frac{{\partial f}}{{\partial t}} + \frac{\partial }{{\partial \phi }}(vf) = 0, $

其中相速度为

$\begin{split} & v(x,t) = \omega \,- \\ & \!\int {{\rm{d}}x'G(\Delta x)\int_{{\text{π}}}^{\text{π}} {{\rm{d}}\phi ' \sin (\Delta \phi + \alpha )f(x',\phi ',t)} }. \end{split}$

复序参量$Z(x, t)$可表示为

$\begin{split} \; & Z\equiv R(x,t)\exp \left[ {{\rm{i}}\Phi (x,t)} \right]\\ =\;& \int {{\rm{d}}x'G(x - x')\int_{{\text{π}}}^{\text{π}} {{\rm{d}}\phi \exp \left[ {{\rm{i}}\phi (x')} \right]f(x',\phi ',t)} } .\end{split}$

利用序参量表达式(70)式, (69)式的相速度$v(x, t)$可表示为

$\begin{split}v(x,t) =\; & \omega - \frac{1}{{2{\rm{i}}}}\left(\bar Z(x,t)\exp \left[ {{\rm{i}}(\phi + \alpha )} \right] \right.\\ & \left. - Z(x,t)\exp \left[{ - {\rm{i}}(\phi + \alpha )}\right] \right). \end{split}$

考虑到相位?的2π周期性, 可将分布函数表为相位?的傅里叶级数形式:

$\begin{split} f(x,\phi,t) =\; & \frac{1}{{2{\text{π}}}}\bigg\{1 + \sum\limits_{n = 1}^\infty \left[{h_n}(x,t)\exp ({\rm{i}}n\phi )\right.\bigg. \\ & \bigg.\left.+ {{\bar h}_n}(x,t)\exp ( - {\rm{i}}n\phi ) \right] \bigg\}, \end{split}$

利用OA拟设^[47], 系数h_n相互之间不独立, 它们满足幂律关系:

${h_n}(x,t) = {h^n}(x,t),\;{\bar h_n}(x,t) = {\bar h^n}(x,t), $

这里$h(x, t) = {h_1}(x, t)$. 将(73)式代入(72)式可得

$\begin{split} f(x,\phi,t) =\;& \frac{1}{{2{\text{π}}}}\bigg\{1 + \sum\limits_{n = 1}^\infty \big[{h^n}(x,t)\exp ({\rm{i}}n\phi ) \\ & + {\bar h^n}(x,t)\exp ( - {\rm{i}}n\phi )\big]\bigg\}. \end{split}$

将(74)式代入连续性方程(68)式, 并比较两边exp?(i?)的不同阶系数可得

$\frac{{\partial h}}{{\partial t}} = - {\rm{i}}\omega h + \frac{1}{2}[\bar Z\exp ({\rm{i}}\alpha ) - Z\exp ( - {\rm{i}}\alpha ){h^2}]. $

利用(74)式也可将序参量重新表达为

$Z(x,t) = \int {{\rm{d}}x'G(x - x')\bar h(x',t)} . $

利用(75)式和(76)式可得到序参量的演化动力学行为. 位于x处振子的相位长时分布函数为

$f(\phi )\sim \sum\limits_{n = - \infty }^\infty {{h^n}(x)\exp ({\rm{i}}n\phi )} .$

当$|h(x, t)| < 1$时, 对(77)式求和并做归一化可得

$f(\phi ) = \frac{{1 - |h{|^2}}}{{2{\text{π}}[1 - 2|h|\cos (\phi - \arg h) + |h{|^2}}}, $

此处arg代表复函数h的幅角. 上述分布函数为与(43)式形式相同的Poisson核函数, 其中心位于argh处, |h|刻画分布的非均匀性. 当|h|=0时, f(?)为均匀分布; 当0 < |h| < 1时, f(?)为单峰分布; 当|h|=1时, 由于分布(78)式分子分母均趋于零, 可求极限得f(?)简并为δ函数 f(?)=δ(? – argh), 它代表所有振子的锁相.
系数h(x, t)也可以看作是在空间x处exp(i?)的统计平均, 即

$h\left( {x, t} \right) =\int \nolimits_{ - {\text{π}}}^{\text{π}} {\rm{d}}\phi f\left( {x, \phi, t} \right){\rm{exp}}\left( {{\rm{i}}\phi } \right).$

与(70)式定义的序参量Z(x, t)相比较, h(x, t)没有对空间的积分, 即没有考虑非局域效应. 因此h(x, t)量度的是空间x处附近的同步相干性, 是一种局域序参量. 若对某一x处有|h(x, t)|=1, 则在$\left( {x - \varepsilon, x + \varepsilon } \right)\left( {\varepsilon \ll 1} \right)$范围内的振子会处于锁相(相干态), 而当|h(x, t)| < 1时, 则在$ \left( {x - \varepsilon, x + \varepsilon } \right) $范围内的振子处于非相干态, 以此可以给出两种态共存的奇异态.
在定态情况下, 进一步利用(76)式可以得到序参量R(x)的自洽方程, 进一步可以确定其空间分布. 当然, 也可用Kuramoto自洽方程理论来进行讨论, 可以得到与上述OA方法在定态时一致的结果. 进一步的计算需要借助于数值方法, 详细讨论可见我们的专著^[51]与综述文章^[59,60].
2 -->

6.3.高阶序参量及其动力学

很长时间以来, 人们聚焦于Kuramoto模型中的简单耦合函数情形开展研究. 在多数情况下, 振子间相互作用并不是简单的一阶正弦耦合函数. 对于一般的耦合函数, 在第4节中已经看到, 不同阶序参量之间是复杂的相互依赖关系. Daido^[72]早在20世纪90年代初就通过对一般耦合函数的Fourier展开提出了广义序参量, 但没有进一步讨论不同序参量之间的关系, 而是使用传统的Kuramoto自洽方程方法进一步讨论. Pikovsky和Rosenblum^[73]讨论了如下具有相同自然频率耦合振子的Kuramoto-Daido模型:

${\dot \theta _k} = \omega + \frac{1}{N}\sum\limits_{j = 1}^N {h({\theta _j} - {\theta _k})}, $

并引入了高阶序参量{α_n}. Pikovsky在该工作中早于Ott等^[47]提出了序参量之间的依赖关系(拟设条件), 即高阶序参量{α_n, |n| ≥ 2}均依赖于α_±1. 把(79)式中的耦合函数写为Fourier级数形式, 则动力学方程可写为

${\dot \theta _k} = \omega + \sum\nolimits_n {{h_n}{\alpha _n}\exp ( - {\rm{i}}n{\theta _k})}, $

其中{α_n}为(27)式定义的广义序参量. 考虑只有其中第l阶Fourier分量耦合的情形, 则动力学方程可写为

${\dot \theta _k} = \omega + \operatorname{Im} [H(t)\exp ( - {\rm{i}}l{\theta _k})], $

其中H(t)依赖于序参量α_l.
下面讨论利用WS变换来进行理论分析^[74]. 动力学方程(81)可以重新写为

$\frac{{{\rm{d}}[\exp ({\rm{i}}l{\theta _k})]}}{{{\rm{d}}t}} = {\rm{i}}l\exp ({\rm{i}}l{\theta _k})\omega + \frac{l}{2}[H - \bar H\exp ({\rm{i}}2l{\theta _k})]. $

注意到M?bius变换不唯一, 而其逆变换唯一, 则利用M?bius逆变换

$\exp ({\rm{i}}l{\theta _k}) = \frac{{\alpha (t) + \exp ({\rm{i}}{\varphi _k})}}{{1 + \bar \alpha (t)\exp ({\rm{i}}{\varphi _k})}}, $

(83)式可变换为{φ_k}, α及其时间导数$\left\{ {{{\dot \varphi }_k}} \right\}$, $\dot \alpha $:

$\begin{split}& \dot \alpha = l({\mathop{\rm i}\nolimits} \omega \alpha + H/2 - \bar H{\alpha ^2}/2),\\& {{\dot \varphi }_k} = l\{ \omega + {\mathop{\rm Im}\nolimits} [H(t)\bar \alpha ]\},\end{split}$

可以看到{φ_k}的方程右边不依赖于k, 说明变换后的相角{φ_k}均以相同角速度运动, 因此可引入新变量$\psi \left( t \right) = {\varphi _k}\left( t \right) - \varphi _k^0$, 其中{φ_k⁰}为一组运动常数. 显然$\dot \psi = {\dot \varphi _k}$, 说明M?bius变换映射可写为

${M_l}:\exp ({\rm{i}}l{\theta _k}) = \frac{{\alpha + \exp [{\rm{i}}(\varphi _k^0 + \psi )]}}{{1 + \bar \alpha \exp [{\rm{i}}(\varphi _k^0 + \psi )]}}, $

与一阶耦合函数相比, 这里形式上只是多出了l因子.
很多具体问题可以应用上述结果加以讨论. 例如近年来得到关注的三体相互作用平均场振子系统^[75,76]

${\dot \theta _k} = \omega + \frac{1}{{{N^2}}}\sum\limits_{i,j = 1}^N {\sin ({\theta _i} + {\theta _j} - 2{\theta _k})},$

该系统对应于方程(81)中l = 2, H(t)=α₁²的情形, 因而可利用序参量方程(84)进行讨论.
再例如具有阻挫的二阶简谐耦合振子系统^[58]

${\dot \theta _k} = \omega + \frac{1}{N}\sum\limits_{j = 1}^N {\sin [2({\theta _j} - {\theta _k} + \gamma )]}, $

其中γ为相移因子, 代表阻挫效应. 对于这种情况, $l = 2, \;H\left( t \right) = {\rm{exp}}\left( {2{\rm{i}}\gamma } \right){\alpha _2}$. 二阶或高阶简谐耦合函数在很多系统中都存在, 因此近年来的研究很多^[58,77-79], 也有相应的一些实验研究^[80-83], 这里不再展开.

综上所述, 大自由度复杂系统的同步是一类典型、基本的涌现行为. 这种涌现行为的出现是大量个体或自由度共同参与、通过相互作用自组织形成的整体有序. 有序意味着复杂系统低维宏观行为的出现, 它往往是以序参量的出现作为标志的, 可借助于统计物理学思想和方法来加以研究. 就序参量的涌现机制来看, 协同学的理论与方法揭示了大量自由度的复杂系统如何通过内部自组织与竞争产生出有序, 其中非常重要的理论任务就是对复杂系统通过集体变量的分析来甄别快慢变量, 并利用支配原理来进行降维, 得到序参量的动力学. 降维方法的精髓在于尺度的可分离性, 特别是时间尺度的快慢可区分性, 这种特性可以其他形式体现出来, 如拓扑流形中中心流形的出现、变量中守恒律的存在等等.
在复杂系统同步的序参量动力学理论方面, 除了Kuramoto的建立和求解序参量满足的自洽方程方法以外, Ott和Antonsen提出的拟设是近年来应用较为广泛的方案, 该方案首先引入集体坐标即各阶广义序参量, 然后通过各阶量之间的非独立性来达到降维的目标. 这里面就隐含了绝热消去快变量的思想. 需要指出的是, OA拟设在很多情况下是一种推测, 且在很多情况下不只是一阶序参量是慢变量, 另外的高阶序参量也可能是慢变量, 这需要从原始的各阶序参量方程出发, 利用协同学基本原理来进行分析.
下面讨论一下对于实际系统同步作为整体动力学的研究分析. 利用序参量动力学对同步的研究既可以是理论建立模型的研究, 也可以是对实际问题利用的降维方法进行研究. 相比于基于模型的理论研究, 对于实际复杂系统有序行为的序参量动力学重构和研究显然有着更为重大的意义和价值. 对于很多实际中的复杂动态系统, 在宏观层面构建其序参量动力学是可能的, 但也富于挑战性. 序参量动力学对实际问题的研究是通过采集微观数据, 然后在此基础上进行宏观动力学重构, 因而是一个跨尺度动力学重构问题, 迄今为止这仍然是一个开放的领域, 与复杂网络重构研究具有同样重大的价值, 是对复杂系统动力学方面的重要探索范式, 值得大力开展研究. 就本文讨论的相位动力学及其同步而言, 其目前在应用方面的主战场之一是神经与脑科学^[84]. 我们不妨以此为例进行说明. 神经元的主体微观动力学是电信号的积累发放, 目前实验方面已经可以采集大量实时数据并进行分析整合. 对于脑电波及其各种动力学协同行为的研究, 可以先抽取不同部位采样的脉冲时间序列, 再定义对应的相位, 就可以建立相位的时间序列信息, 据此计算不同序参量, 通过分析序参量的低维行为进行讨论. 就目前的研究来看, 很重要的一个课题是序参量动力学模型的重构, 这一问题至今无论在神经科学还是其他领域都还没有很大的突破, 但已有一些具有启发性的研究思路. 例如, Zhang等^[85]和Chen等^[86] 近几年提出了主超前相位方法、动力学权重因子方法以及网络动力学重构方法, 以此可以分析复杂系统的动力学并结合序参量层次的动力学描述, 完全有可能提供一个可行的图景. 这些工作还需要大力推进.