Observation of Scaling of Higher-Order Flow in a Unitary Fermi Gas

Collective flow provides a powerful probe of emergent behavior in strongly interacting many-body systems. In relativistic heavy-ion collisions, the expansion of the hot and dense matter produced by colliding nuclei converts initial spatial anisotropies into correlations in the final momentum distribution [1]. The resulting azimuthal flow harmonics, conventionally denoted by $v_n$, have become central observables for understanding the transport properties and collective behavior of the quark–gluon plasma (QGP) [2, 3, 4, 5]. In particular, elliptic flow ($v_2$) and triangular flow ($v_3$) encode the response to spatial deformations and their harmonic dependence offers a sensitive probe of dissipative transport [6, 7, 8], as viscosity preferentially attenuates shorter-wavelength excitations [9]. Higher-order harmonics probe smaller spatial structures [10, 11, 12, 13, 14]; acoustic-scaling analyses organize their relative suppression approximately through $n^2/R$, where $R$ characterizes the initial transverse size [15, 16, 17]. However, the fluctuating initial geometry and the predominantly final-state access to momentum harmonics make it difficult to resolve how a selected harmonic develops and evolves during expansion [18, 19].

Ultracold Fermi gases (UFG) offer a unique experimental platform to address such emergent collective flow [20, 21, 22]. Their initial density geometry can be precisely engineered, and the temporal evolution of the expansion can be fully reconstructed. Furthermore, the well-characterized equation of state and control of temperature and density provide independent thermodynamic scales for testing hydrodynamic predictions [23, 24, 25]. At

a Feshbach resonance, the scattering length of two-component Fermi gases diverges, where interactions are as strong as allowed by quantum mechanics (so-called the unitarity limit) and the system becomes scale invariant. Their thermodynamic and transport properties are then governed by universal functions of density and temperature [26, 27]. Experiments have demonstrated the effects of strong interactions on anisotropic expansion [28, 29], scale invariant dynamics [30, 31, 32], shear viscosity [33, 34, 35], sound and density response [36, 37, 38, 39, 40]. More recently, interaction-induced elliptic flow has been observed in few-fermion systems [41], while comparisons between ultracold atomic gases and heavy-ion collisions have revealed common features of anisotropic expansion [42].

Despite these advances, the dynamical emergence of higher-order flow remains unexplored in ultracold Fermi gases. Most previous expansion studies have focused on global cloud deformation or second-order momentum anisotropy. Higher harmonics are particularly valuable because they probe spatial structures on progressively shorter length scales and therefore provide enhanced sensitivity to viscosity, nonlinear mode coupling and the breakdown of hydrodynamic behavior [12, 43]. Moreover, the influence of the initial density distribution can obscure the interpretation of collective flow in conventional harmonic confinement.

Here we program angular deformations of nearly homogeneous potentials to investigate the expansion of a unitary Fermi gas. We measure anisotropic flow harmonics up to the seventh order. We find that, after appropriate normalization, the flow dynamics exhibit a common early-time evolution across different harmonic orders and system sizes. This universal behavior subsequently breaks down, with the different harmonics developing distinct dynamics. We further identify a quadratic harmonic-order dependence of $n^{2}/R_{0}$ analogous to acoustic scaling in nuclear collisions in the measured response and establish a temperature-dependent dissipative contribution to the expansion. The temperature dependence and time evolution of this scaling provide new constraints on the interplay of pressure-driven conversion and transport. These observations reveal how universal flow generation evolves into harmonic-dependent transport and provide a connection between strongly interacting atomic gases and the acoustic attenuation of anisotropic flow in heavy-ion collisions.

Experiment.— We prepare a balanced mixture of the two lowest hyperfine states of ${}^{6}\mathrm{Li}$ close to the broad Feshbach resonance at a magnetic field $B \simeq832\ \mathrm{G}$, where the system enters the unitary regime. A digital micromirror device (DMD) generates a nearly uniform transverse density bounded by an engineered boundary shape Figure 1(a)

where $A$ sets the deformation amplitude and $\Phi_{n}$ its orientation. We vary $R_{0}$ and harmonic order $n$ from 3 to 7 while fixing $A/R_{0}=0.25$. For $n=2$, we prepare a uniform ellipse with a mean radius of $R_{0}$. Before release, the third direction along $z$-axes is harmonically trapped with frequency $\omega_{z}/(2\pi)=120(5)\ \mathrm{Hz}$. Simultaneous front and side imaging provides thermometry using the equation of state [23, 30]. The samples have typical Fermi energy $E_{\mathrm{F}}\sim k_{\mathrm{B}}\times350\ \mathrm{nK}$ and total atom number $N\sim1.2\times10^{5}$ for $R_{0}\sim30\ \mu\mathrm{m}$. We prepare dimensionless temperatures $T/T_{\mathrm{F}}\simeq0.15-0.35$, where $T_{\mathrm{F}}=E_{\mathrm{F}}/k_{\mathrm{B}}$. All optical confinement is suddenly removed, initiating interaction-driven expansion.

In contrast to harmonically trapped gases, the equilibrium density is nearly uniform and expansion begins through rarefaction waves launched from the boundaries [44], providing a well-defined measurement from spatial deformation to anisotropic collective flow.

The gas evolves for a controlled time of flight (TOF) $t$ at unitarity, and then a standard absorption imaging is performed. The effect from a weak residual magnetic saddle is corrected through a calibrated affine transformation.

For each image we compute the unweighted, signed angular moment of the column density $\Sigma(\mathbf{r},t)$,

with coordinates $\mathbf{r}$ centered on the cloud and a fixed final-state direction $\Psi_{n}^{\mathrm{f}}=\Phi_{n}+\pi/n$. We use the ratio $\mathcal{V}_{n}/\epsilon_{n}$ to quantify the efficiency with which the initial spatial anisotropy is converted into collective flow, where $\epsilon_{n}=\langle r^{n}e^{in\phi}\rangle/\langle r^{n}\rangle$ is evaluated from the programmed initial profile, adopting the conventional radial weight in nuclear-flow analyses [6, 16]. During early expansion, $\mathcal{V}_{n}$ measures spatial anisotropy; after sufficiently long TOF into the ballistic regime, it approaches the corresponding velocity-angle harmonic. At the experimentally accessible TOF, the observed evolution of $\mathcal{V}_{n}$ reflects the coupled density and velocity dynamics.

Geometric scaling and universal early-time response.— We present $\mathcal{V}_{n}(t)$ normalized to its initial magnitude in

Figure 2(a),(b), showing distinct dynamics for different $n$ and $R_0$, respectively. The anisotropy begins negative, crosses zero, and subsequently develops resolved maxima for the higher harmonics. At the same temperature $T/T_F=0.15(1)$, rescaling the TOF $t$ with $q=nt/R_0$ approximately collapses the early trajectories and further to the peak time for higher $n>3$ [Figure 2(c)], identifying a geometric timescale $R_0/n$.

At different temperatures, the trajectories remain distinguishable in $q$ [Figure 2(d)]. For nearly homogeneous unitary Fermi gas, the initial density-weighted sound speed is $c_s^2=5\int P_0\,d^3r/(3mN)$ giving an acoustic scale $\tau=c_sq$, where $P_0$ is the initial pressure and $m$ is the atomic mass. When plotted as functions of $\tau$, we observe that the early flow dynamics obtained for different $n$, sizes, and temperatures all collapse onto a common curve [Figure 2(e)], exhibiting the universal early-time response.

This universal behavior can be understood from the boundary driven hydrodynamics of the expansion. Immediately after the gas is released, the initially near uniform interior remains approximately unperturbed, while rarefaction waves propagate inward from the boundaries. At unitarity, a locally rarefaction wave follows a self-similar solution of the Euler equations. Its density and velocity profiles depend on distance from the released boundary only through the ratio of that distance to $c_st$. For zero initial flow and a weak deformation of a uniform transverse disk, $\mathcal{V}_n(0)\simeq-A/R_0$, $\dot{\mathcal{V}}_n(0)/|\mathcal{V}_n(0)|\simeq(n^2/R_0^2)\int P_0\,d^3r/(mN)$ and $\dot{\mathcal{V}}_n(0)=0$, this yields

Remarkably, all the rescaled data agree well with this universal prediction up to $\tau\approx1$. The subsequent separation of the trajectories, including their inversion times, probes the evolving pressure and flow during three-dimensional expansion.

Harmonic scaling at later times

At longer expansion times, the normalized flow harmonics depart from their common trajectory. The dynamics become harmonic dependent, revealing a crossover beyond the initial universal regime. Figure 3 compares measurements at two later common rescaled times $q_1<q_2$ as marked in Figure 2(c) at $T/T_F=0.15(1)$. Each slice includes two initial radii and harmonic orders up to $n=7$, and is evaluated at $t_j=q_jR_0/n$. Remarkably, the normalized angular moments at these times show an approximately exponential dependence on $n^2/R_0$:

Here $K_j$ is an effective attenuation length and $A_{0,j}$, the logarithmic intercept. Linear fits yield $K_1=1.22(3)\,\mu\mathrm{m}$ and $K_2=1.68(4)\,\mu\mathrm{m}$, with a goodness of fit $R^2$ all above $0.99. The inset of Figure 3(b) fixes $n=3$ and varies $R_0$ from $18$ to $54\,\mu\mathrm{m}$: the logarithmic response again scales approximately linearly with $1/R_0$ at fixed $T_F=249(6)\,\mathrm{nK}$ and $T/T_F=0.26(7)$. This further confirms the radius dependence over a wider range.

To test this scaling independently of the normalization convention [16, 18], we plot the two-time logarithmic ratio in Figure 3(c), $\ln[\mathcal{V}_n(q_1)/\mathcal{V}_n(q_2)]\simeq-\Delta A_0+\Delta K n^2/R_0$, where $\Delta A_0=A_{0,2}-A_{0,1}$ and $\Delta K=K_2-K_1$. The $\epsilon_n$ and any time-independent multiplicative normalization cancel. The observed linearity therefore confirms the $n^2/R_0$ dependence of the late-time attenuation independently of the normalization method.

The observed quadratic harmonic order dependence suggests the acoustic attenuation. For a nonrelativistic unitary Fermi gas, the shear-viscous contribution to the sound dispersion is

where $\eta$ denotes shear viscosity and $n_0$ is the initial atomic density. The amplitude of a weakly damped sound perturbation therefore decays as $A_k(t)\propto\exp[-Dk^2t]$, where $D=2\eta/3mn_0$. For a smooth angular harmonic of order $n$, the characteristic wave number is approximately $k(t)\simeq n/R(t)$. For a scale invariant expansion, $R(t)^2\simeq R_0^2+v_{\mathrm{exp}}^2t^2$ [30], where the expansion speed $v_{\mathrm{exp}}$ is on the order of $c_s$. Therefore, the accumulated attenuation exponent at fixed scaled time $q$ is given by $D_n\sim\int_0^{qR_0/n}dt'\,Dk^2(t')\sim(n^2D/R_0v_{\mathrm{exp}})\arctan(v_{\mathrm{exp}}q/n)$. In the late-time limit $v_{\mathrm{exp}}q/n\gg1$, this yields the shear-acoustic scaling

Wavelength stretching then makes further attenuation small, so the accumulated $n^2/R_0$ factor persists at later common $q$. The precise value of $D_n$ is sensitive to the sound diffusivity, the expansion history, and the duration of hydrodynamic evolution. Quantitative comparison with hydrodynamic calculations can disentangle these contributions and constrain transport coefficients.

The $n^2/R_0$ dependence is particularly noteworthy in light of the acoustic scaling of the final momentum anisotropy in relativistic heavy-ion collisions [15, 16]. In that setting, the damping of flow harmonics is governed by the shear viscosity relative to the relativistic entropy density. For a nonrelativistic Fermi gas, the corresponding inertial quantity is the mass density $mn_0$. However, the expanding potential differs from a small amplitude sound-wave problem. The flow harmonics are generated dynamically from an initially isotropic velocity distribution, and the sharp boundaries excite the dynamical expansion with spatial wavelengths. Thus, the observed quadratic scaling has to be understood as an effective collective response rather than a single-mode damping law.

Temperature dependence of harmonic attenuation.—We prepare various temperatures at $R_0\simeq31\ \mu\mathrm{m}$ for $n=2-7$ at fixed $T_F=352(6)\ \mathrm{nK}$. Figure 4(a) shows the harmonic scalings at $q_1$ with the fitted slopes and intercepts shown in (b),(c) as a function of the dimensionless temperature.

The fits proportional to $T/T_F$ and $(T/T_F)^{1/2}$ in Figure 4(b) illustrate that our data are consistent with a range of temperature scaling of $K$. The linear dependence, on the one hand, has a simple high-temperature argument from Eq. (6): when attenuation dominates the harmonic slope, $K$ is proportional to $D/v_{\mathrm{exp}}$. At temperatures above the superfluid transition, the gas is a strongly interacting normal quantum fluid. Its shear viscosity can be expressed in a universal form. The high temperature unitary-gas transport gives $\eta\propto T^{3/2}$ at fixed density and $c_s\propto T^{1/2}$ [26, 37, 45], hence

where $k_F=(3\pi^2n_0)^{1/3}$ is nearly constant for these preparations. On the other hand, the explored temperature range $T/T_F \simeq0.15$–$0.35$ is nevertheless highly degenerate, so Pauli exclusion is expected to play a role that affects the attenuation through both viscosity, expansion speed, and thermal transport. An effective temperature scaling closer to $(T/T_F)^{1/2}$ in this temperature range is also consistent with previous measurements of the sound diffusivity [37, 38]. Our measurements thus provide new constraints on the interplay of pressure-driven conversion and transport.

Conclusions.— Our measurements revealed the emergence and subsequent evolution of higher-order collective flow in a strongly interacting, scale-invariant, nearly homogeneous quantum fluid. Programmable angular deformations revealed an early universal acoustic response and later $n^2/R_0$ scaling for harmonic orders up to $n = 7$. Ideal hydrodynamics quantitatively explained the initial quadratic evolution, while viscous acoustic attenuation captured the late-time harmonic scaling and its temperature dependence.

Looking forward, harmonic-resolved expansion measurements could provide access to temperature dependent transport coefficients to test string theory limit [21, 46] and nonlinear flow couplings, including the conversion of lower-order modes into higher harmonics. The

approach developed here could also be straightforwardly generalized across the entire BEC–BCS crossover regime, enabling the investigation of how changes in interaction strength and scale invariance affect the emergence of collective flow. This work thereby extends quantum gas experiments from elliptic flow to higher-order anisotropic flows, opening up new opportunities to connect programmable ultracold quantum gases with collective dynamics in strongly interacting nuclear matter.

Acknowledgments

Acknowledgments.— The project was supported by National Key R&D Program of China (No. 2022YFA1404202), National Science Foundation of China (No. 12234008, No. 12688101, No. 12547102), the Shanghai Municipal Science and Technology Major Project (No. 2019SHZDZX01, No. 25LZ2600400), Innovation Program for Quantum Science and Technology (No. 2021ZD0302100), and Shanghai QiYuan Innovation Foundation.

References