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 explicitly treating 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 behavior, arising from the interplay between Brownian motion, solvent hydrodynamics, and potential interactions.^{1–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 that 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 are currently a variety of microscopic approaches. Each of these aims to capture a particular aspect of cooperative particle motion within a limited range of shear-rates and thermodynamic parameters but fails to provide a unified global picture. While exact results can be obtained for low density systems at low shear-rate,^{5–7} systems at intermediate^{7–19} or high densities^{20–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 that correctly capture the coupling between gradients in shear-rate and density remain under development.^{25–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, therefore, 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 that “look inside” the flowing system is provided by developments in the visualization and tracking of particle motion in experiments (confocal microscopy),^{25,31–33} together with the detailed information provided by computer simulations of model systems under flow.^{34–37}

Theoretical studies of inhomogeneous fluids, both in- and out-of-equilibrium, are primarily based on the spatially varying one-body density alone and do not usually directly involve the inhomogeneous two-body correlations (although exceptions do exist^{38–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 inhomogeneous 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, therefore, 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.

## II. SUPERADIABATIC-DDFT

### A. General framework

*D*

_{0}is the diffusion coefficient,

*ϕ*is the interparticle pair potential,

*r*

_{12}= |

**r**

_{1}−

**r**

_{2}|, and

*V*

_{ext}is a time-dependent external potential. The second equation is an approximate equation of motion for the two-body density, given by

*V*

_{ad}, appearing in (2) generates a fictitious external force that stabilizes the adiabatic system. This is obtained from the Yvon–Born–Green (YBG) relation of equilibrium statistical mechanics

^{46,47}

### B. Low shear-rate solutions for the pair distribution function in bulk

*ρ*(

**r**,

*t*) →

*ρ*

_{b}. In this case, the adiabatic two-body density is both isotropic and translationally invariant. Consequently, the integral term in Eq. (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:

*ρ*

^{(2)}will not depend on the absolute position coordinate.

**v**, which, for shear applied in the direction of the

*x*-axis, with shear-gradient in the

*y*-direction, is given by

*D*

_{0}, and where the nonequilibrium pair distribution function is given by

*g*

_{eq}, can be obtained from the interaction pair potential using an appropriate equilibrium theory (we will use an integral equation closure).

*ρ*

_{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 boundary-layer in

*g*as the shear-rate is increased.

^{6}We will henceforth restrict our attention to the special case of steady-shear, $\gamma \u0307(t)\u2192\gamma \u0307$, 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, Eq. (10) reduces to

**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.

*g*

_{sup}to be a linear function of $\gamma \u0307$. (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–50}) Substituting (12) into (13) and neglecting terms quadratic and higher in $\gamma \u0307$ then yields the following linearized steady-state condition:

^{18}

**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.)

*f*(

*r*). These are given by

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

^{3}

*g*

_{sup}contributes to the off-diagonal stress tensor elements.

*η*, is obtained by dividing

*σ*

_{xy}by the shear rate. In the limit of low shear-rate, $\gamma \u0307\u21920$, Eq. (17) gives the low shear viscosity,

*η*

_{0}. Substituting Eq. (15) into Eq. (17) and using the appropriate form for the rate-of-strain tensor in shear flow, given by Eq. (A5), then yields

*f*and the integral in (18) depend on the dimensionality of the system. We then find

_{2D}=

*πρ*

_{b}/4, and the 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 Eq. (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, therefore, an approximate low-shear viscosity.

*f*(

*r*) = 1/(3

*r*

^{3}), and Eq. (19) recovers the quadratic term in the well-known low density expansion (see p. 114 in Ref. 2),

*η*

_{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 Eq. (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 hard-disks in two-dimensions, which presents a phenomenology qualitatively similar to that of hard-spheres in three-dimensions while remaining convenient for visualization of the distorted pair correlations in the *xy*-plane. The solution of the radial balance equation (16) requires 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.

*f*(

*r*) is obtained from the numerical integration of Eq. (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/2

*r*

^{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 increases, 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 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

Knowledge of the nonequilibrium pair correlation function also enables the calculation of the interaction contribution to the low-shear viscosity, *η*_{0}, via Eq. (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 Ref. 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.

*η*

_{0}could be obtained from superadiabatic-DDFT by employing the “test-particle method,” as we will discuss in Sec. IV.

## III. RUSSEL–GAST-TYPE THEORY

In Sec. II, we analyzed 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 Eq. (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 Eq. (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} which constituted one of the earliest theoretical approaches to the microstructure and rheology of colloidal suspensions at a 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

*N*− 2 particle coordinates yields the following, formally exact equation of motion for the two-body density

^{41,44}

^{41,46}

### B. Low shear-rate solutions for the pair distribution function in bulk

*ϕ*

_{mf}(

*r*), defined by $geq(r)\u2261e\u2212\beta \varphi mf(r)$.] In the low density limit,

*g*

_{eq}(

*r*) →

*e*

^{−βϕ(r)}, which yields

*g*

_{sup}is linear in the flow-rate yields the following expression to leading order:

*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

*f*

^{⋆}(

*r*) for the cases of two- and three-dimensions. We refer the reader to Appendix D for additional details regarding the derivation of Eq. (37) and the boundary conditions required to solve them.

### C. Low-shear viscosity

^{46}We emphasize that

*f*

^{⋆}(

*r*) appearing in Eq. (38) is given by solution of the radial balance equations (37) and is distinct from the function

*f*(

*r*) obtained from solution of Eq. (16). In two- and three-dimensions, Eq. (38) becomes

### D. Numerical results

*g*

_{eq}for input to the radial balance equation (37). In Fig. 4, we show scatter plots of the quantity

*f*

^{⋆}(

*r*) is obtained from numerical integration of Eq. (37) and that $r=x2+y2$. The predictions of the Russel–Gast-type closure are generally very similar to those obtained from the superadiabatic-DDFT approach of Sec. II, 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/2*r*^{2}, with accumulation and depletion at contact within the compressional and extensional quadrants, respectively. As the bulk density increases, we observe the emergence of packing oscillations. We recall that superadiatic-DDFT predicted that the radial function *f* appearing in Eq. (23) has a contact value that 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 Eq. (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 Eq. (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 increases, 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 Eq. (40), for which results are shown in Fig. 6. We find that the Russel–Gast-type approximation scheme generates a low-shear viscosity curve with values slightly larger than those 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 Eq. (9). This provides a conceptual link between density 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 Sec. 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 shear-rates considered here, the only relevant transport coefficient is the viscosity, since the lowest-order shear-induced changes to the diagonal stress tensor elements are quadratic in $\gamma \u0307$. For larger shear-rates, we would need to consider the extension of our solution ansatz (15) or (36), respectively, to higher-order in $\gamma \u0307$, 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 balance 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–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 Refs. 19 and 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, is 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 toward 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 Eq. (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 the 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 Eq. (17). We suspect that a test-particle implementation of superadiabatic-DDFT will yield a low-shear viscosity that greatly improves 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–64} It would thus be interesting to investigate whether a similar extension would be feasible 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 two-body correlations, could well provide a source of inspiration for future developments.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**S. M. Tschopp**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). **J. M. Brader**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Resources (lead); Supervision (lead); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.

### APPENDIX A: RHEOLOGICAL QUANTITIES

**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**

*κ*^{65}For the aforementioned case of shear flow, we obtain

^{6,18}This motivates our choice of ansatz for

*g*

_{sup}in Eqs. (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 rate-of-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.

**M**, the following holds:

*a*(

*r*) is a scalar function of

*r*= |

**r**|. For an arbitrary, spatially dependent vector field,

**u**, we have

*a*(

*r*)

**M**is given by

**M**is symmetric, it can be shown that

*f*(

*r*). Assuming incompressible flow, ∇ ·

**v**= 0, using identity (B1) and rewriting the projected velocity as

**r**.

**E**, the coefficients of the quadratic form $(r\u0302\u22c5E\u22c5r\u0302)$ 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.

### APPENDIX C: BOUNDARY CONDITIONS

*r*→ ∞ and

*r*→ 0.

*r*> 1. [Note that

*ϕ*(

*r*> 1) = 0 for the hard-sphere potential.] Taking the scalar product of Eq. (C4) with the radial unit vector yields

*f*(

*r*→ ∞) = 0, fully determine the function

*f*(

*r*), given that the input equilibrium radial distribution function is known.

### APPENDIX D: ALTERNATIVE FORMS OF THE RADIAL BALANCE EQUATION AND BOUNDARY CONDITIONS

In the preceding two Appendixes, we provided all the details required for the derivation and solution of the radial balance equations of superadiabatic-DDFT. Since the derivation of both the alternative radial balance equations (37) and their boundary conditions are very similar, we give here only the main equations to highlight the differences between the two approximation schemes.

*f*

^{⋆}(

*r*). Assuming incompressible flow, ∇ ·

**v**= 0, and using the identities (B1) and (B11), we can re-express Eq. (D1) in the following form:

*r*→ ∞ and

*r*→ 0.

*r*= 1 yields

*f*

^{⋆}(

*r*→ ∞) = 0, then fully determines the solutions of the radial balance equations (37).

## REFERENCES

*Soft and Fragile Matter*

*Colloidal Suspension Rheology*

*Structure and Rheology of Complex Fluids*

*Topics in Chemical Engineering*

*An Introduction to Fluid Dynamics*

*Cambridge Mathematical Library*

*Advances in Chemical Physics*

*Colloidal Dispersions*

*Cambridge Monographs on Mechanics*