Shear-induced deconfinement of hard disks

Using Brownian dynamics simulations, we investigate the response to shear of a two-dimensional system of quasi-hard disks that are confined in the velocity gradient direction by a smooth external potential. Shearing the confined system leads to a homogenization of the one-body density profile. In order to rationalize this deconfinement effect, we split the internal onebody force field into adiabatic and superadiabatic contributions. We demonstrate that the superadiabatic force field consists of viscous and of structural contributions. We give an empirical scaling law that yields results for the superadiabatic force profiles both in the flow and in the gradient direction, in excellent agreement with the simulation data.


Introduction
Applying shear to a complex substance constitutes a means to drive the system out of equilibrium in a well-controlled way, and results in arguably one of the most fundamental nonequilibrium setups [1,2]. In its simplest form, shearing is characterized by the shear rateγ as the only relevant (and scalar) nonequilibrium parameter. The physics at play, even only in steady state, is however fundamentally different from the equilibrium properties of the same material. Hence, shearing is an excellent model situation for systematically studying soft matter out of equilibrium.
For the important material class of colloidal dispersions, Matthias Ballauff and collaborators have performed sterling work, developing and exploiting ingeniously tailored particles that respond to temperature variation. Despite the quite complex internal core-shell structure of these thermosensitive colloids [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], the particles interact via an essentially short-ranged steeply repulsive pair potential. Changing the temperature facilitates systematically changing the effective particle size and hence the typical length scale of the interparticle interactions. The underlying mechanism is the thermoresponse of the polymeric particle shell. Controlling the particle size allows to control accordingly the effective Matthias Schmidt Matthias.Schmidt@uni-bayreuth.de 1 Theoretische Physik II, Physikalisches Institut, Universität Bayreuth, 95447, Bayreuth, Germany colloidal packing fraction in the system. As the response to shear depends very sensitively on the packing fraction, thermosensitive colloids give direct access to this crucial thermodynamic parameter. The shear rateγ , often expressed as a dimensionless Peclet number Pe [1], is an external parameter that controls the degree of nonequilibrium driving that the system is exposed to.
Mode-coupling theory (MCT), as spearheaded by Fuchs and coworkers [5, 6, 8-12, 15-17, 20], has provided a platform for rationalizing and in many cases quantitatively describing the results from experiments, as obtained rheometrically. Much insightful work is based on the "schematic model" [5,6,11,15,17] of MCT. The impressive degree of consistency, in qualitative and in quantitative terms, for important quantities is achieved with a very moderate level of empirical input. This is remarkable, as the considered quantities vary typically over many orders of magnitude. MCT operates on the level of two-body correlation functions. In particular, via Green-Kubo relations, the equilibrium stress-stress autocorrelator is set in relation to the nonequilibrium response under shear. The full theory is hence a complex one, which justifies the need for having the schematic model, which is regarded as describing the essential features of the dynamics of a generic time correlation function, in particular historydependence via a memory integral over previous times. The sole equilibrium input in MCT is the static structure factor S(q) of the fluid. The entirety of the time-dependent nonequilibrium phenomena that occur under driving arises from the dynamical structure of the MCT equations of motion.
On principal grounds, one could expect that additional equilibrium information, besides S(q), could be required or be at least useful. Brader, Krüger, and their coworkers [21][22][23][24][25] have hence gone beyond the static limitation by incorporating ideas from classical density functional theory (DFT) for inhomogeneous fluids [26,27]. DFT is a framework that is genuinely adapted for and capable of describing situations where the locally resolved microscopic density distribution ρ(r) is inhomogeneous in space: ρ(r) = const, where r indicates position. Two important (and useful) relationships for equilibrium systems are the following: (i) The internal contribution to the one-body force profile, f ad (r), is solely dependent on the density profile, but it is independent of the external forces that act in the system (say due to the presence of walls). Hence, f ad (r) = f ad (r, [ρ]), where the brackets indicate a functional relationship of mapping the entire function ρ(r ) at all space points r to the force at the given position r. The internal force field f ad (r) arises from all interparticle forces that the remaining particles exert on a particle located at position r in an equilibrium situation. (ii) Two-body correlation functions are contained in the approach and are accessible, alternatively, by the Ornstein-Zernike or the test particle routes (see Ref. [28] for a recent account of the virtues of the test particle approach when combined with an approximate free-energy functional).
Krüger and Brader's work [21][22][23][24][25] for systems under shear is based on the so-called dynamical DFT [26,29,30] for nonequilibrium dynamics. Nevertheless, the approach described in Refs. [21][22][23][24][25] still genuinely works on the two-body level, which a priori puts high strain both on the physical intuition that is required to devise approximation schemes, as well as on computational demand. Although dynamical DFT operates on the one-body level, for a sheared system with flow direction being orthogonal to the density gradient, the bare DFT [26,29,30] gives a null result, in contrast to Krüger and Brader's more sophisticated theory [21][22][23][24][25]. Dynamical DFT, in its bare version [26,29,30], is hence defunct for the description of arguably the most basic nonequilibrium situation of a colloidal system, or more generally, of an inhomogeneous classical liquid. The reason for the failure is of fundamental nature: the true nonequilibrium dynamics is described as consisting of a sequence of "adiabatic states" that are taken to be at equilibrium. Clearly, this assumption is not true in general, and for the case of shear, it is in quite striking contradiction to reality.
However, there is hope for formulating a complete description on the one-body level, as the superadiabatic contributions that occur above the adiabatic effects (correctly accounted for in the dynamical DFT) are both well-defined and well-characterizable objects [31,32] from an extended, kinematic functional point of view. Here, the microscopically resolved velocity profile v(r, t) is a variable on par with the time-dependent density profile ρ(r, t), where t indicates time. The full nonequilibrium dynamics is hence driven both by adiabatic effects, which are functionally dependent only on the instantaneous density distribution, and superadiabatic effects, which possess kinematic dependence on both ρ(r, t ) and v(r, t ) for times t ≤ t, where t is the time of interest. Hence, the superadiabatic contribution is nonlocal both in space and in time, while the adiabatic contribution is Markovian (instantaneous) in time, but nonlocal in space.
For the (important) case of the internal one-body force field f int (r, t), which is generated from the underlying interparticle interaction potential, the adiabatic-superadiabatic splitting is This formal result [31] has implications that are important, practical, and testable. In particular, the independence of f sup (r, t) from the external forces allows to disentangle intrinsic from external effects and hence it offers the potential for deep insights into the coupled nature of the nonequilibrium many-body physics [31]. Power functional theory (PFT) not only provides the existence of the functional map (ρ, v) → f sup . It also establishes on a microscopic footing a rigorous minimization principle for the velocity profile, or equivalently for the microscopic one-body current J(r, t) = v(r, t)ρ(r, t), implying vanishing functional derivative at the minimum, δR t [ρ, J]/δJ(r, t) = 0. Here, R t is the (total) free power functional, which consists of a sum of ideal, external, adiabatic, and superadiabatic contributions. The latter is the superadiabatic excess (i.e., over ideal) free power functional P exc t . The superadiabatic internal force profile is obtained as a functional derivative f sup (r, t) = −δP exc t /δJ(r, t). As in equilibrium DFT, the adiabatic internal force profile is obtained via f ad (r, t) = −∇δF exc /δρ(r, t), where F exc [ρ] is the intrinsic excess Helmholtz free-energy functional, and nabla denotes the derivative with respect to position r. PFT has been both applied and extended over a growing range of physical systems and situations . We give a brief overview in the following.
Here, we present the first study of superadiabatic forces in a system under simple shear, i.e., flow that is characterized by linear dependence of the velocity on position, and hence spatially constant imposed shear ratė γ . In order to trigger a response of the system on the one-body level, as befits the concepts underlying the PFT framework, we expose the sheared two-dimensional system to an additional confining external potential. We choose the (arguably) simplest possible geometry, where the external potential V ext (y) varies only in the shear gradient direction y. As a consequence, all one-body quantities become functions of (only) the y-coordinate, and the system remains translationally invariant in the flow direction x.
The above setup is interesting, as it allows the study of the influence of shear on a prototypical confined situation for a simple fluid, i.e., that of a periodically varying, oscillating potential (along the y-direction of the shear gradient). We can hence monitor both the structural forces that act in the y-direction and the viscous forces that act in the flow direction x. We present simple phenomenological scaling laws that fit the results quantitatively. We consider moderate fluid densities only and hence deliberately stay away from the intricacies of the glass transition. We do, however, consider cases of high driving, where we find the system to exhibit interesting and apparently universal saturation behavior. The saturation relates the external conservative force −∇V ext and the intrinsic superadiabatic force to each other. This is a remarkable mechanism, as external and superadiabatic forces are very different in origin. As a result of the shear, the density profile homogenizes, and hence, shear acts against the confinement induced by V ext . This paper is organized as follows: In "Description of the system," the physical system considered is described and the Brownian dynamics (BD) simulation algorithm is laid out. The results of the simulations are presented in "Results." In "Conclusions," we conclude and give an outlook on possible future work.

Microscopic dynamics
Our two-dimensional system consists of N (indistinguishable) circular particles suspended in an incompressible implicit solvent. The particles are quasi-hard disks; their pair interaction potential is given by where r is the center-center distance between two particles, σ = const denotes the diameter of the particles, and = const sets the energy scale.
In addition to their internal interactions, the particles are subject to thermal fluctuations and to an external force field. The thermal fluctuations induce a random force f ran i (t) on particle i with the following statistical properties: where the angles denote an average over different realizations of the noise, t and t indicate time arguments, ξ denotes the friction coefficient, k B indicates the Boltzmann constant, T indicates the temperature, δ(·) represents the Dirac delta distribution, is the 2 × 2 unit matrix, and the integers i, j = 1 . . . N label the particles.
The external force field splits according to where y i indicates the y-coordinate of particle i and ∇ i is the derivative with respect to the position r i of particle i. Here, the first term is a non-conservative linear shear field acting in the x-direction and the second term is a conservative compression force along the y-direction, given, respectively, by where the shear rateγ is an inverse timescale that characterizes the strength of the driving andê x is the unit vector in the x-direction. The confining potential (6) compresses the system symmetrically in the y-direction towards y = L/2; the constant V ext 0 sets the depth of the potential well. The maximum is at y = 0, and minimum of the well is at y Figure 1 shows an illustration of the model. We use periodic boundary conditions in both spatial directions; L denotes the size of the (square) simulation box. Fig. 1 Sketch of the geometry of the system. Shown are the external potential V ext (y) that compresses the system towards y = 0 (dasheddotted red line); the resulting conservative force field −∇V ext acting in the y-direction (red arrows indicate direction and the solid red line indicates the magnitude); and non-conservative force field (blue arrows) f shear due to the externally imposed shear flow. The system is translationally invariant in the x-direction. Lees-Edwards boundary conditions render the shear force continuous in the periodic images of the system along the y-direction We consider overdamped dynamics, where the inertia of the particles can be neglected and the dynamics of the system is described by the Langevin equation of motion [27]: withṙ i denoting the time derivative of r i .

One-body correlation functions
In steady state, the system is invariant under translation in the x-direction and all one-body quantities (as described below) only depend on the y-coordinate. In our simulations, we check that a steady state is reached before sampling any quantities of interest.
As the one-body density, will not vary in time in steady state, its time derivative ∂ρ(r)/∂t = 0, and the continuity equation reduces to implying that the particle current J(y) is also constant in time. Here, the one-body current is defined as where the velocity v i of particle i at time t is given by a centered difference of its position vector [45], where t is the time step of the numerical integration routine (as detailed in "Brownian dynamics simulations" below). See Ref. [45] for the derivation of the finite difference expression (12) for the velocity in Brownian dynamics. We split the current into three distinct contributions, corresponding to the force densities due to internal interactions (F int ), thermal diffusion (−k B T ∇ρ), and external influence (ρf ext ), and hence, Here, the one-body internal force density distribution is defined as where f int i , as given by Eq. 8, is the force acting on particle i due to the internal interactions with all other particles in the system. The force fields are obtained by dividing force densities by the density profile, i.e., f int = F int /ρ and f ran = −k B T (∇ρ)/ρ ≡ −k B T ∇ ln ρ. In the results described below, we present data for force fields rather than for force densities. The thermal fluctuations generate an average force caused by the thermal fluctuations, Splitting the current and therefore the forces is useful as it enables us to perform the splitting of the internal force field into its adiabatic and super adiabatic contribution. The velocity profile is then obtained as v(y) = J(y)/ρ(y).

Adiabatic construction
The adiabatic system is defined as having the same equilibrium one-body density as the nonequilibrium system [31]. The Mermin-Evans theorem, which is at the heart of DFT, states that in equilibrium for a system with given internal interaction potential, at temperature T , volume V and chemical potential μ fixed, there is a unique mapping from a given one-body density distribution ρ(r) to a corresponding external potential [26,27]. In order to construct the adiabatic system, we use two different methods which we show below to give consistent results. The adiabatic construction implies to take the nonequilibrium density as an equilibrium density and to calculate the corresponding external potential V ad (y) = −μ ad (y), where μ ad denotes the intrinsic chemical potential in the adiabatic system. We use the scaled-particle theory (SPT) equation of state for the hard disk fluid [27], which implies the following form of the chemical potential: Here, η * = π/4σ 2 BH ρ, where σ BH denotes the Barker-Henderson diameter [27], given by σ BH = ∞ 0 dr(1 − e −φ(r)/(k B T ) ). The thermal wavelength only adds an irrelevant constant to the chemical potential. Using Eq. 18 is the first way we construct the adiabatic system. The notation V SPT ad (y) is used to denote external potentials calculated this way. Note that the dependence of the right hand side of Eq. 18 on y is explicitly known, as it arises from the known density profile ρ(y), taken as an input.
The second way of performing the adiabatic construction also employs an approach based on a local density approximation (LDA), but using simulation data as input. Using the LDA is reasonable for the cases considered in the results section below, as the density profile ρ varies slowly on the scale of the particle size σ . For given nonequilibrium conditions, we base the LDA on simulation data of an equilibrium system (i.e., forγ = 0) with an unchanged value of V ext 0 as compared to the nonequilibrium system. This yields, in an LDA approximation, the bulk chemical potential of the system as a function of density, as both the (imposed) external potential V ext (y) and the (sampled) density profile ρ eq (y) are known as a function of position y. Eliminating the parameter y then yields the desired form of the chemical potential as a function of density. In more detail, let y eq (ρ) be the inverse function of ρ eq (y). We then obtain the (bulk) chemical potential as a function of (bulk) density as Using the equation of state Eq. 19 and the results for the nonequilibrium density profile ρ(y), we can now calculate the external potential in the adiabatic system via V LDA ad (y) = −μ LDA ad (ρ(y)).

Brownian dynamics simulations
To simulate the system described in "Description of the system," Brownian dynamics simulations are used. These are based on the integration of Eq. 7, which can be performed with the Euler algorithm, Here, t is the time step used in the algorithm and δr i is a random displacement with Gaussian distribution of standard deviation √ 2 tk B T /ξ. We use the Box-Muller transformation to generate the Gaussian distribution from random numbers generated by the Mersenne Twister algorithm.
At the start of each simulation run, N particles are randomly distributed in a square simulation box of size V 2D = Nπσ 2 /η, while ensuring that there are no overlaps between any pair of the particles. The side length of the square box is L = √ V 2D . We use Lees-Edwards (sliding box) boundary conditions [58,59].
All our simulations were performed with N = 256 particles and with a time step of t = 5 · 10 −6 τ 0 where τ 0 denotes the reduced time scale τ 0 = ξσ 2 / . The constants , ξ , and σ are set to unity; we only consider the case of the reduced temperature T * = k B T / = 1. The packing fraction η, the shear rateγ , and strength of the confinement V ext 0 are our control parameters. When the Peclet number Pe is defined using the particle diameter (rather than the particle radius) as the unit of length, the value of our dimensionless shear rate is identical to that of the Peclet number, as Pe ≡γ σ 2 /D =γ τ 0 /T * , where the diffusion constant is defined as D = k B T /ξ.
The one-body density ρ(y) and the average instantaneous force due to interparticle interaction f int (y) are sampled. The histograms used for the sampling of these parameters have a resolution of σ/15 but are down sampled to a resolution of σ/5 to reduce the noise. For each presented set of parameters, multiple simulations were run (typically at least 150). The results of these are then averaged to produce the final data. Each of the simulations is at first run for t eq = 0.75 τ 0 , so the system can relax into the steady state. After t eq has passed, the sampling of the physical quantities of interest is started at a rate of 2000 τ −1 0 . The simulation is then kept running for t run = 200 τ 0 .
In order to test whether the steady state is reached after t eq , four sets of parameters have additionally been run with t eq = 100 τ 0 and t run = 100 τ 0 and were then compared to the results of the density profiles produced by simulations with t eq and t run . The results of these runs show that the difference between the datasets is smaller than 1% and shows no bias, indicating that t eq = 0.75 τ 0 is sufficient for the relaxation into the steady state. For smaller values of the shear rate than considered here, using longer relaxation times t eq might be necessary.

Behavior of the density profile
We first present the effects that varying the shear rateγ and the strength of the compression V ext 0 has on the density profile ρ(y) of the system. We have considered moderate values of packing fraction ranging from η = 0.1 to 0.35. As a reference, the hard disk fluid transitions to a glassy state for packing fractions larger than η = 0.699 [60]. In equilibrium, the hard disk fluid is stable for packing fraction below 0.701 ≡ (π/4)0.892 (cf. Ref. [61]). The presence of the confining potential leads to a local increase in density, which potentially could increase the local density above the transition packing fraction for higher values of η. As different sets of the parameters only appear to lead to quantitative changes of the density profile and not to different physical effects, results are presented for selected values of parameters.
In Fig. 2, results of simulations with packing fractions of 0.2 and 0.3 are presented. Starting with η = 0.2 in panel (a), we show the effects of an increase in the strength of the confinement on the density profile, while the shear rate is kept constant. We observe that the density is compressed into the trough of the conservative potential V ext and that an increase of V ext 0 leads to a corresponding increase in the amplitude of the peak in the density profile. Figure 2c shows how the density profile at constant strength of confinement is altered by the shear rate. It is apparent that shearing the system leads to a pronounced flattening of the density profile. Figure 2b and d are shown to evaluate the effect that a change in packing fraction has on the density profiles. One can observe that the magnitude of the deviation of the local density from the bulk density ρ b = N/V 2D for corresponding sets of parameters remains similar.
The effect of variation of a single parameter on the density profile is indicated by an arrow. The symbols denote density profiles ρ LDA ad (y) (crosses) and ρ SPT ad (y) (circles), obtained by the two methods of adiabatic construction described in Section "Adiabatic construction"

Forces in the nonequilibrium system
We next present results for the internal force field f int (y). This is a particularly interesting quantity, as it contains the superadiabatic part of the force in the nonequilibrium system (cf. Eq. 1). We first investigate the y-component, f int y (which acts in the gradient direction). The main panels in Fig. 3 display the results for the same control parameters as the density profiles shown in Fig. 2. Figure 3a and b illustrate the relationship between the degree of confinement of the system and the behavior of the internal force field, in steady state. As expected, Fig. 3a shows that the magnitude of f int y increases as the applied external potential is increased. One can also observe that in the homogeneous sheared system (i.e., without compression, V ext 0 = 0 and therefore with constant density profile), f int y vanishes. The effects that an increase in shear rate has on f int y are presented in Fig. 3c and d. The data shown suggests that the shear rate increases the magnitude of f int y very rapidly for low shear rates, but f int y appears to reach a saturated state foṙ γ τ 0 → ∞. This is intuitively clear considering that f int y is unlikely to be larger then −f ext y in steady state. Recall that V ext is kept constant in Fig. 3c and d. The insets in Fig. 3 show the y-components of the force fields decomposed as in Eq. 14, along with the total force field f total y , i.e., the sum (14) divided by the local density. We show results for four different systems. We expect f total y (y) to vanish because of the symmetries in the system. As throughout, subscripts x and y denote the components of vectors. In Fig. 3 e, f, g and h, it is apparent that f total y indeed vanishes identically, apart from some residual noise in the data. The figure also demonstrates that both f int y and f ran y act in the opposite direction to f ext y (y) and their sum exactly cancels the applied compression.
Comparing the results of the equilibrium system shown in Fig. 3h to the sheared system in Fig. 3f reveals a pronounced difference in the way the force, due to the compression f ext y , is compensated. In Fig. 3h, f int y is not too different in magnitude to f ran y , while in the nonequilibrium case presented in Fig. 3f, the entropic force f ran y is much smaller and f int y almost solely compensates f ext y . This is in line with the results from the density profiles, as the homogenization unavoidably leads to the Fig. 3 Internal force field f int y (y) as a function of y/σ for multiple sets of parameters η,γ , and V ext 0 . The parameters held constant in panels a-d are the same as in Fig. 2. The insets show contributions to the force balance in Eq. 14 for four examples. The internal and the total force are sampled in simulations, while the deterministic external force is imposed and the random force is calculated using Eq. 16  Obviously, this effect becomes more apparent as the shear rate and therefore the homogenization is increased, as is clear from a comparison of Fig. 3g and e. In Fig. 4, the x-component of the internal force field, f int x (i.e., the component that acts in the flow direction) is plotted for the same parameters as f int y above. We leave out the results for a shear rate ofγ τ 0 = 300, as this is very noisy. We first consider the effects of an increase of V ext 0 on the density profile (cf. panels (a) and (b)). One can observe that the relation between V ext 0 and f int x is practically linear. In Fig. 4 c and d, the effects of a variation of the shear rate on f int x are illustrated. The magnitude of the force increases as the shear rate is increased. However, the force does not increase linearly.

Adiabatic construction
We use the two methods outlined in "Adiabatic construction" to perform the adiabatic construction. In order to assess their performance, we have carried out simulations of the adiabatic system, i.e., runs in equilibrium under the influence of an external potential V ad (y), as obtained following the two procedures (SPT and LDA) described in "Adiabatic construction." In the following, we compare the respectively obtained density profiles ρ SPT ad (y) and ρ LDA ad (y) to each other, as well as to the corresponding (target) nonequilibrium density profile ρ(y).
The adiabatic density profiles for representative parameter choices are shown as symbols in Fig. 2. The quality of the agreement demonstrates that both versions of the adiabatic construction reproduce the density profile of the nonequilibrium systems very well. The relative difference in density is smaller than ∼ 1%. Comparing the different methods of adiabatic construction reveals good agreement between the results, showing that our simulations are able to reproduce the behavior of the hard disk fluid as predicted by SPT. However, the approach using only our simulation data and LDA to calculate the chemical potential is slightly superior. Therefore, we use data obtained with these simulations in the analysis to obtain the superadiabatic force field.

Superadiabatic forces
Having obtained results for both the internal force field in the nonequilibrium system, f int , as well as in the adiabatic system, f ad , we proceed by using the force splitting (1) in order to calculate the superadiabatic forces according to As the adiabatic system is translationally invariant in x, we expect the corresponding force component to vanish, f ad x ≡ 0. Our simulation data confirms this expectation. Equation 21 accordingly implies f sup x (y) = f int x (y); hence, in the flow direction, the superadiabatic force is identical Fig. 4 Viscous force field, scaled as f int x (y)σ/ , as a function of the scaled position coordinate y/σ for multiple sets of parameters η,γ , and V ext 0 . The parameters held constant are the same as in Fig. 2. In panels (c) and (d), the circles denote the superadiabatic force calculated using Eq. 22 and the density profiles shown in Fig. 2 to the full internal force field. These results were discussed above in Section "Forces in the nonequilibrium system". The situation in the gradient direction is different, as the density inhomogeneity in y leads to f ad y = 0. In Fig. 5, the y-component of the superadiabatic force field is presented in the main panels. The superadiabatic force behaves in a similar manner as the internal force: It reaches a saturated state as the shear rateγ is increased and it is proportional to the strength of the external confinement.
Comparing the typical magnitude of the superadiabatic force (shown in Fig. 5) to that of the total internal force (shown in Fig. 3) reveals that even for small shear rates f sup y already makes up about one-tenth of f int y . The relative contribution of f sup y increases as the shear rate is increased and surpasses the adiabatic contribution for the intermediate values of shear rate considered. For the highest values ofγ τ 0 presented, the internal force is completely dominated by superadiabatic contributions. We conclude that the superadiabatic effects in sheared systems play a significant role and that any attempt at describing the nonequilibrium dynamics requires an understanding of the superadiabatic contribution.
The saturated state corresponds to the state of the system with vanishing density gradient. Comparing both components of f sup (y) with each other reveals they have a similar shape, but with the sign flipped and that the xcomponent has a larger magnitude. This again implies a direct link to the density gradient, which we investigate in the following section.

Scaling of the superadiabatic force fields
The observations described in "Superadiabatic forces" suggest that the superadiabatic force profile depends linearly on the gradient of the density. It is furthermore apparent that this force also depends on the applied shear rate. Assuming a linear dependence on the shear rate, the simplest possible representation of the superadiabatic force as a scaling law is where the constant c y is dependent on the packing fraction η and on V ext 0 . In order to compare this empirical rule to our data, we need to take the gradient of the densities presented in Fig. 2. As direct numerical differentiation produces very noisy results, we first fit a polynomial of order 4 to the density profiles and then carry out the differentiation analytically.
The results of this calculation (for the y-component) are denoted by the circles in Fig. 5 a 1 and b 1 . We find that for low shear rates up to ≈ 100γ τ 0 , the y-component of the scaling law shows very good agreement with our data for both considered values of the packing fractions. However, the scaling form overestimates the superadiabatic force in case of high shear rates, which implies that the saturation effects discussed above are not captured perfectly.
An alternative way of testing the form Eq. 22 can be obtained by considering both the force balance in the nonequilibrium system, where the prime denotes the derivative with respect to y.
Here, Eq. 23 is obtained by dividing the force density balance (13) by the density profile and using the adiabaticsuperadiabatic splitting (1) to express the internal force field f int as a sum of adiabatic and superadiabatic force contributions. Observing that the y-component of the current vanishes in the sheared steady state then yields (23). Equation 24 is the analogue in the adiabatic system; as the adiabatic state is in equilibrium, there are neither flow nor superadiabatic contributions by construction. Using Eq. 24 to eliminate f ad y from Eq. 23, using the scaling form Eq. 22, and integrating in y, one arrives aṫ where d is an integration constant and c y is the same constant as used in Eq. 22. This balance equation only needs the nonequilibrium density and the adiabatic potential (as obtained in the simulations) as input to carry out a test of the scaling law for the superadiabatic forces. The corresponding results are depicted in panels a 2 and b 2 of Fig. 5. One observes the same level of agreement as in the case of the direct comparison of the superadiabatic forces and the scaling law.
In case of the viscous force, we use a form with a modified (empirical) exponent, where c x is a fit parameter, which again depends on the values of η and of V ext 0 . In Fig. 4c and d, the results of this effective scaling law are compared to the simulation data of the superadiabatic force. We find again very good agreement with the real data. Similarly to the viscous force, the structural force, the structural force is also overestimated in case of high shear rates close to the saturated state. Although we have not performed a systematic error analysis, we are confident that our data is much more consistent with the unusual value 0.8 of the exponent in Eq. 26 than with an exponent of 1, or even 0.9. We expect the power law (26) to be an empirical representation of the data rather than a fundamental relationship. We leave the construction of a corresponding approximation for the superadiabatic excess power functional, which is the generator of the nonequilibrium forces via f sup (r, t) = −δP exc t /δJ(r, t), to future work.

Conclusions
We have systematically investigated the effects that shearing has on the density and force profiles of a two-dimensional system of quasi-hard disks upon which an additional conservative confining force field is applied in the direction normal to the flow. Our results show that the shear flow induces a superadiabatic structural force that acts against the compression of the system. As a result, the density profile approaches the homogeneous bulk density for high shear rates. We have investigated the forces that act in the system. There occurs a dissipative effect from the forces acting against the flow of the particles as well as a structural force. The adiabatic construction was used to identify these superadiabatic force contributions. The concept of superadiabatic forces [31,32,44] allows systematic classification [42] of flow and structural effects, including viscous [39][40][41] and structural [41][42][43] force contributions. Future work could be directed at the behavior at low shear rates. Low shear rates are interesting as our current scaling law is not invariant if the shear direction changes, i.e., our scaling law predicts a stronger confinement for such situations. Additionally by investigation of a wider range of parameters, a more definite understanding of the scaling laws could be achieved. Furthermore, considering confining potentials that are less smooth than the one considered here should be interesting, as would be investigating the laning instability reported in Ref. [62], disorder-order transitions reported for three-dimensional systems in Ref. [63] and the general framework for viscosity of Ref. [64]. The results that we provide could form material for an in-depth comparison of PFT and MCT, possibly along the lines of Ref. [33]. Dispersions of thermoresponsive colloids [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] could be excellent model systems for corresponding experimental work that is aimed at systematically studying superadiabatic forces.
if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.