How dissipation constrains fluctuations in nonequilibrium liquids: Diffusion, structure and biased interactions

The dynamics and structure of nonequilibrium liquids, driven by non-conservative forces which can be either external or internal, generically hold the signature of the net dissipation of energy in the thermostat. Yet, disentangling precisely how dissipation changes collective effects remains challenging in many-body systems due to the complex interplay between driving and particle interactions. First, we combine explicit coarse-graining and stochastic calculus to obtain simple relations between diffusion, density correlations and dissipation in nonequilibrium liquids. Based on these results, we consider large-deviation biased ensembles where trajectories mimic the effect of an external drive. The choice of the biasing function is informed by the connection between dissipation and structure derived in the first part. Using analytical and computational techniques, we show that biasing trajectories effectively renormalizes interactions in a controlled manner, thus providing intuition on how driving forces can lead to spatial organization and collective dynamics. Altogether, our results show how tuning dissipation provides a route to alter the structure and dynamics of liquids and soft materials.


I. INTRODUCTION
Nonequilibrium forces can drive novel and specific pathways to modulate phase transitions and selfassembly in materials. The close connection between the net dissipation of energy, powered by these forces, internal transport and spatial organization is especially apparent in living systems [1][2][3][4]. As an example, the flagella motors of E. Coli exhibit a unique phenomenology combining ultra-sensitive response, adaptation, and motor restructuring as a function of the applied load [5][6][7]. Moreover, in vivo studies of the cellular cytoskeleton, as well as in vitro experiments on reconstituted systems, have also shown that motor-induced forces control a large variety of functionality in the cell [8][9][10][11][12].
To elucidate the role of nonequilibrium forces in materials, it is crucial to examine how dissipation affects the emerging dynamics and structure. While equilibrium features are well established, progress in controlling systems with sustained dissipation has been hampered by a lack of general principles [13][14][15][16][17][18][19]. In this context, minimal models of active and driven systems provide analytically and numerically tractable test beds to investigate the interplay between dissipation and material properties far from equilibrium [20][21][22][23][24]. They have illustrated, for instance, how nonequilibrium driving can induce phase transitions and excite novel collective responses in soft media [15,21,[25][26][27]. Recent theoretical work has proposed extending equilibrium concepts to active media, such as the definition of pressure [14,28], to rationalize their phenomenology [29,30]. Others have striven to obtain stationary properties of active matter through perturbation close to equilibrium [16,31,32], inspired by other approaches on driven systems [33][34][35][36].
To investigate how dissipation controls emerging behavior, yet another approach has focused on introducing a bias in dynamical ensembles. Using large deviation techniques, trajectories are conditioned to promote atypical realizations of the dynamics [37,38]. Such techniques have been used, for instance, to investigate the role of dynamical heterogeneities in glassy systems [39][40][41][42][43][44][45] and soliton solutions in high-dimensional chaotic chains [46,47]. More recently, it has been shown that changing dissipation, through a dynamical bias, strongly affects the internal transport and the density fluctuations of nonequilibrium liquids [48,49], thus confirming that controlling dissipation is indeed a fruitful route to tailoring material properties. In spite of these advances, anticipating the emergent dynamics and structure of biased nonequilibrium systems is still challenging in the presence of many-body interactions [38,50], so that precise control has remained elusive so far in this context. Consequently, any generic principle rationalizing spatial organization in terms of dissipation is still lacking.
In this paper, we explore how dissipation affects the dynamics and structure of many-body diffusive systems. First, we consider in Sec. II two types of assemblies of Brownian particles: one in which only a subset is driven by an external force, and one in which a subset of the particles experience an internal active force. We first focus on instances where the fraction of driven particles is less than the fraction of undriven particles, so that driven and undriven particles are respectively referred to as tracer and bath particles. Using the diffusion coefficient of a tagged tracer particle and the density correlations between tracer and bath particles, we connect dissipation to liquid properties. In contrast with [23], our prediction for diffusion follows from a systematic coarse-graining with explicit dependence in terms of microscopic details [51][52][53].
Next, importantly, we put forward a generic relation between density correlations and dissipation valid for an arbitrary driving force: this is our first main result. We demonstrate that this result holds both for fluids in which a fraction of the particles are driven by a fixed external drive and for fluids in which either a fraction of the liquid or the entire liquid is driven by an internal noise, analogous to the driving used in model active matter systems. This result opens the door to estimating dissipation directly from the liquid structure, in contrast with previous approaches based either on perturbing the system [54][55][56][57][58] or on analyzing trajectories and currents in phase space [3,[59][60][61][62]. We illustrate this with numerical simulations for which dissipation is quantified by the deviation from equilibrium tracer-bath correlations. Using these results as a basis, we also show how various aspects of the pair correlation function of a nonequilibrium liquid are effectively constrained by the energy dissipation. Altogether, this set of results clarify how nonequilibrium forces affect the transport and structure of the liquid, thus showing how liquid properties can be modified at the cost of energy dissipation.
Motivated by these results, and to provide concrete intuition for how particular configurations can be stabilized by nonequilibrium forces, we next investigate in Sec. III the emerging structure of Brownian particles subject to a dynamical bias. The explicit form of the bias is inspired by the results of Sec. II connecting dissipation to many-body interactions. Using analytical calculations and numerical simulations based on the cloning algorithm [46,[63][64][65][66][67][68], we show that biased sampling trajectories can be used to renormalize any specific interparticle interaction in a multi-component liquid. The rare noise fluctuations sampled with dynamical bias effectively drive the system away from typical behavior [38][39][40][41][42][43][44][45]50]. Such noise realizations can then serve as proxies of how to control the dynamics by applying an external force with complex protocols. We also illustrate the generality of our ideas by considering an assembly of aligning self-propelled particles [69]. Specifically, we show how biased energy flows can renormalize interactions between particles and stabilize a flocking transition. Overall, our results lay the groundwork for precise control of the emerging structure and collective dynamics in many-body diffusive nonequilibrium systems.

II. DISSIPATION AND LIQUID PROPERTIES
In this Section, we provide a series of relations between energy dissipation and liquid properties in nonequilibrium liquids. Specifically, we consider interacting Brownian particles where a specific set of particles Ω is driven by a non-conservative force F d,i : where δ i∈Ω = 1 if i ∈ Ω and δ i∈Ω = 0 otherwise. The driven particles belonging to the set Ω are referred to as tracers, and others as bath particles. The fluctuating term ξ i is a zero-mean Gaussian white noise with correlations ξ iα (t)ξ jβ (0) = 2γT δ ij δ αβ δ(t), where γ and T respectively denote the damping coefficient and the bath temperature, with the Boltzmann constant set to unity (k B = 1).

A. Deterministic vs. active drive
In what follows, we consider two types of drive: (i) an external force following the same deterministic protocol for all driven particles, and (ii) an internal force given by a noise term independent for each driven particle. Building on recent work [21,23], we take for drive (i) a time-periodic protocol given in two dimensions by where f and ω are respectively the amplitude and the frequency of the drive, so that the drive persistence reads τ = 2π/ω. The relative strength of the drive is given by the Péclet number Pe = σf /T , where σ is the typical particle size [21,23]. In the absence of interactions (v = 0), the average position of driven tracers follows a periodic orbit, describing a circle in two dimensions. In contrast, drive (ii) corresponds to a random self-propulsion as is often considered in active liquids [70][71][72]. Specifically, we use a set of zero-mean Gaussian noises with correlations where d is the spatial dimension. The parameters f and τ respectively control the amplitude and the persistence of fluctuations. Interestingly, the active force with correlations (3) can be obtained from a generalized version of the deterministic force (2) where each particle i is now subjected to an independent drive. The period of the orbit is determined by a series of n oscillators with identical frequencies for all particles, yet independent amplitudes: The essential ingredient of the mapping into active force is to implement disorder in the drive. This is done by taking the oscillator amplitudes as uncorrelated zero-mean Gaussian variables with unit variance: where · d denotes an average over the disorder. It follows that F d,i is also a Gaussian process with zero mean and correlations given by: cos(ω a t). (6) In the limit of a large number of oscillators (n 1), we express these correlations in terms of the density of driving frequencies φ as This establishes that, in the limit of many oscillators, the deterministic drive (4) with disordered amplitude is equivalent to a noise term with spectrum φ. In particular, by choosing φ(ω ) = 2τ / 1 + (ω τ ) 2 , the drive correlations (7) reproduce exactly the ones of the random force in (3).
To illustrate the relevance of this mapping, we simulate numerically the many-body dynamics (1) where every particle is subjected to a disordered drive of the form (4). We use the potential v(r) where Θ denotes the Heaviside step function, which sets purely repulsive interactions. To implement numerically the disorder in driving, it is sufficient to sample the amplitudes {A ai , B ai } and frequencies {ω i } at initial time.
In the regime of high persistence τ and large average density ρ 0 , we observe the spontaneous formation of clusters up to a complete phase separation at large time, see Fig. 1. This is analogous to the motility-induced phase separation commonly reported in standard models of active particles [13,26]. Interestingly, it appears in our case even in the absence of fluctuations (T = 0), namely for a purely deterministic set of equations.
In short, we thus demonstrate that the disordered drive alone reproduces the emerging physics of active systems. This important result bridges the gap between two main classes of nonequilibrium liquids, where the driving force stems from either a deterministic protocol or a random noise. In what follows, we obtain analytic and numerical results for both drives to illustrate the broad applicability of our framework, ranging from systems driven by deterministic fields to active matter systems.

B. Dissipation controls tracer diffusion
To connect tracer diffusion with dissipation, we first describe the dynamics of undriven particles in terms of a coarse-grained variable. Using standard techniques, the dynamics of the density field ρ(r, t) = i ∈Ω δ[r − r i (t)] can be written as a non-linear Langevin equation [51]. In the regime of weak interactions, the density fluctuations δρ(r, t) = ρ(r, t) − ρ 0 around the average density ρ 0 are Gaussian and captured by the following Hamiltonian [53, Snapshot of particles subjected to a disordered drive. A phase separation emerges which is analogous to the motility-induced phase separation of active particles [13,26]. Simulation details in Appendix A and movie in [73].

74, 75]:
where K(r) = δ(r)/ρ 0 +v(r)/T . Note that density fluctuations remain generally Gaussian even for a homogeneous active liquid [70]. The conserved density dynamics reads ∂δρ(r, t) ∂t = D G ∇ 2ˆK (r − r )δρ(r , t)dr where D G = ρ 0 T /γ and γ G = γ/ρ 0 are respectively the field diffusion coefficient and the field damping coefficient. The term Λ is a zero-mean Gaussian white noise Owing to the linearity of the density dynamics (9), it can be readily written in Fourier space δρ(q, t) = ρ(r, t)e −iq·r dr as so that the field dynamics can be directly solved as (11) Considering the limit of dilute driven tracers, where interactions among them are negligible, their dynamics reads with´q =´dq/(2π) d . As a result, (11) and (12) provide closed time-evolution equations for tracers only. It should only be valid for weak interactions a priori, yet previous works have shown that it remains qualitatively relevant even beyond this regime in practice [76][77][78]. Indeed, Gaussian field theories for density fluctuations provide a very good description of simple liquids [74].
To characterize the transport properties of the liquid in the presence of driving forces, our first goal is to obtain an explicit expression, in terms of microscopic details, for the tracer diffusion coefficient: We aim to explore connections between D and dissipation, which is defined from stochastic thermodynamics as the power of the forces exerted by all tracers on solvent: where · denotes a Stratonovich product [79,80]. Dissipation is directly related to entropy production, as a measure of irreversibility, both when the drive is deterministic [79,80] and when it is a correlated noise [81][82][83][84]. Substituting the dynamics (1), the dissipation coincides with the power of driving forces: J = i∈Ω ṙ i · F d,i . Besides, replacingṙ i by its expression in (1), and using the fact that ξ i and F d,i are uncorrelated, we deduce that the dissipation can be further separated into free-motion and interaction contributions as Given thatẇ is the only non-trivial contribution to dissipation, connecting diffusion and dissipation simply amounts to expressing D in terms ofẇ. Deriving transport coefficients in nonequilibrium many-body systems, whose collective effects result from the complex interplay between interaction and driving forces, is a notoriously difficult task [85][86][87][88][89]. We set up a proper perturbation scheme by scaling the pair potential v with a small dimensionless parameter h 1 which controls the coupling between tracer and bath. In Appendix B, we obtain some explicit expressions for D andẇ to quadratic order in h and in the scaled driving amplitude Pe.
First, we discuss the case of the deterministic drive (2), and we focus on the limits of small and large driving frequency, respectively ωτ r 1 and ωτ r 1, where the relaxation time scale τ r = (D G /σ 2 )K(|q| = 1/σ) is set by density diffusion over the tracer size σ. First, at high frequencies ωτ r 1, the rate of work per particleẇ/N and the deviation from equilibrium diffusion D − D eq , where D eq is the diffusion coefficient for Pe = 0, are given byẇ .
(15) In the opposite limit of low frequencies ωτ r 1, we geṫ Bothẇ/N and D − D eq are now independent of the driving frequency ω. As a result, our perturbation theory shows that the scalings ofẇ and D − D eq are identical, both in terms of the drive amplitude Pe and of its frequency ω, in asymptotic frequency regimes. Note that the scaled rate of work γẇ/(N f 2 ) coincides with the reduced equilibrium diffusion γD eq /T − 1 to this order [52,53], as expected from linear response. The case of the active drive with correlations (3) follows by using the mapping between disordered drive and active forcing in Sec. II A. In practice, we first derive the diffusion coefficient D and the rate of workẇ for the driving force (4) at fixed disorder, as a straightforward generalization of the deterministic driving case, and we then average over the disorder. At small persistence τ τ r , we geṫ In contrast, the large persistence limit τ τ r yields the same results as for the low frequencies regime of deterministic drive, namely the expressions (16). Indeed, the force F d,i has a constant direction in such a limit, and the difference between deterministic and active drives, which respectively correspond to independent or similar directions for each tracer, is irrelevant in the limit of dilute tracers.
When the size a of the bath particles is significantly smaller than the tracer size σ a, which amounts to setting different pair potential v for bath-bath and for bath-tracer interactions, one can safely neglect the variation of K(q) in (15)(16)(17), so that K(q) K(|q| = 1/a). Then, in both regimes ωτ r 1 (τ τ r ) and ωτ r 1 (τ τ r ), the renormalization of the diffusion coefficient D − D eq can be simply written in terms of the rate of work per particleẇ/N for Pe 1 as Thus, the excess rate at which tracers move over their own size compared to equilibrium, set by the lhs of (18), is controlled by the rate at which work is applied on tracers by nonequilibrium forces, set by the rhs of (18). The proportionality factor depends on the details of interactions and of density fluctuations. Interestingly, this result is valid both for deterministic and active drives. It corroborates numerical observations obtained previously in a system where composition-dependent diffusion constants can lead to phase transitions [23].

C. Dissipation sets density correlations
We now explore how dissipation relates with static density correlations of the liquid. To this end, we treat undriven bath particles without any approximation in what follows, instead of relying on the Gaussian density field theory for δρ as in Sec. II B, and we consider an arbitrary set of driving forces F d,i . In equilibrium, the liquid structure can be derived from a hierarchy of equations for density correlations, whose explicit form reflects the steady-state condition on the many-body distribution function [90]. In our settings, steady-state conditions should now provide modified equations for density correlations, which can potentially make apparent the connection with dissipation.
This motivates us to consider the average rate at which the potential U = i∈Ω,j v(r i − r j ) changes, which can be written using Itô calculus as Substituting the dynamics (1) and using ξ i · ∇ i v = 0 within Itô convention, we get In the first line of (20), we recognize the rate of workẇ as defined in (14), and the term γẇ act = {i,j}∈Ω F d,i · ∇ i v(r i − r j ) which quantifies the contribution of interactions among driven particles to dissipation. The latter vanishes exactly when the drive is identical for all particles, since {i,j}∈Ω ∇ i v(r i − r j ) = 0 by symmetry, and it can be neglected for active drive when the fraction of driven particles is small. Then, using the steady-state condition U = 0, we deducė and denotes a sum without the overlap of indices: i = j, k = i and k = j. The power balance (21), valid for an arbitrary driving, either deterministic or active, is our first main result. Importantly, it holds for generic interactions and for any number of driven particles, namely not only in the limit of dilute tracers, in contrast with the results in Sec. II B.
In practice, it reflects how density correlations adapt to the presence of nonequilibrium forces. For a vanishing rate of work (ẇ = 0 =ẇ act ), one recovers the first order of the equilibrium Yvon-Born-Green (YBG) hierarchy, in its integral form, for two-component fluids [90]. At finite rate of work (ẇ = 0), the relation between the two-body correlation g and the three-body terms {g 3a , g 3b } is now implicitly constrained by dissipation. A direct implication is that the rate of work can now be inferred simply by measuring static density correlations, provided that the pair-wise interaction potential is known, for a given driven liquid. Importantly, such an approach does not require any invasive methods based on comparing fluctuations and response [54][55][56][57][58], and it does not rely on a detailed analysis of particle trajectories [60,61] or currents in phase space [59,62], whose experimental implementation can require elaborate techniques [3,4].
However, the power balance (21) is not straightforward to test, either numerically or experimentally, due to the three-body correlations. In equilibrium, where tracer and bath particles are indistinguishable, we get g 3a = g 3b . Assuming that this remains approximately valid in the driven case for a small fraction of tracers, the rate of work can simply be written in terms of the force exerted on a tracer To probe the validity of this result, we simulate the dynamics (1) where 10% of particles are subject to the Parametric plot of the rate of workẇ/T and of the statistics of bath-tracer forces i∈Ω F 2 i + T ∇i · Fi /(γT ) when 10% of particles are driven by either (a) a deterministic force, or (b) an active force. The solid line with slope 2 refers to the approximate relation (23). The satisfying agreement with numerical data indicates that the rate of work can be estimated by only measuring bath-tracer forces. The simulations were performed with N = 4500 particles using the procedure described in Appendix A. Parameters: Pe = 12 ( ), 18 ( ), 24 ( ), 30 ( ), 36 ( ); (a) τ T /(γσ 2 ) = 2 × 10 −1 (black), 3 × 10 −1 (brown), 4 × 10 −1 (red), 5 × 10 −1 (orange); driving force, considering either the deterministic periodic drive (2) or the active noise drive (3). Interactions are set by the WCA potential v(r) = 4v 0 (σ/|r|) 12 − (σ/|r|) 6 Θ(2 1/6 σ − |r|) [91]. Our measurements in Fig. 2 show that (23) is indeed a good approximation at small Pe and small τ , namely when the drive only weakly perturbs the liquid. The discrepancy is higher for the active case compared with the deterministic one, sinceẇ act = 0 in the latter without any approximation. In contrast with previous approaches [3,54,92], which rely on prospecting the whole system, our results demonstrate that the rate of work can actually be evaluated with only a small error by considering solely forces acting on tracer: the contribution of forces on other particles is negligible for a small fraction of driven tracers.
To evaluate further the change in liquid structure induced by dissipation, we measure the deviation from equilibrium pair correlations g − g eq due to the driving forces (left column in Fig. 3). In particular, inspired by the two-body contribution in the power balance (21), we focus on the observable At a given τ , scaling I by Pe 2 reveals that all curves almost collapse into a master curve for our numerical range Pe ∈ [12,36], as reported in the middle column of Fig. 3. In practice, particles overlap more for a stronger drive, so that g departs from zero at smaller interparticle distance. To correct for this, we introduce a shift of the curves I(|r|) as |r| → |r| + a, where a(Pe) is a fitting parameter. Given that the rate of work also scales like Pe 2 , it suggests the existence of an underlying relation between´I(r)dr andẇ. In practice, a linear fitting provides a satisfactory agreement between them, as shown in the right column of Fig. 3: where α is a fitting parameter independent of the Péclet number. This empirical relation demonstrates that, in the limit of dilute tracers, the rate of work can actually be directly estimated by comparing driven and equilibrium pair correlations for both deterministic and active drives. Comparing (21) and (24), we deduce the following integral relation between density correlationŝ (25) Interestingly, it is reminiscent again of the connection between density correlations provided by the YBG hierarchy at equilibrium [90]. Similarly, the relation (25) amounts to a constraint on density correlations, now valid for nonequilibrium liquids, which could guide the search for explicit predictions on the emerging structure. Importantly, it does not rely on any equilibrium mapping, in contrast with previous works [93][94][95], since it remains valid for non-negligible dissipation.
The power balance (21) can actually be extended to the case where all particles in the liquid are driven aṡ where g and g 3 now refer respectively to the two-body and three-body density correlations among all particles. This leads to an exact relation between the rate of work and the forces applied to particles F i aṡ which differs from the relation (23) for driven tracers by an overall factor of 2. A result analogous to (27) was found previously for deterministic drive [36,96]. The main difference is that (27) only features interaction forces F i in the rhs, thus allowing one to evaluate the rate of work without any prior knowledge on the driving force. Besides, it is valid for both deterministic and active drives. Moreover, conducting the same analysis of density correlations as for driven tracers, I exhibits again a scaling with Pe 2 , as reported in Fig. 4. We show thaṫ w and´I(r)dr are also linearly related. Introducing the linear coefficient as αρ 0 /(2γ) is consistent with substituting (25) into (26), where g 3,a + g 3,b is now replaced by 2g 3 . Hence, it demonstrates that the rate of work is also accessible from the nonequilibrium deviation of pair correlations in fully driven liquids.
Overall, the results of this Section illustrate how dissipation affects the transport and structural properties of driven liquids, measured in terms of diffusion coefficient and density fluctuations. These findings motivate the following question: can nonequilibrium forces be tuned to reliably stabilize target configurations? To explore this, we rely in what follows on the framework of large deviation theory. In practice, our strategy amounts to biasing trajectories in terms of dissipation, related to many-body interactions by (21), to mimic the effect of an external drive. Following this route, our analytical and numerical results provide some concrete intuition for how interactions in a multicomponent system can be controllably renormalized by nonequilibrium forces. Hence, we demonstrate the ability to nucleate structures different from those characteristic of the equilibrium Boltzmann distribution to help guide self-assembly [97] and collective motion [25] far from equilibrium. These results further illustrate the interplay between energy dissipation and organization in nonequilibrium many-body settings.

III. INTERACTIONS IN BIASED ENSEMBLES
To investigate how target structures and dynamics can be promoted by means of a dynamical bias, we begin by considering a system of interacting Brownian particles without any driving force where the statistics of the noise term ξ i is the same as the one in (1). The rate of workẇ defined in (14) is zero because of the absence of driving. In Sec. II, to obtain a non-zero rate of energy flow through the system, we consider an explicit driving force F d,i and we explore its effects on the transport and structural properties of the liquid. In practice, different types of driving can lead to the same dissipation. In this section, using the framework of large deviation theory, we take an alternative approach where the dynamics is now conditioned by enforcing a required energy flow without any explicit driving. Thus, exploring how the system adapts to this requirement provides a new insight into the relation between dissipation and organization in driven systems, which is distinct from yet complements the approach in Sec. II.
To this end, we focus on the subset of noise realizations that are conditioned on a non-zero rate of work. In particular, these realizations no longer have zero average, so that one can re-define the noise term in (28) as ξ i → ξ i + F aux,i by introducing an auxiliary force F aux,i [38,50]. Hence, the stochastic dynamics given by (28) with added force F aux,i provides an explicit case which ensures a non-zero energy flow rate. In practice, this dynamics can be drastically different from the original one, thus opening the door to stabilizing unexpected structure and to promoting novel collective effects. Interestingly, such a dynamics can actually be regarded as the optimal strategy to effectively enforce a target condition on rate of work [98].
Formally, to study the dynamics conditioned by dissipation, we bias the probability of trajectories. This is done by introducing an exponential weighting factor exp κ´t 0 E(s)ds where E is the observable which conditions the dynamics, e.g. energy flow rate, and κ is a conjugate field. In practice, the relative importance of biasing in the dynamics is controlled by κ, which in turn controls the average value of E [37]. Before deriving the central results of this section, namely relations between biased energy flow rates and organization, we first introduce a simple example in which the connection between auxiliary forces and exponentially biased ensemble can be clearly seen.

A. Dynamical bias and external forces
To introduce pedagogically our methods, we first show how biasing trajectories can lead to effectively introduce a driving force. Inspired by the role of dissipation in emerging liquid properties, as discussed in Sec. II, we bias the equilibrium dynamics (28) with the sum of the dissipation and the rate of work, scaled by T , that would be produced by applying a constant force F d to a subset Ω of particles: where V = (1/2) i,j v(r i − r j ). The path probability P ∼ exp − i´t 0 A i (s)ds corresponding to this biased ensemble is obtained with standard methods [99,100]: where the two first terms correspond to the unbiased dynamics (28), and the third one to the bias in (29). It can also be written as As a result, given that the last term in (31) can be absorbed in a normalization factor, we deduce that the trajectories biased by (29) can be generated, at leading order, in a physical dynamics where the external force 2κF d is applied to every particle in Ω. In particular, it does not feature any long-range interactions which are usually found in auxiliary dynamics [101].

B. Dynamical bias and modified interactions
To go beyond the case of applying a constant force, we now seek for a dynamical bias which regulates particle interactions in a controlled manner. In particular, we examine cases where the control parameters κ ij are specific to particle pairs {i, j}, so that the biasing factor in path probability now reads exp i,j κ ij´t 0 E ij (s)ds . Now, our choice for the biasing function E ij is informed by the connection between the rate of work and many-body interactions in driven liquids, as detailed in Sec. II C. Specifically, we observe that the power balance (21) for deterministic drive (ẇ act = 0) can be written aṡ w = − i∈Ω,j Lv(r i − r j ) in terms of the evolution operator of the equilibrium dynamics (28) defined by γL = i T ∇ i −∇ i V ·∇ i . This motivates us to consider the following bias In the unbiased ensembles of Sec. II C, E ij provides a measure of the rate at which driving forces pump energy into or extract energy from the specific interaction between the i th and j th particles. Here, instead of driving the system with a specific driving force, trajectories are driven by atypical realizations of the noise generated by biased sampling.
To explore how this bias modifies interactions, we first employ a derivation different from the path integral approach in Sec. III A. Based on the procedure in [38,50], the auxiliary physical dynamics, which has the same statistical properties as in the biased ensemble, can be constructed by solving the eigenvalue equation where the eigenvalue λ, parametrized by κ ij , is the scaled cumulant generating function appropriate to E ij . The auxiliary dynamics is then defined by replacing the interaction potential in (28) by the following auxiliary potential:Ṽ In practice, computing G is a highly non-trivial procedure for many-body systems. The explicit solutions considered so far concern either exclusion processes [102][103][104] or particle-based diffusive systems restricted to small noise regimes [105,106] and non-interacting cases in some specific potentials [107][108][109].
In our case, a simple expression can be obtained for the auxiliary potentialṼ i by solving (33) perturbatively at small bias parameter κ. Specifically, we expand where G 0 is the uniform eigenvector associated with the zero eigenvalue. Given that E ij = 0 in steady state, which follows from the vanishing current condition in the unbiased dynamics ( v = 0), the leading non-trivial order of (33) reads Substituting the explicit expressions for the biasing function in (32), we then deduce that 4T G (1) ij = −G 0 v(r i −r j ) is a solution of the eigenvalue problem to order κ. The auxiliary potential follows as Therefore, biasing with (32) amounts to changing the strength of particle interaction by a factor κ ij specifically for any pair {i, j}. This is the main result of this section. While energy flows were sustained by explicit nonequilibrium forces in Sec. II, we now maintain a non-zero average for E ij by a biased sampling of trajectories. The corresponding noise realizations can be thought of as an external protocol, which leads to modifying the energy landscape sampled by the biased system as given in (37). Note that tuning interaction strength between targeted pairs is qualitatively consistent with the effect of external driving. Indeed, phase separation in mixtures of driven and undriven particles, reported both experimentally and numerically, can be rationalized in terms of an effective decrease of specific interactions between these particles [21,23].
Moreover, the techniques in Sec. III A allow one to anticipate the trajectories generated at higher-order when now biasing with (32). To this end, we consider the ensemble where the first-order dynamics, given by the potential (37), is biased with exp ´t 0 ε(s)ds defined in terms of As detailed in Appendix C, this ensemble is equivalent to biasing the original dynamics (28) with (32). Thus, the effect of higher-order bias on trajectories amounts to maximizing the squared forces in the integrand of (38), which effectively tends to cluster particles for both signs of κ ij . Finally, the decomposition between first-order auxiliary potential and higher-order symmetric bias can be extended to a generic class of biases of the form T E ij = LA(r i − r j ) for an arbitrary observable A: the corresponding first-order auxiliary potential V + 2 i,j κ ij A(r i −r j ) is now complemented with the higherorder bias (38) where A replaces v. Such a bias is reminiscent of, yet qualitatively different from, the escape rate used to promote dynamical heterogeneity in glassy systems [41,110]. In our case, clustering is favored for both positive and negative bias parameters κ ij . In particular, this is in contrast with the emergence of a hyperuniform phase, where large scale fluctuations are suppressed, reported when biasing some hydrodynamic theories of diffusive systems [111].

C. Numerical sampling of biased structures
To illustrate the potential of our bias to control liquid properties, we focus in what follows on the specific case κ ij = κδ i∈Ω δ j ∈Ω where all pairs between a subset Ω and other particles are biased with the same strength κ. Here, the set Ω could for instance refer to some tracer particles immersed in the liquid, to connect with the settings in Sec. II.
To confirm numerically the validity of our approach, we first probe the range of the first-order auxiliary dynamics where interactions are predicted to be simply renormalized. We compare measurements of i∈Ω,j ∈Ω E ij κ , where · κ denotes an average in the biased ensemble, obtained from simulations with the renormalized potentials in (37) and from a direct sampling of the biased ensemble. The latter is implemented with a cloning algorithm which regularly selects and multiplies rare realizations for efficient sampling [46,[63][64][65][66][67][68]. For convenience, interactions are now given by the soft-core potential v(r) = v 0 exp − 1/(1 − (|r|/σ) 2 ) Θ(σ − |r|). For weak interactions (v 0 = 4T ), we observe a very satisfying agreement between the two measurements for a finite range of κ, as reported in Fig. 5(a), which supports the validity of our perturbation up to interaction change between −20% and +40%. The range of validity decreases as v 0 /T increases, as shown in Fig. 5(b), and we expect a  5. (a-b) Average biasing observable i∈Ω,j ∈Ω Eij κ = i∈Ω,j ∈Ω Lv(ri − rj) κ/T as a function of bias parameter κ, where L and v respectively denote the evolution operator and the pair potential of the equilibrium dynamics (28). Results from the first-order auxiliary dynamics (solid lines) and from a direct sampling of the biased ensemble (circles) coincide for a finite range of κ. (c-d) Biased density correlation gκ as a function of interparticle distance r/σ obtained from auxiliary dynamics (solid lines) and direct sampling (dotted lines). At leading order, our dynamical bias effectively renormalizes the potential v by a factor κ for specific pairs of particles {i ∈ Ω, j ∈ Ω}, in satisfying agreement with direct sampling. This illustrates the control of liquid structure at small κ and weak interactions. Simulation details in Appendix A.
similar trend when also increasing the number of biased pairs. To explore further the features of this biased ensemble, we now compare the density correlations of biased pairs g κ (r) ∼ i∈Ω,j ∈Ω δ(r − r i + r j ) κ obtained from both direct sampling and first-order auxiliary dynamics. For κ = ±0.1, we observe that the structural modification induced by the bias becomes more dramatic as v 0 /T increases. The agreement between the cloning and auxiliary dynamics is good for the whole curve when v 0 /T = 4, whereas a deviation appears beyond r σ when v 0 /T = 12, as shown in Figs. 5(c-d). In both cases, the region of particle overlap r < σ is well reproduced. These results corroborate the ability of the first-order auxiliary dynamics to capture interaction changes as a simple renormalization of potential strength. In contrast, the tendency for particles to cluster, manifest numerically in the increased peak value at r σ, is a higher-order ef-fect missed by this auxiliary dynamics when v 0 /T = 12. Yet, note that the peak value is comparable for κ = ±0.1, in agreement with (38) being symmetric in κ. Altogether, these results demonstrate that our bias modulates the liquid structure in a controlled manner for small bias and weak interactions as predicted by (37).
Finally, we probe numerically the effect of large bias (|κ| > 1) using direct sampling, to explore configurations significantly distinct from the one of the equilibrium dynamics (28). The particles spontaneously tend to cluster for both positive and negative κ, as shown in Fig. 6 and movies in [73]. This confirms the propensity of trajectories to maximize interaction forces at high bias, as captured by (38). Importantly, the shape of clusters differs depending on the sign of κ: a micelle-like structure featuring the particles in Ω at the core (blue) surrounded by others (red) appears for κ = −3, whereas clusters have a random composition for κ = 3. Again, this agrees  6. Configurations obtained from a direct sampling of the biased ensemble where the pair interactions between red and blue particles are selectively modified. In the unbiased dynamics (κ = 0), interactions are purely repulsive with a soft core which has a similar strength for all particles, either red or blue, so that the system is homogeneous. The dynamical bias promotes clustering for both signs of κ, yet it changes interaction selectively for either sign. The repulsion is increased between red and blue particles for κ = 3, and their interactions become effectively attractive for κ = −3. As a result, the clusters which emerge spontaneously have different structures: either a random composition of mixed reds and blues (κ = 3) or a micelle-like structure with a blue core (κ = −3). This illustrates how biasing specific pairs leads to supervised spatial organization. Simulation details in Appendix A and movies in [73].
with the renormalized interactions being either increased (κ > 0) or decreased (κ < 0). In practice, the interaction strength changes sign when κ < −1 according to (37), so that the red-blue pairs are effectively attractive for κ = −3. To optimize the overall energy, the most favorable configuration then consists in maximizing (minimizing) overlap of particles in Ω (not in Ω), which in turn stabilizes a cluster of blues surrounded by reds. In general, two types of configurations should generically be stabilized for a given interaction potential v, depending on the sign of the bias. Overall, this establishes a reliable proof of principle for the design of tailored self-assembled structures with our specific choice of biased ensembles.

D. Bias-induced collective motion
As a final illustration of how collective effects can be controlled by dynamical bias, we consider a model of selfpropelled particles where interactions are now only mediated via the angular dynamics [69]: (39) where V 0 denotes the self-propulsion velocity, u(θ) = (cos θ, sin θ) is the unit vector, and µ r is the rotational mobility. The term η i is a zero-mean Gaussian white noise with correlations η i (t)η j (0) = 2D r δ ij δ(t) given in terms of the rotational diffusion coefficient D r . To promote alignment between neighboring particles, we choose the pair-wise torque as T (θ, r) = Θ(σ − |r|) sin θ/(πσ 2 ).
This dynamics was originally introduced as a generalization of the Vicsek model to continuous time [25]. Thus, it exhibits a transition between a isotropic state for small density ρ 0 and large noise D r , and a polar state for large density ρ 0 and small noise D r . In practice, the linear stability analysis of the corresponding hydrodynamic equations predicts the transition to occur when 2D r = µ r ρ 0 [69]. In what follows, our aim is to show that such a transition can also be mediated by tuning interactions with a dynamical bias.
To this end, we take the biasing factor in path probability as exp κ´t 0 E θ (s)ds], where the biasing observable E θ reads The average value E θ is proportional to d dt i,j cos(θ i − θ j )Θ(σ−|r i −r j |) , which vanishes in steady state. Then, following the procedure detailed in Sec. III B, we deduce that biasing the dynamics (39) with E θ amounts to considering renormalized interactions of the form Thus, by promoting a non-zero average for E θ , aligning interactions can be tuned in a controlled manner at first order in κ. Higher order bias leads to maximize the squared torque, as presented in Appendix C. We test this prediction numerically using a direct sampling of biased trajectories. We consider values of Isotropic state Polar state Configurations obtained from a direct sampling of the biased ensemble for aligning self-propelled particles. The color code refers to the orientation of particles. In the unbiased dynamics (κ = 0), we observe isotropic and polar states respectively at large noise (Dr > D * r ) and small noise (Dr < D * r ). Here, the critical noise is D * r = 8 and we take the noise values Dr = {7, 9} for respectively the polar and isotropic regimes. The dynamical bias leads to renormalize interactions in a controlled manner, which effectively changes the transition threshold as D * r → D * r (1 + κ) at leading order. As a result, one can stabilize either isotropic or polar states, respectively for κ < 0 and κ > 0, thus illustrating the ability to trigger or inhibit collective effects in nonequilibrium systems. Simulation details in Appendix A and movies in [73].
{D r , µ r , ρ 0 } above and below the threshold D * r = µ r ρ 0 /2, where the system exhibits either isotropic or polar states in the unbiased dynamics (κ = 0), as shown in Fig. 7. Specifically, when the original system is isotropic (D r > D * r ), we observe a transition to polar for κ > 0, and, conversely, when it is polar (D r < D * r ) a transition to isotropic for κ < 0. This confirms our result (41) where the bias amounts to changing the angular mobility as µ r = (1 + κ)µ r + O(κ 2 ) for weak κ, so that the linear instability is either triggered or suppressed by solely tuning κ, all other parameters being held the same. Thus, these results demonstrate how biasing the dynamics with an appropriate observable leads to control the emergence of spontaneous organization, with potential interest for other nonequilibrium dynamics.

IV. CONCLUSION
Developing techniques to characterize and control the behavior of systems operating far from equilibrium remains a central and outstanding problem. Despite the apparently complex interplay between internal dissipation and emerging properties, we have demonstrated that tracer diffusion and density correlations can simply be connected to dissipation in driven liquids. We have also constructed a mapping between deterministic and active drives for a specific active matter model, thus showing how our approach can potentially be extended to a broad class of systems. Importantly, our results open promising perspectives to evaluate dissipation simply from the structure of the system. Inspired by recent works [18,112], one could also introduce a map of dis-sipation, directly related to the statistics of interaction forces, to resolve spatially where energy is released in the thermostat. Though the corresponding integrated map would not cover the total dissipation, it would already provide insightful information about locations of low and high dissipation with respect to a constant background set by the squared driving amplitude.
In practice, monitoring dissipation with a well-defined parameter remains an open challenge for many-body systems. To this end, biased ensembles enable one to specify the statistics of dissipation by introducing an additional control parameter, analogously to the change from microcanonical to canonical ensemble in equilibrium thermodynamics [38,50]. This is done by selecting rare noise realizations which drive the system away from typical behaviour, without introducing any driving force. Pioneering works were focused on favoring dynamical heterogeneities, without affecting the structure, of kinetically constrained models [39][40][41][42][43]. Yet, more recent studies have shown the potential to also modify density correlations in diffusive systems [48,49,113].
Using these large-deviation techniques, we have put forward a particular set of biased ensembles which allows one to regulate the liquid structure in a controlled manner. The explicit form of the bias is motivated by the relations between dissipation and structure that we have derived for driven liquids. At leading order, any bias in this class simply leads to introducing additional interactions in the dynamics. Furthermore, higher-order bias systematically constrains the trajectories to favor the formation of clusters. Based on minimal case studies, we have sampled the biased configurations, using stateof-the-art numerics [46,[63][64][65][66][67][68], to illustrate the ability to stabilize specific structures and collective effects in a controlled manner.
Since dynamical bias consists in favoring rare noise fluctuations, the corresponding dynamics effectively provides useful insights on how to promote atypical configurations with an external drive. In practice, the driving protocol should simply mimic the biased noise realizations. This line of thought has already been exploited for efficient sampling of the biased ensemble [45,65,114,115], where control forces make rare events become typical. Moreover, since our analytic framework encompasses the case of a specific bias for each pair of particle, it could potentially be regarded as a fruitful route to promote the spontaneous self-assembly of complex structures at the cost of energy dissipation. For instance, inspired by recent works [116,117], one might consider our approach to design energetic landscapes, in terms of the pair-specific bias parameters, which selectively stabilize some target molecules.
Overall, these results illustrate how specifying the amount of energy dissipated by nonequilibrium forces allows one to constrain the dynamics and structure of driven liquids. This paves the way towards controlling the emerging properties of such systems by tuning dissipation accurately. It remains to investigate whether similar results can be obtained in more complex systems which could, for instance, potentially include anisotropic building blocks, such as driven chiral objects or active liquid crystals [27,118,119] In Sec. II A, a custom code of molecular dynamics, based on finite time difference, is used to perform the simulations in a two-dimensional box 10 2 σ × 10 2 σ with periodic boundary conditions. The time step is δt = 10 −4 and the initial condition is homogeneous. Parameter values: ρ 0 = 0.7, T = 0, γ = 1, f = 3 × 10 −2 , τ = 10 3 , v 0 = 5, n = 10 2 , σ = 1.
In Secs. II B and II C, numerical simulations of the dynamics (1) are performed using the LAMMPS simulation package in a two-dimensional box 10 2 σ ×10 2 σ, where σ is the particle diameter, with periodic boundary conditions at average density ρ 0 = 0.45. Our custom code actually implements overdamped Langevin equations of motion with finite time difference. It simply utilizes the efficient force computation routines that are built as a part of the Molecular dynamics package. The system is first relaxed for 10 3 conjugate gradient descent steps, and later equilibrated during 50τ . We evaluated average values over 10 independent trajectories with duration 150τ . The density pair correlations were constructed using 10 independent trajectories each with duration 50τ . We performed error analysis from the independent simulations and obtained negligible errors for all the data in Figs. 2-3-4. The time step is 5 × 10 −4 and the bin size for computing the pair correlations is 0.01σ. We performed simulations at other values of the time step {10 −4 , 10 −5 } and of the bin size {5 × 10 −3 σ, 2 × 10 −2 σ} to confirm that our calcula-tions of the α coefficients are well converged. Parameter values: T = 1, γ = 10 2 , v 0 = 1, σ = 1.
In Sec. III C, a custom code of molecular dynamics, based on finite time difference, is used to perform the simulations in a two-dimensional box 10σ × 10σ with periodic boundary conditions. We bias the pair potential between 8 blue particles and 16 red particles. To sample the biased ensemble, we use the cloning algorithm described in Appendix A of [65]. The time interval for cloning is ∆t = 10δt and the number of clones is 1600. The time step is δt = 10 −4 , the initial relaxation time is 10 4 ∆t, and the total simulation time is 10 6 ∆t. Parameter values: T = 1, γ = 1, v 0 = 4 (Fig. 6), σ = 1.
In Sec. III D, a custom code of molecular dynamics based on finite time difference is used to perform the simulations. N = 128 particles are simulated in a two dimensional box of size 4σ × 4σ with periodic boundary conditions. To sample the biased ensemble, we use the cloning algorithm described in Appendix A of [65]. The time interval for cloning is ∆t = 10δt and the number of clones is 200. The time step is δt = 10 −3 , and the total simulation time is 500∆t. Parameter values: V 0 = 2, µ r = 2, ρ 0 = 8, σ = 1.

Appendix B: Dissipation and diffusion
This appendix is devoted to the derivation of the dissipation rate J and the diffusion coefficient D of a driven tracer, as defined in Sec. II. We employ a perturbative treatment at weak interactions, originally introduced for a particle driven at constant force in [52,53]. To this aim, the tracer-bath interaction potential v is scaled by a small dimensionless parameter h 1 in what follows. Besides, we focus on the regime of dilute tracers, so that interactions among them, either direct or mediated by the bath, can be safely neglected.
The dynamic action associated with the tracer dynamics (11)(12) follows from standard path integral methods [99,100]. It can be separated into contributions from the free tracer motion and from interactions, respectively denoted by A 0 and A int : where D 0 = T /γ is the tracer diffusion coefficient in the absence of interactions (v = 0), andr 0 is the process conjugated with the tracer position r 0 . For weak interactions h 1, any average value can be then expanded in terms of h as · = · 0 −h 2 A int · 0 +O(h 4 ), where · 0 is the average taken with respect to A 0 only. As a result, determining the first correction from interactions in any observable amounts to computing the corresponding average A int · 0 .
Considering the dissipation rate per particle J /N = ṙ 0 ·F d , the leading order is ṙ 0 0 ·F d = |F d | 2 /γ = f 2 /γ, and the first correction reads −h 2 A intṙ0 0 ·F d . Given the explicit form of A int in (B1), the correlations of interest are where we have used that the tracer statistics is Gaussian in the absence of interactions, following [52,53]. From this result, we get where we have used γ G = γ/ρ 0 and D G = ρ 0 D 0 . Expanding at small f , we deduce (B4) Substituting the explicit expression of the deterministic drive (2) in (B4), and then integrating over u and w, we obtain (B5) For the case of active drive with correlations (3), we exploit the equivalence with a disordered drive detailed in Sec. II A. Substituting the explicit drive (4) in (B4), and then averaging over disorder in the limit of many oscilla-tors (n 1), we get The asymptotic results for the rate of workẇ = f 2 /γ−J , presented in Sec. II B for both deterministic and active drives, follow directly.
We now turn to deriving the diffusion coefficient D. It is defined in terms of the mean-squared displacement (MSD) ∆r 2 0 (t) = r 0 (t) − r 0 (t) 2 as D = lim t→∞ ∆r 2 0 (t) /2dt. At leading order, the MSD reads ∆r 2 0 (t) 0 = 2dD 0 t. To obtain the first order, we need to compute the following correlations where we have used again that A 0 is Gaussian in terms ofr 0 , yielding Ensemble (b) is associated with the first-order auxiliary dynamics, whose potential reads V +2 i,j κ ij A(r i −r j ), γṙ and k in (C4) and comparing with A (b) k in (C3), it appears that A (a) and A (b) are indeed equal up to a boundary term proportional to i,j κ ij A(r i (t)−r j (t))− A(r i (0) − r j (0)) which can be neglected at large t: this establishes the equivalence between ensembles (a) and (b).
We now turn to demonstrate the equivalence between two ensembles related to the Vicsek-like dynamics (39). Ensemble (c) is biased with the factor exp κ´t 0 E θ (s)ds , where E θ is defined in (40). Ensemble (d) corresponds to the first-order auxiliary dynamics, with torque given by (1 + κ)T , biased with exp ´t 0 ε θ (s)ds where The dynamic actions for each ensemble, denoted by B (σ) (t) = i´t 0 B and Expanding B in (C7), it appears that B (c) and B (d) only differ by a term proportional to i,jθ i T (θ i − θ j , r i − r j ) which can safely be neglected at large t, thus proving the equivalence between ensembles (c) and (d).