活性物质^[1]是由自驱动单元组成的非平衡态体系, 这些单元通过持续消耗能量产生定向运动、力学响应或形变, 从而涌现出团簇、涡旋、条带等丰富的集体行为^[2,3]. 活性物质在自然界中广泛存在, 从微观尺度的分子马达与细菌悬浮液; 到宏观尺度的鸟群和鱼群, 均是典型代表.

在复杂环境中, 活性物质不可避免地与边界及障碍物相互作用, 这会显著改变其输运特性与集体模式. 例如免疫细胞在致密组织中的迁移^[4], 以及细胞质内微丝微管的自组织^[5]. 大量研究采用障碍物点阵或单个大型障碍物来模拟此类复杂环境^[6–9], 结果表明, 障碍物可以作为成核中心诱导粒子聚集, 触发运动诱导的相分离(motility-induced phase separation, MIPS)^[7,8], 并在非对称扭矩驱动下形成旋转团簇. 类似现象也在活性聚合物系统中被观察到^[10]. 此外, 障碍物的几何形状、密度与空间分布, 以及粒子的持久性和手性等内在属性, 都会显著影响自驱动粒子的输运行为, 使其呈现堵塞、快速迁移或定向输运等不同状态^[4,11,12].

值得注意的是, 传统模型通常将活性粒子简化为形状和力学性质不变的刚性球体或棒状^[13–16]. 这一简化极大地推动了ABP(active Brownian particle)、Vicsek等经典理论的发展^[12], 但也忽略了一个关键事实: 真实生物体系中的许多个体具有柔性与可变形性. 例如迁移中的真核细胞^[17]、螺旋体细菌^[5]以及细胞骨架驱动的囊泡^[18], 均能在运动过程中发生显著形变. 近年来的研究表明, 柔软性能够显著改变活性物质的相行为, 促进动态团簇^[19]和拓扑缠结网络的形成^[20], 并与内部活性物质分布耦合, 诱导非对称形变^[21]、旋转^[22]及定向迁移^[23].

在受限环境中, 活性粒子的柔软性同样发挥着关键作用. 当柔性活性粒子穿越几何约束(如窄通道、孔隙)时, 其形变能力使其能够通过比自身尺寸更小的开口, 呈现出与刚性粒子截然不同的输运行为^[19,22]. 例如, 可变形微生物能够通过“爬行”或“蠕动”模式在复杂环境中高效迁移^[24–26]; 柔性活性链则可以在障碍物周围形成动态团簇^[10], 表现出反常扩散^[27]和定向输运^[28]. 这些现象表明, 柔软性不仅是活性粒子的被动特性, 更是调控其运动行为的关键因素.

先前工作中发现, 在没有对称性破缺因素的前提下, 当活性足够高且障碍物尺寸适当时, 球形自驱动粒子可在圆形障碍物周围形成持续涡旋^[9]. 而热噪声影响粒子的堆积与重排, 存在一个最佳噪声可最大化涡旋速度与极性取向^[29]. 本文在上述基础上, 将刚性粒子替换为可形变粒子, 重点探讨粒子的软硬程度如何影响涡旋形成. 研究结果表明, 粒子柔软度是调控障碍物附近涡旋行为的关键参数. 随着粒子柔软度变化, 体系从随机态经过渡态转变为涡旋态. 在该体系内, 受极性选择性捕获与粒子形变能力之间的竞争的影响, 可变形粒子的涡旋态内部存在一个最佳柔软度, 使最内层粒子的平均涡旋速度与极性切向分量达到峰值, 涡旋对柔软度呈现非单调依赖. 结构分析进一步显示, 障碍物附近粒子呈现明显的分层结构. 同时, 六重取向序参量和边界厚度涨落随柔软度单调变化, 表明最佳涡旋流动并非对应于最有序的结构, 而是出现在团簇极性有序与粒子柔软度相平衡的中间区域. 本文研究结果不仅为理解可形变活性物质在复杂环境中自组织行为提供新的依据, 也有望为活性材料在靶向药物输送^[30]、微流控器件^[31]及软体机器人等领域的应用提供参考.

在二维模型中, 整个体系由半径为$ R $的圆形障碍物与质量为$ m $、尺寸为$ \sigma $的小的圆形自驱动粒子(self-propelled particles, SPPs)组成, 处于同一模拟盒子中. 障碍物的中心固定在模拟盒子的中点, 设为原点. 障碍物是由一系列尺寸为$ \sigma $的珠子模拟, 这些珠子紧密连接, 为使障碍物表面光滑, 相邻珠子之间的间隔约为$ 0.2\sigma $. 所有的圆形自驱动粒子与边界的珠子都是通过纯排斥的Weeks-Chandler-Anderson (WCA)势相互作用, 表达式为

其中$ r $为粒子与珠子之间的距离, $ \varepsilon $为相互作用强度, $ {r}_{\text{cut}}={2}^{1/6}\sigma $. 不同的是, 圆形自驱动粒子之间采用伸缩势, 表达式为

其中$ {r}_{ij} $为自驱动粒子之间的距离; $ k $为自驱动粒子之间的相互作用强度, 其可以表示为粒子的柔软度; $ {r}_{0}=\sigma $, 这使得自驱动粒子具有一定柔性, 在相互作用时能够发生形变(图1(a)). 第$ i $个自驱动粒子的运动遵循朗之万方程:

方程(3)控制粒子的平动运动, $ {r}_{i} $是粒子$ i $的位置矢量, 第1项表示粒子间的相互作用力, 第2项为摩擦力, 第3项为活性力, 第4项为随机力, 而方程(4)则决定与推进方向$ {\hat{\mu }}_{i}=(\cos {\theta }_{i}, \sin {\theta }_{i}) $耦合的扩散转动动力学. $ {U}_{i} $是粒子$ i $的总势能, $ \gamma $和$ {D}_{\text{r}} $分别为平动摩擦系数和旋转扩散系数, $ {F}_{\text{a}} $为每个粒子$ i $上的推进强度, $ {\eta }_{{{}_{i}}} $和$ {\xi }_{i} $为单位方差高斯白噪声.

使用自主修改的LAMMPS软件进行模拟, 模拟盒子大小为$ 200\sigma \times 200\sigma $, 采用周期性边界条件. 为了减轻高持久性下的有限尺寸效应, 当粒子穿过模拟盒子边界时, 粒子的推进方向将随机化^[32](图1(b)). 黑色方框表示模拟盒子, 红色虚线选定的黄色区域为随机化粒子自推进方向的区域, 其中带有蓝色箭头的红色实心粒子表示进入该区域前的推进方向, 带有绿色箭头的红色空心粒子表示进入该区域随机化后可能的推进方向. 本文采用约化单位, 设置$ m=1 $, $ \sigma =1 $, $ \varepsilon =1 $, 因此时间单位$ \tau =\sqrt{m{\sigma }^{2}/\varepsilon } $. 这里保留了方程(3)中的惯性项, 以便可以使用velocity-Verlet算法和较大的时间步长. 设置摩擦系数$ \gamma =10 $, 使得自驱动粒子的运动基本上是过阻尼的. 在模拟中, 固定$ {F}_{\text{a}}=5 $, 即漂移速度$ {v}_{\text{a}}={F}_{\text{a}}/\gamma =0.5 $, 并且设置$ {D}_{\text{r}}=0.001 $以实现高持久性. 先前的研究发现^[29], 低噪声或者无噪声时, 障碍物容易诱导涡旋形成, 而较大的噪声对涡旋形成具有抑制效果, 并且模拟了低噪声下柔软度对涡旋行为的影响, 与$ {k}_{\text{B}}T=0 $时的结果定性一致. 为此本文设置$ {k}_{\text{B}}T=0 $忽略热噪声的影响, 专注于研究不同柔软度下粒子的自驱动行为. 体系中粒子浓度保持不变, $ {\phi }_{}=\rho \text{π} {\sigma }^{2}/4=0.15 $($ \rho $为自由空间中自驱动粒子的数密度), 该浓度较低, 在体相中不会自发形成大团簇或相分离. 所有统计数据至少对6个独立模拟结果进行平均, 每次计算得到的结果运行持续15000τ.

当粒子活性较高且障碍物尺寸处于适当范围内时, 自驱动粒子倾向于在障碍物附近积聚, 并形成持续的涡旋运动, 图2(a)给出了长时间平均的粒子流场$ \overline{\displaystyle\sum \nolimits_{i}{\boldsymbol{v}}_{i}\delta ({\boldsymbol{r}}_{i}-\boldsymbol{r})} $($ \overline{\cdot \cdot \cdot } $表示长时间平均). 如果障碍物的尺寸过大或过小, 或者粒子的持久性($ 1/{D}_{\text{r}} $)不足, 则障碍物附近的粒子将经历从动态有序的涡旋态(vortex state, VS)到过渡态(transition state, TS)或随机态(random state, RS)的转变. 本文通过调节自驱动粒子之间的相互作用强度$ k $, 也能诱导障碍物附近粒子在上述3种状态之间转变(图2(b)).

图3是$ R=35 $时, 3种状态下障碍物附近最内层的粒子的角速度$ \omega $和累积转过角度$ \Delta \theta $(涡旋在长时间累积下转过的角度和)的时间演化图像. 当相互作用强度$ k $较小时, 体系处于RS态, 粒子整体处于无序状态, $ \omega (t) $和$ \Delta \theta (t) $在零附近波动. 当相互作用强度$ k $处于适当的大小时, 体系进入TS态, 粒子可沿着某一个方向绕障碍物持续运动一段时间, 但在受到大涨落作用时, 绕行方向发生反转, $ \omega (t) $和$ \Delta \theta (t) $均表现出围绕零值的大幅度的波动与反转运动. 当相互作用强度$ k $进一步增大至较大值时, 体系转变为VS态, 粒子在障碍物表面形成稳定的定向旋转运动, $ \omega (t) $围绕某一非零值进行小幅波动, 而$ \Delta \theta (t) $则随时间呈现出单调递增(逆时针旋转)或者降低(顺时针旋转)的趋势. 也就是说, 随着活性粒子逐渐接近于刚性粒子, 粒子将在障碍物表面形成定向旋转运动.

尽管通过调节粒子的柔软度可以使体系从无序的RS态向有序的VS态进行单向转变, 但有趣的是, 与噪声对涡旋行为的影响类似, 在VS态内部, 柔软度对涡旋结构的影响也呈现出非单调特性. 图4(a)展现了在不同柔软度条件下, 障碍物表面最内层粒子的平均涡旋速度$ {\overline{v}}_{1}= \left\langle \dfrac{1}{N}\left| \displaystyle\sum\nolimits_{i}\boldsymbol{r}_{i}\times {\boldsymbol{v}}_{i}/r_{i}^{2}\right| \right\rangle $, 其中$ N $为障碍物最内层粒子的数目, r_i表示最内层粒子$ i $沿着径向方向的单位矢量. 结果表明, 对不同尺寸的障碍物, 均存在一个最佳柔软度, 使得最内层的涡旋速度达到峰值, 并且随着障碍物半径增大, 这个最佳柔软度相应地向更大值偏移. 涡旋的形成主要是源于粒子的取向整流, 为了验证这一观点, 我们进一步计算了障碍物表面最内层粒子极性的平均切向分量$ {\overline{\mu }}_{1}=\left\langle \dfrac{1}{N}\left| \displaystyle\sum\nolimits_{i}\boldsymbol{r}_{i} \times {\boldsymbol{\mu}}_{i}\right| \right\rangle $ (图4(b)). 结果显示, 极性切向分量与平均涡旋速度的变化趋势一致, 在不同柔软度下都表现出非单调依赖关系. 值得注意的是, 障碍物半径较小时, $ {\overline{v}}_{1} $对应的最佳柔软度$ k $比$ {\overline{\mu }}_{1} $对应的值要小. 这一偏差主要是取向有序性与粒子柔软度竞争的结果. 半径较小时, 障碍物具有较大的曲率, 迫使粒子在运动到障碍物表面时频繁改变运动方向. 当粒子较软时, 粒子间相互作用容易诱导发生形变, 削弱取向的协同性; 当粒子较硬时, 尽管取向稳定性增强, 但粒子之间容易相互阻塞, 影响整体旋转运动. 因此, $ \overline{v} $和$ \overline{\mu } $的两者峰值对应的最佳柔软度存在一定偏移.

与此同时, 模拟中粒子在障碍物表面聚集时会呈现明显的分层结构. 为了分析不同柔软度下层间粒子的动态特性, 计算了障碍物表面粒子的平均涡旋速度$ \overline{v} $和极性的平均切向分量$ \overline{\mu } $随径向距离的变化关系(图5), 其中$ r-R $为距离障碍物表面的径向距离. 结果表明, 在靠近障碍物表面的位置, $ \overline{v} $和$ \overline{\mu } $都呈现线性衰减趋势; 但距离障碍物表面较远区域则呈现指数衰减(见黑色虚线拟合曲线), 两者之间由一个小的肩状区域过渡连接. 在靠近障碍物表面, $ \overline{v} $线性衰减是由于粒子在这个区域累积形成团簇. 而在团簇区域之外, 粒子表现为稀疏的气相行为, 粒子难以维持聚集状态, 取向相对随机化, 因此$ \overline{v} $和$ \overline{\mu } $都随距离表现为指数衰减. 而中间出现肩状特征, 是由于粒子在不断接近障碍物表面时, 原本在障碍物表面的粒子被不断挤压, 粒子被推动加速, 导致最里层粒子的$ \overline{v} $和$ \overline{\mu } $较大. 而此区域外面几层, 由于粒子在交换与重排, 粒子径向运动较强而切向运动较弱. 与之相比远场区, 粒子重排效应减弱, 切向运动有所恢复, $ \overline{v} $和$ \overline{\mu } $值有所上升. 进一步观察发现, 第1层的$ \overline{v} $和$ \overline{\mu } $与图4结果一致, 呈现非单调性. 然而, 线性衰减区域的斜率随粒子的柔软度变化而单调变化的: 粒子的$ k $值越大, 速度随时间下降越平缓, 且线性范围越广. 这说明粒子越硬, 层间协同性越强, 更有利于维持稳定的涡旋运动.

为了分析不同柔软度下, 涡旋形成后能够稳定存在的原因, 统计了进入距离障碍物表面一定范围内, 两种极性取向粒子对应的概率. 其中, 切向极性分量与涡旋旋转方向一致的粒子被定义为pSPPs, 相反取向的粒子则被定义为nSPPs. 相应地, 不同极性取向的粒子进入距离障碍物表面$ r-R $范围内的概率, 则分别表示为P(pSPPs)和P(nSPPs). 结果表明, 越靠近障碍物表面, 体系对pSPPs的接受偏好越显著(图6(a)). 在障碍物附近, 随着柔软度改变, P(pSPPs)和P(nSPPs)呈现单调变化趋势, 粒子的$ k $值越大, 涡旋结构对pSPPs的选择偏向性越明显. 这种现象主要是由于软粒子具有更强的可形变能力和层间滑移特性, 在缓慢的碰撞与挤压过程中, 粒子朝向可以不断调整并相互适应, 即使是nSPPs也较容易进入最内层区域. 相较之下, 硬粒子形成的堆积结构更加紧密, pSPPs一旦进入团簇, 其推进方向不容易发生改变, 而nSPPs容易在挤压过程中被弹出. 从最内层的P(pSPPs)统计结果, 如图6(c)也能明显看出, $ k $值越大, 障碍物表面的团簇对pSPPs的选择偏向性越大.

尽管不同柔软度下粒子对pSPPs均表现出一定的极性选择偏向性, 并且这种偏好对涡旋的形成与维持产生一定影响, 但由于不同$ k $值下粒子发生形变的难易程度存在差异, 粒子在运动过程中容易重新取向, 因此, 障碍物表面团簇内粒子的极性分布特征能够更直接决定涡旋的性质. 图6(b)给出了在障碍物表面$ r-R $范围内两种粒子的比例Q(pSPPs)和Q(nSPPs). 在障碍物表面最内层明显发现, Q(pSPPs)和Q(nSPPs)随粒子的柔软度变化表现出非单调性关系, 图6(d)进一步计算了障碍物最内层的Q(pSPP)关于柔软度的变化, 更加直观地体现了这一非单调行为. 这种非单调关系的峰值位置与图5(a), (b)中$ {\overline{v}}_{1} $和$ {\overline{\mu }}_{1} $的变化是相对应的.

不同柔软度的粒子在障碍物表面聚集时, 会形成高密度的团簇区域, 并且呈现明显的层状结构(图7). 为了进一步探究不同柔软度的粒子所形成涡旋的结构与特征, 计算了障碍物表面粒子的径向分布函数(图8(a)). 当粒子较软时($ k $值较小), 容易发生形变, 导致层状结构较为模糊. 随着粒子逐渐变硬, 障碍物表面分层结构越来越清晰, 堆积厚度显著增大. 同时密集堆积的结构也增大了团簇中粒子排列的有序性. 为此, 计算了障碍物表面团簇内粒子的六重键取向序参量$ {{\varPsi }}_{6} $, 表达式如下:

其中$ j $为团簇内粒子, $ {q}_{6}(j)=\displaystyle\sum \nolimits_{k=1}^{{Z}_{{{}_{j}}}}{\text{e}}^{\text{i}6{{\theta }_{jk}}}/{Z}_{j} $是第$ j $个粒子的取向参量, $ {Z}_{j} $是与第$ j $个粒子相邻的粒子的数目, $ {\theta }_{jk} $是第$ j $个粒子与其第$ k $个相邻粒子之间连线形成的矢量与参考方向x轴所构成的夹角. 从图8(b)可以看出, 随着$ k $值增大, 粒子在障碍物边界的聚集程度提高, 堆积结构变得更加有序.

与此同时, 计算了聚集在障碍物表面的团簇厚度的角向涨落$ \Delta h $, 其表达式如下:

其中$ h(\theta ) $表示方位角$ \theta $处的团簇厚度, $ \overline{h} $表示平均厚度. 如图8(c)所示, 在粒子较硬的情况下, 团簇外围粒子与团簇内部粒子进行交换与重排. 然而, 由于粒子硬度较大, 形变能力较弱, 难以像软粒子一样通过形变调整自身取向, 因此这种重排过程往往会伴随着团簇的局部破裂, 从而导致角向涨落$ \Delta h $出现大幅度变化. 随着$ k $值增加, $ \Delta h $呈现单调增长趋势.

关于活性物质与障碍物之间的研究, 一直是一个值得探索的方向. 本文研究了具有可调柔软度的可形变自驱动粒子在光滑的圆形障碍物表面的集体行为. 结果表明, 随着粒子的硬度增大, 体系会发生从随机态(RS)到过渡态(TS)最后表现为涡旋态(VS)的转变. 进一步分析发现, 涡旋态内部的动力学特性对粒子的柔软度呈现非单调依赖: 存在一个最佳$ k $值使得涡旋最内层的涡旋速度与极性取向达到峰值, 这主要是因为障碍物表面团簇对不同取向粒子的偏向性捕获. 值得注意的是, 涡旋的形成不仅受粒子极性取向的影响, 还与粒子柔软度有关, 当粒子较软时取向容易发生改变, 二者之间相互竞争导致涡旋速度与极性取向峰值出现偏移. 与此同时, 结构特征分析表明, 随着粒子变硬, 障碍物表面的堆积结构有序度单调性变化. 这也再次表明, 最有效的涡旋的建立在团簇极性取向与粒子柔软度的微妙平衡上. 总之, 本文的研究不仅增进了对可形变活性粒子在复杂介质下的集体运动的认识, 也有望为活性物质在生物、医学、材料等领域的设计与应用提供参考.