Superadiabatic dynamical density functional theory for colloidal suspensions under homogeneous steady-shear

The superadiabatic dynamical density functional theory (superadiabatic-DDFT) is a promising new method for the study of colloidal systems out-of-equilibrium. Within this approach the viscous forces arising from interparticle interactions are accounted for in a natural way by treating explicitly the dynamics of the two-body correlations. For bulk systems subject to spatially homogeneous shear we use the superadiabatic-DDFT framework to calculate the steady-state pair distribution function and the corresponding viscosity for low values of the shear-rate. We then consider a variant of the central approximation underlying this superadiabatic theory and obtain an inhomogeneous generalization of a rheological bulk theory due to Russel and Gast. This paper thus establishes for the first time a connection between DDFT approaches, formulated to treat inhomogeneous systems, and existing work addressing nonequilibrium microstructure and rheology in bulk colloidal suspensions.


I. INTRODUCTION
Colloidal suspensions exhibit a rich variety of rheological behaviour, arising from the interplay between Brownian motion, solvent hydrodynamics and potential interactions [1][2][3].For example, the phenomena of shear thinning, shear thickening and yielding are relevant for many commercial products and industrial processes.In order to control and tune the rheological properties of a suspension for any specific application it is necessary to have an understanding of how the microscopic interactions between the constituents influence the macroscopic response [4].The challenge for nonequilibrium statistical mechanics is to formulate robust and accurate first-principles theories based on tractable approximation schemes which capture the essential physics while remaining sufficiently simple for concrete calculations to be performed.
For the most commonly studied situation, namely bulk suspensions under homogeneous shear, there currently exist a variety of microscopic approaches.Each of these aims to capture a particular aspect of the cooperative particle motion within some limited range of shear-rates and thermodynamic parameters, but fail to provide a unified global picture.While exact results can be obtained for low density systems at low shear-rate [5][6][7], systems at intermediate [7][8][9][10][11][12][13][14][15][16][17][18][19] or high densities [20][21][22][23][24] invoke a diverse range of approximate closure relations to account for the correlated motion of the particles.Inhomogeneous systems, for which the density and shear-rate vary in space, are more challenging to treat theoretically and appropriate closure relations, which correctly capture the coupling between gradients in shear-rate and in density, remain under development [25][26][27][28][29][30].
The clearest path to a first-principles theory of suspension rheology is to focus on the particle correlation functions and then integrate these to obtain the macroscopic rheological properties of interest.This makes it possible to connect the macroscopic constitutive relations to the nonequilibrium microstructure of the system and thus gain microscopic insight into the physical mechanisms at work [2,5].For systems with pairwise additive interparticle interactions the two-body correlations are the primary objects of interest.In the absence of hydrodynamic interactions, these enable full calculation of the stress tensor which is the key quantity of interest in rheology.Additional motivation to develop theories which 'look inside' the flowing system is provided by developments in the visualization and tracking of particle motion in experiments (confocal microscopy) [25,[31][32][33], together with the detailed information provided by computer simulations of model systems under flow [34][35][36][37].
Theoretical studies of inhomogeneous fluids, both inand out-of-equilibrium, are primarily based on the spatially varying one-body density alone and do not usually involve directly the inhomogeneous two-body correlations (although exceptions do exist [38][39][40][41][42]).This is a major difference between the standard dynamical density functional theory (standard DDFT) of inhomogeneous fluids [43,44] and the aforementioned approaches to bulk colloidal rheology.When applied to bulk systems subject to homogeneous shear flow the standard DDFT does not present, as is, a useful framework, since the one-body density is not affected by the shear and remains constant in time.There have nevertheless been attempts to supplement the one-body equation of standard DDFT with additional empirical correction terms to avert this issue [26,27,45].
However, a new DDFT framework has recently been developed, the so-called superadiabatic-DDFT [41], which treats explicitly the dynamics of the two-body correlations.This scheme presents a significant improvement in describing the dynamics of inhomogenous fluids in the presence of time-dependent external potentials by accounting, via the two-body correlations, for the structural rearrangement of the particles as the system flows [41,42].The superadiabatic-DDFT (composed of a pair of coupled equations for the one-and two-body density) per construction avoids the shortcomings of standard DDFT, since its equation for the two-body density is affected by shear and can thus be directly used to study bulk homogeneous systems.By predicting the shear-induced distortion of the pair correlations, superadiabatic-DDFT allows to obtain the viscosity as an output of the theory.This is not a quantity reachable in any standard DDFT treatment of the problem, but is a very relevant bridge to the world of rheology.
In this paper we apply superadiabatic-DDFT to bulk systems under steady-shear flow and investigate its predictions for the shear-distorted pair distribution function and the low-shear viscosity.In addition, we show that a variation of the superadiabatic-DDFT reproduces, in the bulk limit, an early theory due to Russel and Gast [10] and thus provides a generalization of their approach to the case of inhomogeneous fluids in external fields.This establishes a clear connection between existing microscopic theories of bulk colloidal rheology and DDFT approaches to the dynamics of inhomogeneous fluids.The numerical predictions of this Russel-Gast-type approach are compared with those of supradiabatic-DDFT.

A. General framework
The superadiabatic-DDFT, presented in detail in reference [41], consists of a pair of differential equations for the coupled time-evolution of the one-and two-body densities.It is applicable to systems with pairwise interparticle interactions.The first equation of superadiabatic-DDFT is the exact expression for the one-body density where The second equation is an approximate equation of motion for the two-body density, given by where the superadiabatic contribution to the two-body density is defined according to 1. Sketch of the geometry.We choose particle 1 as the origin of our cartesian coordinate system (r1 is thus implicitly fixed equal to zero).
The adiabatic two-body density, ρ ad , is obtained by evaluating the equilibrium two-body density functional at the instantaneous one-body density The adiabatic potential, V ad , appearing in (2) generates a fictitious external force which stabilizes the adiabatic system.This is obtained from the Yvon-Born-Green (YBG) relation of equilibrium statistical mechanics [46,47], applied to the nonequilibrium system.We note that the approximate equation for the two-body density (2) becomes exact in the low density limit.

B. Low shear-rate solutions for the pair distribution function in bulk
We now specialize to a bulk system, for which the external potential is set equal to zero and the one-body density becomes constant, ρ(r, t) → ρ b .In this case the adiabatic two-body density is both isotropic and translationally invariant.Consequently, the integral term in equation ( 5) vanishes and the adiabatic potential becomes constant in space.This constant can then be set equal to zero without loss of generality.To exploit the translational invariance of the system we now define the following relative and absolute coordinates that then allow us to make the replacements where we henceforth drop the subscripts on the nabla.Due to translational invariance ρ (2) will not depend on the absolute position coordinate.Since we wish to investigate systems under homogeneous flow we introduce the affine velocity field, v, which, for shear applied in the direction of the x-axis, with sheargradient in the y-direction, is given by where γ is the shear-rate (as illustrated in Fig. 1).From equation (2) we thus obtain the following equation of motion for the pair distribution function in which we note the emergence of the pair diffusion constant, 2D 0 , and where the nonequilibrium pair distribution function is given by Its superadiabatic component, which encodes the flowinduced distortion of the microstructure, is defined as g sup (r, t) = g(r, t) − g eq (r), (12) where the equilibrium radial distribution function, g eq , can be obtained from the interaction pair potential using an appropriate equilibrium theory (we will use an integral equation closure).Finding a solution to equation (10) for arbitrary values of γ and ρ b is a difficult task.Even in the low density limit considerable effort must be expended to obtain numerical solutions, due to the emergence of a boundarylayer in g as the shear-rate is increased [6].We will henceforth restrict our attention to the special case of steady-shear, γ(t) → γ, for which the time-independent, steady-state pair distribution function, g(r, t) → g(r), can be obtained (almost) analytically for arbitrary values of ρ b .In the steady-state, equation (10) In equilibrium, v(r) = 0 and the pair-current, ∇g sup (r) + g sup (r)∇βϕ(r), vanishes.This condition leads trivially to the solution g sup (r) = 0, which implies that g(r) = g eq (r), as expected.
To obtain a low shear-rate solution of the steady-state equation (13) we assume g sup to be a linear function of γ. (This will be sufficient for our present purposes.However, we note that care should be exercised when assuming linearity, since boundary-layer formation can significantly complicate the picture [19,[48][49][50].)Substitution of ( 12) into (13) and neglecting terms quadratic and higher in γ then yields the following linearized steadystate condition 2D 0 ∇ • ∇g sup (r) + g sup (r)∇βϕ(r) In order to capture correctly the anisotropy induced by the flow we make the following ansatz [18] where we have introduced the rate-of-strain tensor, E, and the isotropic radial function f (r).The rate-of-strain tensor contains all relevant information about the affine flow field, while the function f (r) depends only on the interaction potential, the bulk density and the system dimensionality.(Additional information regarding the rate-of-strain tensor is provided in appendix A.) Working through the substitution of expression (15) into equation ( 14), as described in appendix B, finally yields the radial balance equations, required to determine f (r).These are given by dg eq (r) dr in two-and three-dimensions, respectively.The boundary conditions required to solve equations ( 16) depend on the pair potential under consideration and are discussed in detail in appendix C.Although we have chosen to focus on shear, we note that the expressions presented in this subsection remain valid for any incompressible, translationally invariant flow field and so could be applied directly to, e.g.extensional flows.

C. Low-shear viscosity
For a bulk system in a steady-state, the nonequilibrium pair distribution function is related to the stress tensor according to the following exact expression [3] where I is the unit tensor.Due to the isotropy of the equilibrium pair distribution function only g sup contributes to the off-diagonal stress tensor elements.
The interaction part of the viscosity, η, is obtained by dividing σ xy by the shear rate.In the limit of low shearrate, γ → 0, equation (17) gives the low shear viscosity, η 0 .Substituting equation ( 15) into equation ( 17) and using the appropriate form for the rate-of-strain tensor in shear flow, given by equation (A5), then yields Both the function f and the integral in (18) depend on the dimensionality of the system.We then find in two-and three-dimensions, respectively.We have introduced the two-dimensional area fraction, Φ 2D = πρ b /4, and three-dimensional volume fraction, Φ 3D = πρ b /6, where all lengths have been non-dimensionalized using the characteristic diameter of the particles.If the function f were known exactly, then equation ( 19) would yield the exact low-shear viscosity.Within the present approach the adiabatic approximation employed to close the two-body equation of motion (2) yields approximate radial balance equations ( 16) for f , and thus an approximate low-shear viscosity.
For the well-studied case of low density hard-spheres in three-dimensions, equation ( 16) has the analytic solution, f (r) = 1/(3r 3 ), and equation ( 19) recovers the quadratic term in the well-known low density expansion (see page 114 in reference [2]), which applies in the absence of hydrodynamic interactions between the particles.The first term in (20) is the solvent viscosity, η s .The second term arises from the drag of the solvent on the surface of each individual sphere [51].The third term, which is exact for a system without hydrodynamic interactions, represents the influence of direct potential interactions between the particles and comes from evaluating the integral in equation (19) and then employing the Stokes-Einstein relation k B T /D 0 = 3πη s .

D. Numerical results
To investigate the predictions of the superadiabatic-DDFT approach we will focus on the special case of harddisks in two-dimensions, which presents a phenomenology qualitatively similar to that of hard-spheres in threedimensions, while remaining convenient for visualization of the distorted pair correlations in the xy-plane.Solution of the radial balance equation ( 16) requires FIG. 2. Superadiabatic-DDFT g sup .For hard-disks in two-dimensions we show the superadiabatic contribution to the pair correlation function, gsup, in units of γd 2 /2D0.Since we consider only low shear-rates gsup exhibits quadipolar symmetry, see equation (15).Increasing the bulk density leads to packing oscillations.
as input the hard-disk equilibrium radial distribution function, g eq .We obtain this quantity by employing the famous Percus-Yevick closure of the homogeneous Ornstein-Zernike equation [52,53], known to be accurate for the hard-disk system at low and intermediate bulk densities.
In Fig. 2, we show scatter plots of the rescaled superadiabatic contribution to the distorted pair correlation function, namely where r = x 2 + y 2 and f (r) is obtained from the numerical integration of equation (16).Results are shown for four different values of the bulk density, ρ b , one per panel, and are valid at low shear-rates.At the lowest considered bulk density, ρ b = 0.1, the radial function f (r) obtained numerically is very close to its low-density limit value, f (r) = 1/2r 2 , and no packing oscillations are visible in the resulting scaled g sup .In this panel we observe an accumulation of particles at contact, r = d, in the 'compressional' quadrants, defined as the region where sign(y) = −sign(x), and an opposite depletion effect in the 'extensional' quadrants, where sign(y) = sign(x).As the value of the bulk density is increased packing oscillations develop at larger values of r.These oscillations in g sup take both positive and negative values and are a nontrivial prediction of the radial balance equation (16).
The packing structure of the full pair correlation function, g(r) = g eq (r) + g sup (r), will thus differ from that  23), with a fit to simulation data, taken from reference [54].The low-shear viscosity from simulation diverges as the system approaches random close packing, whereas the superadiabatic-DDFT retains the low density limiting form for all values of the area fraction.
observed in bulk as a consequence of the applied shear.However, our numerical solutions reveal that the contact value of the radial function remains unaffected by changes in ρ b and is given by in each of the panels shown in Fig. 2. Knowledge of the nonequilibrium pair correlation function also enables calculation of the interaction contribution to the low-shear viscosity, η 0 , via equation (18).Results obtained by solving the radial balance equation (16) and evaluating the integral (18) for each chosen value of the bulk density are shown in Fig. 3 and compared with a fit to Brownian dynamics simulation data taken from reference [54].While both curves overlap for low values of ρ b , the theoretical superadiabatic-DDFT curve strongly underestimates η 0 relative to the simulation at higher densities.We note that the simulation curve diverges as the system approaches random close packing, located at Φ 2D ≈ 0.82 for monodisperse hard-disks.
Due to the density independence of the contact value (22), the above described route to obtaining the interaction contribution to the low-shear viscosity recovers its low-density form, namely Accounting for the influence of shear-flow on the bulk three-body density would give corrections to higher-order in area fraction in equation (23).However this mechanism is neglected per construction within the current approximation.The fact that we recover only the leadingorder contribution (23) to the low-shear viscosity is thus not surprising and is consistent with application of the adiabatic approximation at the two-body level.It seems to us very likely that improved predictions for η 0 could be obtained from superadiabatic-DDFT by employing the 'test-particle method', as we will discuss in section IV.

III. RUSSEL-GAST-TYPE THEORY
In the preceding section we analysed the properties of superadiabatic-DDFT for bulk systems subject to homogeneous shear.All predictions of the theory for the pair distribution function and viscosity arise from the adiabatic closure employed to arrive at equation (2).In this section we show that there exists an alternative way to implement the adiabatic closure and that this results in a new approximate equation of motion for the inhomogeneous two-body density, different from equation (2).In the bulk limit the new approximation reduces to an equation of motion for the pair distribution function first introduced in 1986 by Russel and Gast [10] and which constituted one of the earliest theoretical approaches to the microstructure and rheology of colloidal suspensions at finite volume fraction.It thus seems appropriate to call our adiabatic closure of the equation of motion for the inhomogeneous two-body density a 'Russel-Gast-type' approximation.

A. Alternative adiabatic closure
Integrating the many-particle Smoluchowski equation over N −2 particle coordinates yields the following, formally exact equation of motion for the two-body density [41,44] This equation can also be rewritten in the following alternative form which remains fully equivalent to equation (24), since there is no approximation involved at this stage.
Using the standard form (24) of the exact equation of motion and applying the adiabatic approximation solely on the three-body density, i.e. ρ (3) where the adiabatic three-body density is defined as in analogy to equation (4).Substitution of the secondorder Yvon-Born-Green (YBG2) equation [41,46], into ( 26) yields the second equation of superadiabatic-DDFT, namely equation ( 2).This is the procedure followed in reference [41].
An alternative method to implement the adiabatic closure, is to start with the (reformulated) exact equation of motion (25) and to apply the approximation to the integral term, ad (r 1 , r 2 , t) We note that equation ( 29) is no longer equivalent to the previous equation (26).The essential difference between these two options is that in equation ( 26) we approximate the joint probability density, whereas in equation (29) we approximate the conditional probability density.Substitution of the rewritten YBG2 equation, ad (r 1 , r 2 , t) into expression (29), then yields the following alternative equation of motion for the two-body density ad (r 1 , r 2 , t) .
We refer to this approximation strategy as 'Russel-Gastlike', since equation (31) reduces to the Russel-Gast equation (see reference [10]) for the pair distribution function in the bulk limit, as we will demonstrate in the following subsection.

B. Low shear-rate solutions for the pair distribution function in bulk
Following the same scheme as in subsection II B, the bulk limit of equation ( 31) is given by ∂g(r, t) ∂t = −2D 0 ∇• −∇g sup (r, t)+ g sup (r, t) g eq (r) ∇g eq (r) This provides an alternative to equation (10).(As a sidenote we mention that the second term on the right handside of equation ( 32) can be rewritten using the following rearrangement g sup (r, t) g eq (r) ∇g eq (r) = −g sup (r, t)∇ − ln g eq (r) = −g sup (r, t)∇ βϕ mf (r) , where we have introduced the 'potential of mean force', ϕ mf (r), defined by g eq (r) ≡ e −βϕ mf (r) .)In the low density limit, g eq (r) → e −βϕ(r) , which yields g sup (r, t) g eq (r) ∇g eq (r) → −g sup (r, t)∇βϕ(r), such that equations ( 32) and (10) then become equivalent (which is not the case in general).
In the steady-state at finite densities, equation (32) reduces to the condition 2D 0 ∇ • ∇g sup (r, t) − g sup (r, t) g eq (r) ∇g eq (r) Using the definition (12) and again assuming that g sup is linear in the flow-rate yields the following expression to leading order 2D 0 ∇ • ∇g sup (r, t) − g sup (r, t) g eq (r) ∇g eq (r) − ∇ • g eq (r) v(r) = 0. (35) Since equations ( 14) and ( 35) are non-equivalent the solution of equation ( 35) now requires a different ansatz than used previously.The appropriate choice in the present case is given by where the ⋆ notation makes it explicit that the (still unknown) radial function, f ⋆ (r), is different from the function, f (r), previously encountered.Substitution of ansatz ( 36) into the linearized steady-state equation ( 35) then yields the alternative radial balance equations dg eq (r) df ⋆ (r) dr g eq (r) + 4 r 3 f ⋆ (r)g eq (r), dg eq (r) to determine f ⋆ (r), for the cases of two-and threedimensions.We refer the reader to appendix D for additional details regarding the derivation of equations ( 37) and the boundary conditions required to solve them.

C. Low-shear viscosity
Within the alternative adiabatic closure the low-shear viscosity is given by d dr e −βϕ(r) y eq (r)f ⋆ (r), (38) where we have introduced the so-called 'cavity distribution function', y eq (r) = g eq (r) e βϕ(r) , familiar from liquid-state theory [46].We emphasize that f ⋆ (r) appearing in equation ( 38) is given by solution of the radial balance equations (37) and is distinct from the function f (r) obtained from solution of equations ( 16).
In two-and three-dimensions equation (38) becomes respectively.For the special case of low density hardspheres in three-dimensions, the analytic solution of For hard-disks we show the superadiabatic contribution to the pair correlation function, gsup, in units of γd 2 /2D0, see equation (36).Since the results are qualitatively similar to those in Fig. 2 and thus difficult to compare, we refer the reader to Fig. 5 for a more detailed analysis.
equation ( 37) is the same as that given by equation ( 16), namely such that the Russel-Gast-type closure reproduces the exact low density expansion of the low-shear viscosity of hard-spheres, given by equation (20).

D. Numerical results
We now investigate some numerical predictions of the alternative closure, where we again consider the harddisk system in two-dimensions and use the Percus-Yevick closure to obtain g eq for input to the radial balance equation (37).In Fig. 4 we show scatter plots of the quantity for four different values of the input bulk density.Note that the radial function f ⋆ (r) is obtained from numerical integration of equation (37) and that r = x 2 + y 2 .The predictions of the Russel-Gast-type closure are generally very similar to those obtained from the superadiabatic-DDFT approach of the previous section, although deviations become apparent on closer inspection.The results for g sup at the lowest bulk density considered, ρ b = 0.1, are very close to the known low-density limit, for which f ⋆ (r) = 1/2r 2 , with accumulation and depletion at contact within the compressional and extensional quadrants, respectively.As the bulk density is increased we observe the emergence of packing oscillations.We recall that superadiatic-DDFT predicted that the radial function f appearing in equation ( 23) has a contact value which remains independent of the bulk density.The analogous quantity within the present approximation is the product f ⋆ (r) g eq (r) appearing on the right hand-side of equation (42).We find that the contact value of this product does exhibit a nontrivial dependence on ρ b and, as we will see, this causes the low-shear viscosity generated by equation (40) to deviate from that predicted by superadiabatic-DDFT.
In Fig. 5 we show the pair distribution function on both the extensional axis, along which the particles get pulled apart by the shear, and on the compressional axis, along which the particles get pressed together.In order to show more clearly the difference between the superadiabatic-DDFT and Russel-Gast-type approximations we show in the inset the difference between the nonequilibrium and equilibrium pair distribution functions.On the extensional axis both approximation schemes predict a reduction in the height of all maxima, including the first peak at particle contact.The reduction in amplitude of the peaks is slightly more pronounced within the Russel-Gast-type approximation.In both cases the radial position of the second and higher-order peaks shifts to larger values of r, consistent with the fact that particles are being pulled away from each other by the shear flow.In contrast, on the compressional axis the height of all maxima increase and the radial locations of the second and higher-order peaks are shifted to smaller values of r, since the particles are being pushed closer together by the shear flow.The two theories again make similar predictions, although the Russel-Gast-type approximation yields a slightly larger increase in peak height than the superadiabatic-DDFT.The general similarity of the predictions from the two approximation schemes for the nonequilibrium microstructure is reassuring, as it appears that the results are robust with respect to the details of how the adiabatic approximation is implemented.
From the nonequilibrium pair correlation function we can calculate the zero-shear viscosity using equation (40), for which results are shown in Fig. 6.We find that the Russel-Gast-type approximation scheme generates a lowshear viscosity curve with values slightly larger than that predicted by the superadiabatic-DDFT.We can attribute this effect to the difference in the contact value of the pair distribution function between the two approximations, see Fig. 5. Within the Russel-Gast-type approach the difference between the contact value on the extensional axis and the contact value on the compressional axis is larger than for superadiabatic-DDFT.This larger extensional/compressional asymmetry has the consequence that for any given area fraction η 0 from the Russel-Gast-type approach is larger than for the superadiabatic-DDFT.Although the Russel-Gast-type closure produces values of the low-shear viscosity closer to the simulation results, it has the undesirable feature of changing curvature as the area fraction increases beyond about 0.45.This questions the physicality of the alternative approximation and its rescaling potential for higher densities [12].

IV. DISCUSSION
In this paper we have investigated the predictions of superadiabatic-DDFT for the nonequilibrium pair correlation function and low-shear viscosity of bulk systems subject to the homogeneous shear flow defined by equation (9).This provides a conceptual link between density FIG. 6. Low-shear viscosity.Analogous plot to Fig. 3, but now including the prediction of the Russel-Gast-type approximation.This additional result predicts slightly larger values of η0 than the superadiabatic-DDFT, but exhibits an unphysical curvature as the area fraction is increased.functional-based approaches, which are focused on the inhomogeneous one-body density [43,44,55], and microscopic theories of homogeneous bulk rheology, of the type pioneered by Brady, Russel, Wagner and others [2,3,18].A clear connection between these two substantial bodies of literature, which have coexisted for decades with little to no interaction, is established by our Russel-Gast-type closure of the equation of motion for the two-body density.As detailed in section III, this new scheme generates a fully inhomogeneous DDFT, which reduces to the known Russel-Gast theory for bulk systems subject to homogeneous shear [10].
The application of shear distorts the bulk pair distribution function and, for systems interacting via a pairwise additive potential, enables the interaction part of the stress tensor to be calculated.For the low shearrates considered here the only relevant transport coefficient is the viscosity, since the lowest-order shearinduced changes to the diagonal stress tensor elements are quadratic in γ.For larger shear-rates we would need to consider the extension of our solution ansatz (15) or (36), respectively, to higher-order in γ, which would then generate nonvanishing normal stress differences and shear-induced modifications to the system pressure [6,29].It would be interesting to investigate the extent to which numerical solutions of the superadiabatic-DDFT can account for nonlinear rheological and microstructural changes as the shear rate is increased.For example, it is known from both simulation [37] and experiment [56] that the nonequilibrium bulk three-body distribution function deviates significantly from its equilibrium form, which could be difficult to capture using an adiabatic closure on the two-body level.
Focusing on the low shear-rate regime enabled a largely analytic investigation of the nonequilibrium pair distribution function, with only the solution of the radial bal-ance equation and evaluation of the low-shear viscosity demanding (relatively simple) numerical integration.At larger values of the shear-rate more sophisticated techniques must be employed, since boundary-layer formation prevents the application of straightforward perturbation theory [6,[48][49][50]57].An interesting extension of our study to obtain results for higher shear-rates could be to employ either multipole methods [58], specialized numerical schemes [6,59] or intermediate asymptotics [19,49].Although the theories presented here are valid for an arbitrary pair interaction potential in two-or three-dimensions, we chose to perform numerical calculations for the special case of hard-disks, as this presents the simplest continuum model for the study of shear, while still retaining much of the essential phenomenology of three-dimensional systems.In references [19,49] it is shown how the pair distribution function under shear can be calculated analytically for more general interaction potentials.
The two routes presented in this work are built on two variations of the same approximation.The numerical output of those theories is broadly similar but, as clearly pictured by the final viscosity plot, Fig. 6, are different in more demanding situations at higher bulk density.One can thus ask which scheme is more likely to succeed in further investigations.The unphysical shape of the viscosity curve obtained with the Russel-Gast-like closure suggests that this variation of the approximation may be less robust than that of superadiabatic-DDFT and could prove more difficult to improve and refine.
A natural first step towards better predictions for the viscosity at higher bulk densities would be to apply Brady's semi-empirical rescaling method [12].For the present Brownian dynamics, without hydrodynamic interactions, this consists of multiplying the low density limiting solution for the interaction contribution to η 0 (given by equation (23) for hard-disks) by the contact value of the equilibrium radial distribution function.The resulting low-shear viscosity has been shown to agree well with simulation data for hard-spheres [16] and we find that this is also the case for hard-disks.Since the superadiabatic-DDFT predicts that the low-density limiting result (23) holds for all values of ρ b , a rescaled version of superadiatic-DDFT would exactly reproduce Brady's theory for η 0 .In contrast, a rescaled Russel-Gast-type theory would continue to exhibit an unphysical curvature, as the contact value of g eq is a monotonically increasing function of bulk density.Investigating how such a rescaling factor could emerge from a systematic, first-principles extension of the superadiabatic-DDFT closure scheme will be an interesting topic for future study.
While the above considerations may provide a path to structural improvements of the superadiabatic-DDFT equations (i.e.following from a new two-body closure) there remains a great deal to be learned from the current version of the theory.In this paper we have considered only the direct application of superadiabatic-DDFT to the bulk system, meaning that we work from the outset with a constant one-body density and focus on the properties and predictions of the resulting two-body equation of motion.Although this sheds light on the internal structure and physical content of the theory, we should keep in mind that the superadiabatic-DDFT is still a density functional theory and, as such, aims primarily to predict the dynamics of the one-body density.If we are interested solely in bulk systems under shear, as in the present work, then the full power of superadiabatic-DDFT for the one-body dynamics can be harnessed by employing a more sophisticated implementation scheme: the test-particle method.Fixing a test-particle at the coordinate origin (i.e.setting V ext (r) = ϕ(r)) induces a spatially varying steady-state one-body density, ρ(r), as particles around the test-particle accumulate in the compressional quadrants and are depleted from the extensional quadrants [54].This inhomogeneous one-body density can be related to the bulk pair distribution function according to g(r) = ρ(r)/ρ b , from which the viscosity can be calculated using equation (17).We suspect that a test-particle implementation of superadiabatic-DDFT will yield a low-shear viscosity which much improves on the simple quadratic expression (23), obtained by direct application of the two-body equation of motion to bulk.The proposed test-particle calculation would require explicit treatment of the inhomogeneous two-body density in the presence of a test-particle and would thus capture, to some extent, the shear-induced distortion of the bulk three-body correlations.Investigations in this direction are underway.
Finally, we mention the issue of hydrodynamic interactions.Although the majority of standard DDFT studies do not consider solvent hydrodynamics, this aspect has been incorporated into the formalism [60][61][62][63][64].It would thus be interesting to investigate whether a similar extension would be feasable for superadiabatic-DDFT.Many of the existing theories of homogeneous bulk rheology mentioned in the introduction, including the bulk Russel-Gast theory, were formulated to include hydrodynamic interactions to some level of approximation.These works, all of which are focused on the two-body correlations, could well provide a source of inspiration for future developments.
Appendix A: Rheological quantities In a system undergoing translationally invariant, homogeneous flow, the velocity field can be conveniently expressed in the following form where κ is the spatially constant velocity gradient tensor.As an example, choosing flow in the x-direction, with shear-gradient in the y-direction enables κ to be represented in matrix form as for a three-dimensional system with shear-rate γ.The right hand-side of (A1) can be decomposed into a sum of two terms, where are the (symmetric) rate-of-strain tensor and the (antisymmetric) rate-of-rotation tensor, respectively.The rate-of-strain tensor describes 'pure straining motion' leading to relative motion between any two particles, whereas the rate-of-rotation tensor yields pure rotational motion, which does not affect their relative separation [65].For the aforementioned case of shear flow we obtain The anisotropy of the nonequilibrium pair distribution function in a system subject to a slow translationally invariant flow is determined solely by its straining motion and does not involve its rotational component [6,18].This motivates our choice of ansatz for g sup in equations ( 15) and (36).
Appendix B: The radial balance equation Substitution of the ansatz (15) into the linearized steady-state equation ( 14) yields the radial balance equations (16).Since this procedure is not straightforward we outline here the main steps of the calculation.Although we are primarily interested in the case of homogeneous shear, we formulate the problem using the rateof-flow tensor E, which can also describe other flow types.This not only extends the generality of our results, but also presents the technical advantage that certain vector/tensor identities can be employed to make the calculation as clean as possible.Unless otherwise stated, all expressions given in this appendix are valid in arbitrary dimensionality.
We start by collecting a few useful identities.For a general, spatially varying, second-rank tensor M the following holds where r = r/r is a unit vector, rr is a dyadic product and where a(r) is a scalar function of r = |r|.For an arbitrary, spatially dependent vector field, u, we have where the double dot notation in the second term indicates a full contraction, i.e. a scalar product followed by a trace operation.The divergence of the product a(r)M is given by The advantage of this representation is the appearance of the dyadic tensors, rr and (I − rr), which project either along or perpendicular to the relative position vector r.
On the right hand-side of (B12) we have to calculate the divergence of the scalar product between the vector r • E and a second-rank tensor.We can thus exploit relation (B2) to obtain (B13) We will consider separately the two terms appearing on the right hand-side of (B13), which we henceforth refer to as (i) and (ii).
To simplify (i) in two-dimensions we use (B3), (B4) and (B5).This yields Since these expressions remain valid for all choices of the traceless tensor E, the coefficients of the quadratic form (r•E•r) must be equal.We thus obtain the radial balance equations ( 16) stated in the main text.We note that, in addition to the explicit appearance of the pair interaction potential in (B17), respectively (B18), both the pair potential and the bulk density enter these equations implicitly via the equilibrium radial distribution function.
which is different from equation (C6) The boundary condition at particle contact, equation (D5), together with f ⋆ (r → ∞) = 0, then fully determines the solutions of the radial balance equations (37).

FIG. 3 .
FIG.3.Low-shear viscosity.A comparison of η0 from superadiabatic-DDFT, see equation(23), with a fit to simulation data, taken from reference[54].The low-shear viscosity from simulation diverges as the system approaches random close packing, whereas the superadiabatic-DDFT retains the low density limiting form for all values of the area fraction.

FIG. 4 .
FIG. 4. Russel-Gast-type closure for g sup .Analogous plot to Fig. 2, but this time employing the alternative closure.For hard-disks we show the superadiabatic contribution to the pair correlation function, gsup, in units of γd 2 /2D0, see equation(36).Since the results are qualitatively similar to those in Fig.2and thus difficult to compare, we refer the reader to Fig.5for a more detailed analysis.

FIG. 5 .
FIG. 5. Pair distribution function.We show the full nonequilibrium pair distribution function of hard-disks along the extensional axis, y = x, and the compressional axis, y = −x.Results are given from both the superadiabatic-DDFT and the Russel-Gast-type approximations, for ρ b = 0.7.The additional black curve indicates the equilibrium radial distribution function for comparison.The dimensionless shear-rate is chosen to be γd 2 /2D0 = 1.The inset shows the difference between the nonequilibrium and equilibrium pair distribution functions to highlight the difference between the two approximations.Due to the symmetry of the low-shear pair distribution function for a given approximation the curves in the compressional and extensional directions are identical (up to a sign).