Consider a classical system of stochastic degrees x t at time t (weakly) coupled to a heat bath at temperature T . The system’s potential energy U depends on x t and on a time-dependent deterministic protocol X t , i.e., U ( x t , X t ) . The system is prepared in equilibrium at time t → − ∞ , with a protocol value X − ∞ = X t , for simplicity [62]. The time dependence of the protocol drives the system away from equilibrium.

The derivative of U in terms of X t , F t : = ∂ X t U ( x t , X t ) is the observable conjugate to X t . For example, if X t is a position as in our experiments, F t is (minus) the corresponding force. F t can be micro- or macroscopic; for example, let X couple linearly to an observable A ( x ) ; U ( x t , X t ) = U ( x t , 0 ) − X t A ( x t ) , i.e., F t = − A ( x t ) . Thus, if A ( x ) is a macroscopic field, such as macroscopic magnetization, F t is macroscopic. If A ( x ) is the position of a molecular particle, F t is microscopic. The following remains valid if X t enters U nonlinearly.

The statistical properties of F t in this nonequilibrium situation are encoded in its cumulants and in its correlations with another state observable, B t = B ( x t , X t ) , which we aim to study here. The well-known FDT connects the covariance in the unperturbed system and the first moment of B under weak driving [1, 2], β 2 ∫ − ∞ t d s   X ̇ s ⟨ B t ; F s ⟩ eq = β ⟨ B t ⟩ − ⟨ B t ⟩ eq + O ( X ̇ 2 ) , (1) with β = 1 / k B T and Boltzmann constant k B . ⟨ … ; …   ⟩ denotes the second cumulant, and similarly for higher orders below. Expectation values ⟨ …   ⟩ eq are evaluated using equilibrium trajectories with a protocol value fixed at X t . Expectation values ⟨ …   ⟩ are measured in the driven system, i.e., under the given time-dependent protocol [62].

In Ref. [62], we derive identities connecting the nonequilibrium cumulants of F t and B t to different orders. These give rise to the following series involving the mentioned cumulants [62], β 2 ∫ − ∞ t d s   X ̇ s   ⟨ B t ; F s ⟩ = β ⟨ B t ⟩ − ⟨ B t ⟩ eq + β 3 2 ∫ − ∞ t d s ∫ − ∞ t d s ′   X ̇ s X ̇ s ′ ⟨ B t ; F s ; F s ′ ⟩ − β 4 6 ∫ − ∞ t d s ∫ − ∞ t d s ′ ∫ − ∞ t d s ″   X ̇ s X ̇ s ′ X ̇ s ″ × ⟨ B t ; F s ; F s ′ ; F s ″ ⟩ + O ( X ̇ 4 ) , (2)

As presented in Ref. [62], this series expansion can also be obtained from a known fluctuation theorem [33, 41], albeit with an open question regarding cumulants versus moments.

Equation 2 is, as indicated, correct up to the fourth order in driving X ̇ t , under the assumption of local detailed balance [63]. Expanding ¹ Equation 2 to the first order yields the FDT in Equation 1 so that it is included in Equation 2. To higher orders, first and second cumulants do not fulfill the FDT, and Equation 2 quantifies their difference in terms of third and fourth cumulants of F and B . Notably, the second to fourth lines of Equation 2 vanish for purely Gaussian distributed F and B so that first and second cumulants obey the FDT to the given nonequilibrium order. It is important to note that, in Equation 2, the protocol X ̇ appears as pre-factors and in the nonequilibrium cumulants themselves; i.e., the latter are evaluated under application of driving. As Equation 2 only requires measurement of F and B , we use the notion of model free.

We exploit Equation 2 with experiments of Brownian particles interacting with micellar fluid. In particular, we use silica particles of diameter ∼ 1 µm suspended in a 5-mM equimolar solution of cetylpyridinium chloride monohydrate (CPyCl) and sodium salicylate (NaSal). At concentrations above the critical micellar concentration ( ≳ 4 mM), this fluid is known to form giant worm-like micelles, leading to a viscoelastic nonlinear behavior at ambient temperatures [64]; see SM. At 5 mM, we determine the relaxation time of the fluid from microrheological recoil experiments, where a particle is first driven with a constant external force which is then suddenly removed, to be ∼ (3 ± 0.2)s [65]. A small amount of silica particles is added to the micellar solution, which is contained in a rectangular capillary with 100 µm height and kept at a temperature of 25 °C. This sample is placed on a custom-built optical tweezer setup that uses a Gaussian laser beam of wavelength 532 nm and a 100 × oil immersion objective (NA = 1.45). The laser beam yields a potential U ( x t − X t ) , as shown in Figure 1a, centered at X t , trapping one of the silica particles with coordinate x t . F t = ∂ X t U ( x t − X t ) is thus the force acting on the particle by the trapping potential (or vice versa). As the micellar degrees do not couple to X , they do not enter F t explicitly, and knowing how they enter U is not required to apply Equation 2. Thus, the use of Equation 2 does not require the detection of the positions of micellar particles, and applying it to such a complex fluid demonstrates its strength. We consider B ≡ F , and made the potential asymmetric, to obtain a finite second-order response of β ⟨ F t ⟩ . This allows to test Equation 2 to the second order in our experiments, and it is achieved by a controlled lateral displacement of the vertically incident laser beam from the center of the objective lens (see SM). Notably, Equation 2 could also be tested in a purely viscous fluid, which also shows a second-order response due to the nonlinear potential. However, to demonstrate that Equation 2 is valid beyond simple systems, we have chosen the more challenging case of a micellar fluid.

To apply the driving protocol, the sample cell is moved, whereas the optical trap remains stationary in our experiments. This is achieved using a piezo-driven stage, on which the sample is mounted and translated in an oscillating manner relative to the trap. In the fluid’s rest frame, this yields a periodic motion of the potential minimum X t , i.e., the protocol, X t = X ̂ sin ( ω t ) , (3) with the amplitude X ̂ and the frequency ω . Particle trajectories are recorded with a frame rate of ∼ 150 Hz using a video camera, and particle positions are determined using a custom MATLAB algorithm. To yield sufficient statistics, each protocol ( X ̂ , ω ) was measured over 1400s. We allowed the system to reach a steady state by recording trajectories only after at least five oscillation periods had passed. Thereafter, no further equilibration was visible in the data. Prior to each nonequilibrium protocol, we recorded particle trajectories for another 1000s with X t at rest. These equilibrium data were used to check that the form of U does not vary between measurements and also to obtain the force cumulants under equilibrium conditions.

With the protocol of Equation 3, Equation 2 takes, expanded to the second order, the form β 2 X ̂ ω ⟨ F ̃ t ; F t ⟩ = β ⟨ F t ⟩ + β 3 X ̂ 2 ω 2 2 ⟨ F ̃ t ; F ̃ t ; F t ⟩ + O ( ( X ̂ ω ) 3 ) , (4) where the tilde denotes the cosine transform, i.e., F ̃ t ≡ ∫ − ∞ t d s   cos ( ω s ) F s . We restrict the analysis to the lowest nontrivial, i.e., second, order and expand Equation 2 accordingly, also using ⟨ F t ⟩ eq = 0 . Notably, to extract the second-order contribution from the last term in Equation 4, the third cumulant in Equation 4 is replaced by its equilibrium version. This is because of the pre-factor X ̂ 2 ω 2 and will be done in the following analysis.

The cumulants in Equation 4 depend on time t in a periodic manner, as shown in Figures 1b–d, for ω = 8.4 rad/s and driving ² amplitudes ranging from X ̂ = 0.03 µm (light) to X ̂ = 0.14 µm (dark). Figure 1b shows the mean force, which, as expected for a driven oscillator, is a periodic function with period T p = 2 π ω . For the smallest amplitude shown, the mean force is nearly harmonic with frequency ω , as expected from linear response. With growing amplitude, higher harmonics occur, as expected from nonlinear response. This asymmetric system shows the second-order response with expected frequencies of 2 ω and 0 ω .

Figure 1c shows the force covariance for the same parameters and color code. For small amplitude X ̂ , the curves in Figures 1b,c are equal within the experimental accuracy, as analyzed in detail below. For larger driving amplitude, the force covariance develops higher harmonics with signatures of second order. Very little is known about the properties of such nonequilibrium fluctuations, and quantifying these is difficult. It is notable that the curves in Figure 1c, for larger amplitudes, deviate from Figure 1b, a deviation that we claim to be quantified by Equation 4.

Figure 1d shows the third cumulant of force for the same parameters and color code. We have here restricted to the equilibrium cumulant as it appears in Equation 4, multiplied by ( X ̂ ω ) 2 . The curves in Figure 1d thus differ only because of the factor ( X ̂ ω ) 2 . They thus scale quadratically in driving velocity and only show frequencies of 2 ω and 0 ω .

Equation 4 states that in the shown range of amplitudes, the curves in Figure 1c are given by the sum of the curves in Figures 1b,d. For X ̂ = 0.14 µm, the respective summed curve is shown as a dashed line together with the mean force and the force covariance in Figure 1e. The agreement is convincing and a confirmation of Equation 4.

To test this prediction systematically, we dissect the curves in Figures 1b–d into the contributions from harmonics with frequencies 0 ω , ω , and 2 ω , respectively, i.e., we expand the cumulant of order n into harmonics with the frequency m ω , β n ( X ̂ ω ) n − 1 ⟨ ( F ̃ t ; ) n − 1 F t ⟩ = ∑ m = 0 ∞ A m ( n ) ⁡ sin ( m ω t + ϕ m ( n ) ) , (5) where the coefficients A m ( n ) depend on X ̂ ω . We set ϕ 0 ( n ) ≡ π / 2 for consistency. Equation 4, projected on the harmonic of order m , yields relations between coefficients and phases for each m , which we can test.