在本节中，我们解释了我们如何拍摄密集的人群以及如何将原始图像转换为定量数据。

　　我们在潘普洛纳（Pamplona）拍摄的人群由圣费尔米（SanFermín）节的与会者组成。他们中的大多数人都是当地人，他们非常了解活动的发展方式：他们知道着装要求以及何时举起红色的手帕，在节日开幕式演讲和烟火之后，他们热切地等待乐团。许多外国游客也存在，人群主要由所有性别的年轻人组成。由于Chupinazo是在西班牙国家网络上播出的，并在国际上广受欢迎，因此我们可以假设大多数参与者都知道他们在聚集之前将面临的高人群密度。结果，Chupinazo的参与者充满了喜庆的心情，并有明显的饮酒迹象。当人群稀疏时，可以看到一小群人唱歌，鼓掌和轻轻地推动彼此而不会产生重大的动作。但是，当人群密度超过ρ*时，我们会观察到高振幅运动，在不同的上下文中被称为人群地震21。例如，在重金属音乐会上，没有人在人群茂密时积极推动他人的迹象。相反，人们试图保护自己免受人群的侵害。最后，那些无法进入入口广场的人继续访问广场并向内推动广场。两个入口被警察部队阻止，直到中午前几分钟。在这里，数百人（主要是当地人）尝试聚集在已经拥挤的广场上。广场上的人数在中午占据了高峰，可能大约有6,000人。

　　我们在2019年，2022年，2023年和2024年拍摄了四场Chupinazo活动。2020年和2021年的版本由于190限制而被取消。天气在2019年和2023年是晴天。在2022年和2024年，天气多云和多雨。我们没有被告知，也没有观察到仪式组织的任何变化。

　　我们从图1B中指示的两个观察点中录制了视频。在2022年，2023年和2024年，我们在两个景点上使用了Sony FDR-AX33和Sony FDR-AX43A摄像机。帧速率为每秒25帧，分辨率为3,840×2,160像素（4K视频）。2019年，我们使用了安装在尼康D500摄像机上的远摄镜头（Sigma 150-600 mm F5-6.3 DG OS HSM S），帧速率为30 fps，分辨率为3,840×2,160像素。We illustrate the observation set-up in Extended Data Fig. 6. The cameras were mounted on standard tripods on the balcony of the leftmost building (fifth floor of the tourism office of Pamplona,​​ 2023 and 2024), and on custom mounts on the rightmost observation point (private balcony on the fourth floor of the highest building surrounding the Plaza Consistorial, 2019, 2022, 2023 and 2024).

　　我们使用了与参考文献相同的协议。12纠正透视扭曲。简而言之，我们将每个帧的像素坐标（u，v）映射到方形框架（x，y，z）中的坐标。目标是生成具有固定像素尺寸的图像。这是使用平面同构图完成的，假设针孔摄像机并忽略了Z轴的高度。为此，我们首先选择了广场的参考图像（扩展数据图7a），当远离行人时。我们选择了在广场上分布的四个点，并在Google Earth Pro的卫星视图上与它们匹配（扩展数据图7b）。然后，我们使用两个视图中的这四个点计算了平面同构矩阵，并相应地重新缩放了所有图像。最后，通过以米和像素为单位的卫星视图上的距离，我们确定了以米为单位的像素的大小。在重新制定的图像中，所有像素具有相同的维度（扩展数据图7C）。在图7的扩展数据中，我们在原始图像，卫星视图和所得重新制定的图像上显示了2023年使用的点集。校正图像上的像素几乎是正方形，它们的尺寸分别在x，y方向对应于0.028×0.029 m，0.046×0.045 m，0.012×0.012×0.012 m和0.012 m和0.013×0.012 m，分别在2019，2022，2022，2023，2023和2024和2024中。

　　为了估计平面同型协议引起的相对误差，我们将原始参考点之间的距离与校正图像中测得的距离进行了比较。我们在x，y方向上分别发现了以下相对错误：在2019年，2022年，2023年，2023年和2024年分别分别为1％×12％，1％×1％，1％×3％和1％×3％。

　　如参考。12，我们还校正了零高近似值。以前，我们假设图像中的所有可见像素都位于与地面同一平面上。但是，我们观察到不在地面上的行人的头。在这里，我们回想起该公式，并给出为零高近似校正的值测量的值。可以在参考文献中找到更多详细信息。12。令ΔX'表示在校正的图像中测量的两个头之间测量的距离。由于行人的高度，该距离被高估了。令ΔX表示同一图像中同一行人的脚之间的距离。这些距离与以下公式相关

　　其中H是观察点的高度，H是行人的典型高度。这种关系在各个方向上都是有效的，并且独立于图像中的位置。我们根据西班牙人的平均身高设置H = 1.73 m，并通过使用带有行人用作长度参考的行人的建筑物图片估计H = 16±1 m。

　　最后，在2019年，天气晴天，建筑物在我们的粒子图像速度法（PIV）测量中引入了伪像的阴影（图1A）。为了解决此问题，我们执行了一个附加的图像处理步骤。我们使用了Github上可用的Python脚本来匀浆Image32,33的亮度。此过程删除了建筑物阴影边缘出现的所有PIV伪影。

　　每年，我们专注于大约1小时的时间间隔。我们将时间的起源定义为节日参加者抬起红色手帕以表示节日的开放，并根据Chupinazo的传统改变其行为。这发生在中午左右，在音乐节正式开始之前，没有外部指导决定了人群的动态。没有声学，视觉或机械信号可以协调人群的运动。

　　节日开幕之前，庆祝活动的一部分包括握着大国旗并在广场上玩大气球（扩展数据图7d）。当我们测量密度和速度场时，我们需要确保仅考虑没有阻碍我们视野的这些物体的空间区域。

　　为此，我们使用了标准的机器学习工具Yolov834。它旨在检测图像上的指定对象。从2019年和2022年开始，我们通过手动注释多边形，将其封闭在我们的图片上的气球和旗帜上，从而在2019年和2022年对其进行了训练。实际上，我们使用原始图像上的LabelStudio软件注释了图像。然后使用平面同构参数校正检测到的对象的面积和位置，并将其转换为应用于每个帧的动力掩模（扩展数据图7E）。密度和速度场均在掩盖图像上计算（扩展数据图7F）。

　　为了检测人群在人群中的头寸，我们利用了机器学习算法P2PNET35的优势。它旨在根据图片中的卷积神经网络估算图片中的人数，这些神经网络在300张图像上备有密集的人群（Shanghaitech DataSet36）。我们在2019年，2022年和2023年录像中拍摄的4张手动注释图像中进一步训练了神经网络。为了注释图像，我们使用了一个自制的Python脚本和经过训练的P2PNET，用于3,500个4批次的时代。图片的整个视野​​中的自动计数程序和手动计数过程之间的差异小于5％。我们检测到原始图像上的头部（图1C，D）。

　　我们的目的是通过最大程度地减少误报和遗漏检测的影响来清理数据。鉴于行人不断出现在框架中，与误报或错过的发生相比，他们更可能被正确检测到。我们预测检测错误表现为高频“闪光”。为了抵消这一点，我们通过在检测点上应用高斯卷积，将检测数据转换为密度场，其中3σ大约等于头的大小，σ是二维高斯的典型宽度。然后将该密度场平均在0.4 s以上，并应用33％的阈值以仅保留高点的点。随后采用了最大值找到算法来获取清洁的检测数据。

　　为了定义密度字段，我们使用了相同的窗口（尺寸和位置）与速度场（补充表1）。在每个窗口中，大约为1.5×1.5 m2，我们计算了人数，并将该值划分为窗口区域。鉴于图像的质量，我们能够检测到个人并在图1B中仅针对2022年，2023年和2024年图像定义了局部密度场。

　　我们使用Trackpy37测量图3D所示的轨迹。

　　我们将Chupinazo人群描述为Continua。在流体力学中经常使用，我们使用常规PIV12测量了人群中的Eulerian速度场。我们使用了约1.5 m的典型窗口大小，典型的时间步长约为0.5 s。补充表1总结了表征观察到的动力学的精确不同值。单个PIV盒中速度场的两个组件的测量误差在2019年和2022年为0.1 m s -1，在2023年和2024年为0.05 m s -1。

　　在本节中，我们解释了数据的起源以及将原始图像转换为定量数据的过程。

　　2010年的Love Parade是一个在7月在德国Duisburg举行的音乐节。在活动期间，在节日区的主要入口处，有21人死亡，500多人受伤（图8扩展数据）。这一创伤事件已经成为众多科学研究的主题，也是因为许多视频显示了与灾难的融合，已公开发布4。

　　互联网上有许多手机和相机的录音。组织者发布了一些灾难发生的主要条目的摄像机镜头（https://www.youtube.com/watch?v=qpzisdbe3da）。我们选择了两个剪辑（从14:04到15:46剪辑1，从19:24到19:49剪辑2），在那里我们观察到人群的高振幅运动。

　　我们无法对2010年爱情游行的数据进行任何图像校正。

　　由于低分辨率，我们无法测量夹子1中的密度，这阻碍了检测和手动注释。在剪辑2中，我们注释了一个图像。我们通过用像素来测量节日参与者的肩膀宽度，然后将其转换为米，估计了视野覆盖的区域。我们对人群密度的估计值为ρ= 8±1 m -2。

　　我们为Chupinazo和Love Parade人群使用了相同的方法。我们使用常规PIV12测量了人群中的Eulerian速度场。由于图像上没有足够的参考点可以执行透视校正，因此我们使用原始图像来测量近似速度场。我们总结了所选的不同参数的值，以表征补充表1中观察到的动力学。

　　我们使用与Chupinazo人群相同的方法。

　　我们主要使用Python数值语言的Numpy软件包开发了所有数值工具。

　　我们使用了来自参与者检测的数据（请参见上文）和static.pair_correlation_2d从Trackpy37考虑到非周期性边界条件。径向对相关函数是在没有标志和气球的矩形中计算的，并在最密集的人群中平均约7分钟（见下文）。

　　为了计算速度V（R，T）的功率谱，速度取向以及图3A – C的平方速度V2（R，T）字段，我们首先使用Numpy.fft.fft函数计算了每个位置R的离散时间傅立叶变换。为了一致性，我们通过信号的总持续时间将这些值归一化。具有脉动ω0的完美余弦信号将具有离散的傅立叶变换，显示振幅为±ω0的幅度为1/2的峰值。局部光谱由每个离散时间傅立叶变换信号的平方标准给出。最后，我们进行了局部光谱的空间平均值来计算SV（ω）和。这些光谱是使用2023数据的持续时间30 s的信号进行计算的，2022数据为7分钟30 s，2019年数据为30 s。在2019年，由于过热，我们的相机的安全系统将它们重置大约每30 s。因此，我们在较短的时间间隔内计算了光谱，但是在我们所有30米长的收购中平均它们。

　　为了准确比较2019、2022和2023光谱，我们将光谱除以它们的平均值在ω= -2 rad s -1和ω= 2 rad s -1之间，并在5个后续数据点上进行最终平滑后获得了最终光谱。

　　为了对图3F进行时间与频率分析，我们应用了上述方法，以时间间隔[T，T+ΔT]使用ΔT= 3分钟计算功率谱，并重复相同的操作以增加t的值，以1 s的步骤增加。

　　为了计算速度，方向和速度场的径向相关函数，我们首先计算了它们的二维空间相关函数。向量函数f（r，t）的相关cf（r）的值由

　　其中Δf（r，t）= f（r，t）-f（r，t）r，t，r，t，t平均为r和时间t的平均位置。实际上，我们使用numpy在点r + r上计算ΔF的值。滚动以移动形成向量ΔF（r）的数据。然后，我们使用numpy.nanmean来计算平均值。我们使用多处理。pool并行了此计算。为了计算径向相关函数，我们创建了与径向位移相对应的垃圾箱。对于从R0到R0+ΔR不等的bin，我们平均R0≤r的二维相关函数C（R）的所有值< r0 + δr. We defined the instantaneous correlations by performing the time average over 3-min-long intervals. We defined the correlation lengths as the distances at which the radial correlation functions reach the value 0.1.

　　Here we focus on the oscillatory component of the crowd motion. We first applied a temporal band-pass filter to the velocity field. The filter width was set to the full-width at three-quarter maximum height of the power spectra shown in Fig. 3a. In practice, we performed a time Fourier transform on the velocity field, we applied the band-pass filter and performed an inverse Fourier transform.

　　We defined the local spin field ϵ(r, t) as the sign of the instantaneous increment of the angle θ made by the velocity with the x axis. For each position in space, we computed the orientation of the filtered velocity field . We healed the discontinuities of θ(r, t) so that . We performed a moving average of time window T0/2, with T0 = 2π/ω0 and ω0 the pulsation of the maximum of the power spectrum. At each position r and time t, we measure the sign ϵ(r, t) of the time variations of θ(r, t), which defines the spin field ϵ(r, t) = ±1. We illustrate this protocol in Extended Data Fig. 9 and show the resulting spin fields for the 2022 event, which reveal the same phenomenology as in Fig. 4. We estimated the error bars of the spin field distribution using the jackknife method over 10,000 samples38. For the correlation length of the spin field, we ran the band-pass filter on time windows of 3 min and computed the correlation length as detailed above.

　　In this section, we show that confinement and odd frictional forces conspire to yield spontaneous chiral oscillations at an incoherent speed in dense crowds.

　　Our starting point is the mean-field description of the crowd mechanics, which takes the form of Newton’s second law ∂tv = f/ρ, where ρ is the mass density and f is the sum of all the body forces originating from the interactions between the pedestrians and the ground or between the pedestrians and the confining walls. We decompose it into four terms: f(t) = −ργv(t) + ρp(t) − ρku + ρσζ(t), where the first drag term classically models the damping of the velocity fluctuations via the momentum exchange with the ground. In the spirit of a Landau expansion, we retain the lowest-order term in v and consider a linear drag force. The second term is what we refer to as an active friction with the ground. It is a vector that models the conversion of the body deformations into propulsive forces. The third term stems from the confinement of the crowd by walls. Again, in the spirit of a Landau expansion, we retain the lowest-order term in the displacement variable u (with v = ∂tu). The last term in the force definition is a Langevin noise source that classically accounts for the coupling to all the fast degrees of freedom ignored within our mean-field picture. ζ(t) is a Gaussian random noise of zero mean and covariance .

　　In Supplementary Information, we show that the constitutive relation that defines the propulsive force p cannot take the form p(t) = f(v(t), u(t)) (see also Extended Data Fig. 5a–c). In others words, we must take into account the proper dynamics of p. Because the crowd is globally isotropic (Extended Data Fig. 4) and has no reference position in space, the constitutive equation must be independent of u and invariant upon rotations. Furthermore, our experimental observations show that the parity symmetry of the crowd is not explicitly broken. We thus proceed to a systematic Landau expansion to find the most general constitutive equation. Restricting ourselves to first-order derivatives in time and keeping only the leading-order nonlinearities, corresponding to the six third-order terms, we find

　　where ηs are material parameters and ζp(t) is another Gaussian random noise of zero mean, covariance and uncorrelated from ζ(t). We discuss in Supplementary Information the effect of all terms and all nonlinearities in Newton’s second law and equation (6). We study them one at a time to single out their impact on the crowd dynamics. The conclusion of this thorough investigation is that our experimental findings are nicely captured by a minimal model where inertia plays a negligible role and where only two nonlinearities rule the dynamics of the p variable. The α2 term which we refer to as a ‘weathercock’ term in the main text is essential to yield chiral orientational oscillation. However a second nonlinearity is required to stabilize the oscillatory dynamics against fluctuations. More accurately, we show that the dynamics is unstable when all the ηi vanish. However, stable chiral states can emerge when η1 ≠ 0, η3 ≠ 0 or η5 ≠ 0. As the stability domain of the chiral states is much larger when stabilized by a nonlinear windsock effect (η3 and η5), we therefore focus on a minimal model where we disregard inertia and retain only the α2 and η3 nonlinearities.

　　The mechanics of dense crowds is then captured by a minimal set of two equations:

　　where we have set all ηi to zero, except η3 ≡ η/γp. The first equation corresponds to Newton’s second law and the second equation is the constitutive equation of the crowds. We provide a physical interpretation of this constitutive relation in the main text. The nonlinear η term implies that the amplitude of the windsock effect decreases as the amplitude of the propulsive force increases. We also note that this equation is strikingly similar to the mechanical description of active solid metamaterials derived in ref. 27 from a microscopic model. only the sign in front of the double cross-product is different.

　　Although equation (7) has a clear physical meaning, it is unpractical to investigate the fixed points, the limit cycles and the stability of the dynamics. To characterise the deterministic dynamics of our nonlinear system (σ = σp = 0), we write , with u the norm of u and its orientation. We then inject this form into equation (7), and express the dynamics in terms of the three variables u (the displacement magnitude), v = ∂tu (the ‘speed’) and (the square of the angular velocity of ). We then find

　　with βc = γ + k/γp. We stress here that the dynamical variable Ω2 does not prescribe the direction of rotation of u but only its rate of change.

　　We now look for the fixed points of the dynamics (u, v, ) by setting the partial derivatives to zero in equation (8). The first fixed point is trivial and corresponds to quiescent crowds where u = 0 and v = 0, whereas is ill defined as the displacement vanishes (u = 0). There also exists a second pair of non-trivial fixed points for β >βC为V = 0的特征

　　和

　　第二对固定点对应于U（和V）的两个圆极限循环。极限循环的半径为u（u手），并且它们在相反的角频率±ω上被扫描。为了确保这两个限制周期解释了我们人群的振荡动力学，我们必须首先证明它们对小扰动稳定。

　　现在，我们研究两个固定点的线性稳定性。对于微不足道的固定点，我们在无噪声的限制（u = 0，p = 0）的无噪声限制中线性化等式（7），然后找到

　　基质的痕迹是γP（β -βC）/γ。当β=βc时，它会改变符号。决定因素由KγP/γ给出，并且始终保持正。因此，我们得出的结论是，如果β，小固定点是稳定的< βc, as the matrix has two eigenvalues with a negative real part. Instead, the matrix has two eigenvalues of positive real part when β >βC和静态固定点因此变得不稳定。该标准与极限周期的出现相吻合。

　　现在，我们解决了两个极限周期的稳定性。我们设置u = u+ΔU，v =ΔV，并相对于ΔU，ΔV和ΔΩ2线性化上述方程。我们发现

　　Jacobian矩阵的特征值μ是立方方程的溶液-μ3+τμ2+νμ+δ= 0，τthe trace the trace和δ确定剂。特征多项式的三个系数取决于β并读取

　　如果三个特征值具有负实际部分，则极限周期是稳定的。这相当于强加τ< 0, Δ < 0 and ν + Δ/τ < 0. In Supplementary Information, we analyse extensively the conditions under which these three inequalities are verified. We show that there are ranges of parameters for which the limit cycles are stable for all values of β, leading to stable chiral states above βc. Instead, there are some other ranges of parameters for which the limit cycles can destabilize for some values of β. Our results are summarized in the stability diagram = (γγp/k, γη/α2) shown in Extended Data Fig. 10.

　　We have therefore shown that our mechanical model of confined dense crowds evolves according to a dynamics that cannot be captured by a steepest descent in an effective potential. This so-called non-reciprocal and nonlinear dynamics is here characterized by a mean-field phase transition (a bifurcation) from a trivial state to a chiral state where v chases p while p runs away from v along circular trajectories swept at constant rate. This stationary chiral dynamics can operate along two possible directions selected by initial conditions. The parity symmetry of the dynamics is spontaneously broken.

　　Our stability analysis shows that the propulsive force must relax faster than the displacement field to observe spontaneous oscillations: γp > k/γ. This finding a posteriori justifies the relevance of the asymptotic analysis γp  k/γ discussed in the main text, which we now detail (see also Supplementary Information). In this limit, we can ignore the time variations of the fast variable p, which is instantly slaved to u (in the absence of noise). We can then generically decompose p as p = (−K + k)u + Ku, and solve for K and K by using this form in the algebraic equation satisfied by p:

　　where we have used Newton’s second law to express v = ∂tu, as a function of u and p. We are then left with a single equation for the degree of freedom u that takes the form ∂tu = −Veff(u) ± (K/γ)u, where Veff(u) = Ku. The dynamics of the crowd is akin to that of a particle falling in an effective potential Veff and couples to a nonlinear odd spring K. Again the sign of K is not a priori specified and can take both values with equal probability. We sketch the dynamics of the particle in Fig. 5e. Deep in the chiral phase (β  βc), Veff has the shape of a ‘Mexican hat’. once the particle reaches the minimum of the potential, u takes a finite value. In turn, the particle dynamics orbits at constant speed along the ground-state circle under the action of of the odd spring force either in the clockwise or anticlockwise direction depending of the sign of K. This minimal model nicely captures the chiral dynamics of the dense crowds and provides an intuitive explanation for its nonlinear dynamics.

　　We solve equation (7) using a fourth order Runge–Kutta scheme39, with a time step δt (see below). We compute the velocity, orientation and speed spectra for nrun different realizations of the random noise forces and random initial conditions. All spectra are measured in a steady-state regime after nst preliminary integration steps specified below.

　　The definition of the spin ϵ(t) requires filtering the velocity field. To do so, we computed the mean velocity spectrum Sv over the nrun runs and performed a moving average of size nw = 10. We then defined our band-pass filter by considering only the angular frequencies ω for which the value of the smoothed spectra satisfies (quarter-difference between the maximum and the value of the spectrum at zero angular frequency). We finally used inverse fast Fourier transform to compute the filtered signals and . From these filtered signals, we then determined the angle that the filtered velocity field makes with the x axis: . The resulting discontinuous time series was defined on the unit circle. We then unwound this time series to define the orientation angle over . Finally, we defined the spin ϵ(t) as the sign of the difference between two consecutive values of the orientation angle.

　　We estimated the error bars associated with the power spectra and histograms of ϵ using the jackknife method over nrun measurements38.

　　Without loss of generality, we can set γ = 1 and α = 1, which amounts to defining our units of length and time. The model has then six control parameters, which are β, γp, k, η, σ and σp. To integrate the equations of motion, we set δt = 0.001, so that the damping coefficients in equations (1) and (2) verify k/γ  20δt and 1/γp  100δt for the typical values of γp and k that we consider. This value of δt also ensures that the angular frequency Ω always remains in our simulation window. The number of steps n is set so that the system explores exhaustively its phase space. Finally, the number of steps nst = 1.5 × 105 is set so that the system reaches the limit cycle when β >βC。

　　在没有噪声的情况下，当β>βc（即φ均匀地绘制在该范围内，而对于ω符号均等的概率）时，NRUN = 100个模拟在没有噪声的情况下初始化了一个极限循环。否则，模拟将从u = 0和p = 0初始化。

　　对于图3a – c，我们使用k = 0.027，γ= 1.00，γp= 18.00，β/βC= 1.10，η= 0.45，σ= 0.00和σp= 2.00。对于图3F，我们使用了γ= 1.00，γp= 18.00，η= 0.45，σ= 2.00和σp= 2.00。