Heterogeneous Solvent Dissipation Coupled with Particle Rearrangement  in Shear Thinning Non-Brownian Suspensions

Tomoharu Terayama
[Present Address:] Department of Physics, The University of Tokyo, Hongo 7-3-1, Tokyo 113-0033, Japan.
Akira Furukawa
Institute of Industrial Science, University of Tokyo, Meguro-ku, Tokyo 153-8505, Japan.

Dense non-Brownian suspensions exhibit significant shear thinning, although a comprehensive understanding of the full scope of this phenomenon remains elusive.
This study numerically reveals intimate heterogenous coupled dynamics between many-body particle motions and solvent hydrodynamics in shear-thinning non-Brownian suspensions. We demonstrate the spatially correlated viscous dissipation and particle motions; they share the same characteristic length, which decreases with increasing shear rate. We further show that, at lower shear rates, significant particle density changes are induced against the incompressibility of the solvent, suggesting the cooperative creation and annihilation of gaps and flow channels. We discuss that hydrodynamic interactions may substantially restrict particle rearrangements even in highly dense suspensions, influencing the quantitative aspects of macroscopic rheology.

§ INTRODUCTION

Understanding and controlling the rheological properties of non-Brownian suspensions, especially the shear-thickening and shear-thinning behavior, are crucial in many industrial applications and are fundamental problems in engineering and basic sciences [1, 2, 3, 4, 5].
The shear-thickening mechanism is currently attracting considerable attention, and recent studies [6, 7, 8] proposed a convincing mechanism underlying  the essential role of interparticle contact friction in the steep increase in viscosity.
In contrast, the shear-thinning mechanism is still under debate [4, 5, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18].
There are two shear-thinning regimes separated by a shear-thickening regime, and recent rheological experiments combined with mechanical force measurements indicate the following [9, 19]:
In the high-shear-rate post-thickening regime, contact friction still dominates the dynamics but is weaker with increasing shear rate, causing shear thinning [20]. In contrast, in the low-shear-rate pre-thickening regime, constituent particles are not as strongly compressed against each other such that instead of the frictional contact between particles, the lubrication hydrodynamics should be relevant to particle motion.
Under such conditions, as the shear rate and the resultant particle pressure increase, the effective particle gaps (volumes) are increased (decreased), reducing the overall flow resistance [9, 10, 11, 12].  Although these physical pictures are intuitively appealing, further investigations are required to clarify the physical mechanism of shear thinning at the same level as that of shear thickening.

When the frictional contact forces are irrelevant to particle motion for lower shear rates, the viscous dissipation due to the solvent dynamics dominates the net dissipation. Furthermore, the incompressibility of the solvent and the resultant hydrodynamic interactions are expected to strongly constrain particle dynamics [21, 22]. A very recent experimental study certainly indicates that a fluid component (solvent) determines the quantitative aspects of the flow and rheological properties of dense suspensions [23, 24]. However, understanding how the solvent hydrodynamics are involved in suspension rheology remains elusive. In this study, using a two-dimensional model non-Brownian suspension, we explore the dynamic coupling between the particle motion and solvent hydrodynamics and its link with the shear-thinning behavior through a direct numerical simulation.

§ MODELS AND SIMULATION METHODS

In this study, our model of a two-dimensional non-Brownian suspension consists of a mixture of large (L) and small (S) disk-shaped particles in a $50:50$ ratio. The size ratio between large and small particles is $1.4$.
This setup can prevent crystallization, similar to the model system used in previous studies of supercooled liquids [25, 26].
The $i$-th and $j$-th particles interact via the following very steep soft-core or inverse power law potentials:

\begin{eqnarray}
U(r_{ij})=\kappa \biggl(\dfrac{s_{ij}}{r_{ij}}\biggr)^{36}\Theta\biggl(\dfrac{r_{ij}}{s_{ij}}-1.4\biggr),  \label{potential}
\end{eqnarray}

where $s_{ij}=(s_{i}+s_{j})/2$, with $s_i$ being the $i$-th particle size, $r_{ij}$ is the distance between the two particles, and $\Theta(r_{ij}/s_{ij}-1.4)$ is the Heaviside step function that provides a cut-off length.
Here, the soft-core sizes of these particles are denoted as $s_{\rm L}$ and $s_{\rm S}$ for large and small particles, respectively.

Our simulations employ the smoothed-profile method (SPM) [27, 28, 29] to deal with many-body hydrodynamic interactions (HIs) among the constituent particles.
For this purpose, the $i$-th particle is represented through the following field variables

\begin{eqnarray}
\psi_i ({{\mbox{\boldmath$r$}}})=
\begin{cases}
{1 ~~~~ (|{{\mbox{\boldmath$r$}}}-{{\mbox{\boldmath$R$}}}_{i}|-R < -{\xi}/{2})}\\\\
{\dfrac{1}{2}\biggl\{1-\sin\biggl[\dfrac{\pi(|{{\mbox{\boldmath$r$}}}-{{\mbox{\boldmath$R$}}}_{i}|-R)}{\xi}\biggr]\biggr\}} \\
(-{\xi}/{2} \le |{{\mbox{\boldmath$r$}}}-{{\mbox{\boldmath$R$}}}_{i}|-R \le {\xi}/{2})} \\\\
{0 ~~~~ (|{{\mbox{\boldmath$r$}}}-{{\mbox{\boldmath$R$}}}_{i}|-R > {\xi}/{2})}
\end{cases},
\end{eqnarray}

where ${{\mbox{\boldmath$R$}}}_{i}$ is the position of the $i$-th particle, $R=s_i/2$,  and $\xi$ is the interface thickness controlling the degree of smoothness.

The working equations for the velocity field  ${{\mbox{\boldmath$v$}}}({{\mbox{\boldmath$r$}}},t)$ are given as

\begin{eqnarray}
\rho\biggl(\dfrac{\partial}{\partial t}+ {{\mbox{\boldmath$v$}}}\cdot {{\mbox{\boldmath$\nabla$}}}\biggr){{\mbox{\boldmath$v$}}} &=&  {{\mbox{\boldmath$\nabla$}}}\cdot {\stackrel{\leftrightarrow}{\mbox{\boldmath$\Sigma$}}_{vis}}-{{\mbox{\boldmath$\nabla$}}}p+ {{\mbox{\boldmath$f$}}}_{H},
\label{Navier_Stokes} \\
{\stackrel{\leftrightarrow}{\mbox{\boldmath$\Sigma$}}_{vis}} &=&  \eta_s \bigl[{{\mbox{\boldmath$\nabla$}}} {{\mbox{\boldmath$v$}}}+ ({{\mbox{\boldmath$\nabla$}}} {{\mbox{\boldmath$v$}}})^T\bigr], \label{stress_tensor} \\
{{\mbox{\boldmath$\nabla$}}}\cdot {{\mbox{\boldmath$v$}}}&=&0.  \label{incompressibility}
\end{eqnarray}

Equation (<ref>) is the usual Navier-Stokes equation [30].
Here, $\rho$ is the mass density and ${\stackrel{\leftrightarrow}{\mbox{\boldmath$\Sigma$}}_{vis}}$ given in Eq. (<ref>) is the viscous stress tensor, with $\eta_s$ being the solvent viscosity.
Additionally, the hydrostatic pressure $p$ is determined by the incompressibility condition, Eq.  (<ref>), and ${{\mbox{\boldmath$f$}}}_{H}$ is the body force required to satisfy the rigid body condition [29].

As described in the main text, periodic boundary conditions are imposed in the $x$-direction with the linear dimension $L$, and planar top and bottom walls are placed at $y=H/2$ and $-H/2$, respectively. Shear flow is imposed by moving the top and bottom walls in the $x$-direction at constant velocities $V/2$ and $-V/2$, respectively, whereby the mean shear rate is given as $\dot\gamma=V/H$.
We impose no-slip boundary conditions at the top and bottom walls, namely, ${{\mbox{\boldmath$v$}}}(x,y=H/2)=(V/2){\hat {\mbox{\boldmath$x$}}}$ and ${{\mbox{\boldmath$v$}}}(x,y=-H/2)=-(V/2){\hat {\mbox{\boldmath$x$}}}$.
Furthermore, large and small particles, which are identical to those in the inner region, are embedded in the walls. We let ${{\mbox{\boldmath$R$}}}_{i}^{(w)}=(x_{i}^{(w)},\pm H/2)$ be the position of the $i$-th particle comoving with the walls.
The wall particles are randomly distributed with the distance of the nearest neighbor pair, $i$ and $j$, being larger than $s_{ij}$ and smaller than $2s_{ij}$.
This treatment prevents the penetration of the particles through the boundary walls and causes the walls to transmit the forces from the particles.

Regarding the particles that are not on the walls, the equations of motion of the $i$-th particle velocity, ${{\mbox{\boldmath$V$}}}_i$, and the angular velocity,
${{\mbox{\boldmath$\Omega$}}}_i$, are

\begin{eqnarray}
M_i \dfrac{d{{\mbox{\boldmath$V$}}}_i}{dt} &=& {{\mbox{\boldmath$F$}}}_{i, {H}} + {{\mbox{\boldmath$F$}}}_{i, {int}}
+ {{\mbox{\boldmath$F$}}}_{i, {ex}},  \label{VG}
\\
J_i \dfrac{d{{\mbox{\boldmath$\Omega$}}}_i}{dt} &=&  {{\mbox{\boldmath$N$}}}_{i, {H}} + {{\mbox{\boldmath$N$}}}_{i, {ex}},   \label{OG}
\end{eqnarray}

where $M_{i} = \rho \pi s_{i}^2/4$ and $J_i =M_i s_i^2/8$ are the mass and moment of inertia of the $i$-th particle, respectively.
In this study, the particle density is assumed to be the same as the solvent density.
In Eqs. (<ref>) and (<ref>),  ${{\mbox{\boldmath$F$}}}_{i, {int}}$
is the force acting on the $i$-th particle due to the particle-particle interactions:

\begin{eqnarray}
{{\mbox{\boldmath$F$}}}_{i, {int}}= -\sum_{j\ne i} \dfrac{\partial}{\partial {{\mbox{\boldmath$R$}}}_{i}^{}}U(|{{\mbox{\boldmath$R$}}}_{i}^{}-{{\mbox{\boldmath$R$}}}_{j}^{}|),
\end{eqnarray}

where $U$ is given in Eq. (<ref>).

In Eqs. (<ref>) and (<ref>),
${{\mbox{\boldmath$F$}}}_{i, {ext}}$ and ${{\mbox{\boldmath$N$}}}_{i, {ext}}$ are the force and torque exerted on the $i$-th particle due to the external field, respectively,  which are absent in the present study.
Finally, ${{\mbox{\boldmath$F$}}}_{i, {H}}$ and ${{\mbox{\boldmath$N$}}}_{i, {H}}$ are the force and torque exerted on the $i$-th particle due to HIs:

\begin{eqnarray}
{{\mbox{\boldmath$F$}}}_{i, {H}}&=&  -\int {\rm d}{{\mbox{\boldmath$r$}}} {{\mbox{\boldmath$f$}}}_{H},
\label{PFrb}\\
{{\mbox{\boldmath$N$}}}_{i, {H}} &=&
-\int {\rm d}{{\mbox{\boldmath$r$}}} ({{\mbox{\boldmath$r$}}}-{{\mbox{\boldmath$R$}}}_i) \times {{\mbox{\boldmath$f$}}}_{H}.
\label{PTrb}
\end{eqnarray}

Equation (<ref>) enforces the momentum exchange between the solvent and the particle.
The body force ${{\mbox{\boldmath$f$}}}_{H}$ can be determined to approximately fulfill the rigid body condition inside the particle region. Its explicit form is given in the discretized equations of motion (for details, please refer to the literature of the SPM [27, 28, 29]).

In our simulations, we make the equations dimensionless by adopting $\ell$ and $t_0$ as the units for measureing space and time, respectively. Here, $\ell$ represents the discretization mesh size used in solving Eqs. (<ref>)-(<ref>), and $t_0=\rho \ell^2/\eta_s$ represents the momentum diffusion time across the unit length.
As a result, the scaled solvent viscosity is set to $1$, and the units for velocity, stress, force, and energy are defined as $\ell/t_0$, $\rho \ell^2/t_0^2$, $\rho \ell^4/t_0^2$ and $\rho \ell^5/ t_0^2$, respectively.
In our simulations, we set $\kappa=1$ and $H=L=512$.
The sizes of the large and small particles are $s_{\rm S}=$8.0 and $s_{\rm L}=$11.2, respectively, and $\xi=1$.
The area fraction of the particles is defined as

\begin{eqnarray}
\phi=\dfrac{\pi}{8{\mathcal A}}(2N_{\rm L}s_{\rm L}^2+2N_{\rm S}s_{\rm S}^2+N_{\rm L}^{(w)}s_{\rm L}^2+N_{\rm S}^{(w)}s_{\rm S}^2),
\end{eqnarray}

where $\mathcal A=H\times L$ represents the area of the system.
Here, $N_{\rm L}$ and $N_{\rm S}$ represent the number of large and small particles not located on the walls,  respectively, with $N_{\rm L}/N_{\rm S}=1$.
On the other hand, $N_{\rm L}^{(w)}$ and $N_{\rm S}^{(w)}$ denote the number of large and small particles embedded in the walls, respectively, with $N_{\rm L}^{(w)}/N_{\rm S}^{(w)}=1$.
Additionally, we ensure that the numbers of particles on the top and bottom walls are identical.

The hydrodynamic interactions (HIs) among disks, satisfying the solvent incompressibility, are incorporated by adopting the smoothed-profile method (SPM) [27, 28, 29],
which can accurately reproduce even the near-field HIs [29]:  In the present study, the ratio of the simulation mesh size ($\ell$) to the particle size ($s$) is approximately 0.1, determining the spatial resolution of HIs. With this setting, the SPM can quantitatively replicate the near-field HIs (lubrication interactions) to a satisfactory degree up to a closer distance as $h/s\sim 0.1$ (or even less) [29], where $h$ represents the gap distance, measured as the separation distance between the interfaces of adjacent particles.
However, to fully reproduce the singular divergence of lubrication forces, which become crucial at much closer distances of $h/s\lesssim 0.01$, finer resolutions are required. This situation is similar to other hybrid simulation methods [31, 32, 33, 34, 35, 36].
Nevertheless, in systems of particles with soft or rough surfaces, achieving such conditions is difficult.
Furthermore, we do not include frictional contact forces in the present model.
They should be crucial for hard suspended particles at very close contact, while they are expected to be less involved in the actual shear-thinning behavior at lower shear rates [9, 19].

§ RESULTS

§.§ Shear thinning behaviors

(a) $\phi$ dependence of the viscosity $\eta$ for several shear rates $\dot{\gamma}$. The inset shows a closer view for $0 \le \phi \le 0.2$. For $\phi\lesssim 0.4$, $\eta$ almost collapses onto a single curve, while for $\phi\gtrsim 0.4$, $\eta$ increases more steeply for smaller $\dot\gamma$.
(b) Viscosity vs. shear rate for various area fractions. The viscosity measured at the walls, $\eta$, is plotted by open circles, while that measured through the dissipation rate, ${\bar \eta}$, is plotted by filled squares. The thinning exponent exhibits a slight dependence on $\phi$ and $\dot\gamma$, and it appears to vary around -0.6 in the thinning regime.

In this study, the macroscopic viscosity is defined as $\eta=(1/{\dot \gamma}L)\int {\rm d}x \langle \Sigma_{xy}(x,\pm L/2)\rangle$, where $\Sigma_{xy}(x,\pm L/2)$ is the $xy$ component of the stress tensor at the walls and $\langle \cdots \rangle$ hereafter denotes taking the time average in a steady state.
In Fig. <ref>(a), we show the area fraction ($\phi$) dependence of $\eta$ for several shear rates.
For $\phi\lesssim 0.4$, $\eta$ almost collapses onto a single curve. In contrast, for $\phi\gtrsim 0.4$, $\eta$ increases more steeply for smaller $\dot \gamma$.
For such higher $\phi$, the force due to the interaction potential predominantly contributes to the net stress, but it increases less linearly with $\dot\gamma$ on average. Therefore, the resultant viscosity is a decreasing function of $\dot\gamma$.
This situation is more clearly exhibited in Fig. <ref>(b), where $\eta$ is plotted against $\dot\gamma$: for the three larger $\phi$ values$(\gtrsim 0.5)$, $\eta$ exhibits marked shear-thinning behavior with a clear plateau at higher shear rates, but we do not find any plateau at lower shear rates in the present ranges of $\dot\gamma$ and $\phi$.
Note that, in our simulations, the Reynolds number defined as $\rho{\dot \gamma}s_{\rm S}^2/\eta_s$ [37], with $\rho$ and $\eta_s$ being the mass density and the solvent viscosity, respectively, spans a range from approximately $10^{-4}$ to $10^{-1}$ for the present range of $\dot\gamma$.

§.§ Heterogeneous solvent dissipation dynamics

(Color online) Typical snapshots of $i_2/\dot\gamma^2$ at $\phi=0.27$ (a)-(c), and those at $\phi=0.64$ (d)-(f) are shown.
In the insets of (a) and (d), we show the expanded images of small areas. In these insets, although the particles are represented as red and blue circles, we do not distinguish the particle and solvent regions in the main images.
At $\phi=0.64$, some regions transiently show very steep velocity gradients, whereas those gradients in the other regions are gentle, resulting in the spatially heterogeneous viscous dissipation dynamics. The characteristic size of such heterogeneous dynamics is reduced as $\dot\gamma$ is increased.

As mentioned above, in the present study, frictional contact forces are not included, and therefore, the dissipation is entirely due to the solvent viscosity.
The average dissipation rate density, ${\dot \epsilon}_{diss}$, can be written as

\begin{eqnarray}
{\dot \epsilon}_{diss} :=\dfrac{1}{2}{\bar \eta} {\dot\gamma}^2 = \dfrac{1}{\mathcal A}\int_{\mathcal A} \dfrac{1}{2}{\rm d}{ {\mbox{\boldmath$r$}}} \eta_s \langle \bigl(\nabla { {\mbox{\boldmath$v$}}} + \nabla { {\mbox{\boldmath$v$}}}^{\rm T} \bigr)^2 \rangle.
\end{eqnarray}

Then, the effective viscosity, ${\bar \eta}$, is given by

\begin{eqnarray}
{\bar \eta}= \dfrac{\eta_s}{{\dot\gamma}^2 {\mathcal A}}\int_{\mathcal A} {\rm d}{ {\mbox{\boldmath$r$}}}  {\langle i_2 \rangle}.   \label{effective_viscosity}
\end{eqnarray}

Here, we introduce the variable $i_2$ as

\begin{eqnarray}
i_2 &:=& \bigl(\nabla { {\mbox{\boldmath$v$}}} + \nabla { {\mbox{\boldmath$v$}}}^{\rm T} \bigr)^2 \nonumber \\
&=& 2[(\nabla_x v_x -\nabla_y v_y)^2 +(\nabla_x v_y + \nabla_y v_x)^2],
\end{eqnarray}

which is positive in the fluid regions and zero in the particle regions. Note that $(i_2/2)^{1/2}$ represents the rotationally invariant shear gradient and $\eta_s i_2/2$ is the local dissipation rate density. In Fig. <ref> (b), we plot the effective viscosity $\bar\eta$ defined in Eq. (<ref>) in addition to the viscosity measured at the walls, $\eta$. In principle, $\bar\eta$ should coincide with $\eta$, but $\eta$ is significantly larger than $\bar\eta$ at higher $\phi$. Such a deviation between $\eta$ and $\bar\eta$ may be attributed to the intrinsic properties of the SPM: in the SPM framework [27, 28, 29], the equations of motion are constructed to precisely maintain the momentum exchange between the degrees of freedom of the solvent and particles, but this is not the case for the energy exchange.
This apparent violation of energy conservation is introduced by the smoothness of the solid-liquid interface of the particles.
Our preliminary simulations show that when more sharp fluid-particle interfaces are imposed with the particle sizes and total number of particles remaining unchanged, the deviation is significantly reduced (not shown here). Hence, we expect that the quantitative discrepancy between $\eta$ and $\bar\eta$ does not affect at least the qualitative aspects of our arguments and results presented below.

(Color online) Probability distribution function of $i_2/(\eta_s \dot\gamma^2)$ at $\phi=0.27$ (a) and $0.64$ (b) for various $\dot\gamma$.

In Fig. <ref>, we plot typical snapshots of the spatial patterns of $i_2/\dot\gamma^2$ for two different $\phi$ with varying $\dot\gamma$. At $\phi=0.27$, where the shear-thinning behavior is hardly observed, the spatial patterns of $i_2/\dot\gamma^2$ remain apparently unchanged. In contrast, at $\phi=0.64$, where marked shear-thinning behavior is observed, we can see a significant $\dot\gamma$ dependence.
At higher $\dot\gamma(=10^{-4})$, a rather homogeneous spatial pattern of $i_2/\dot\gamma^2$ is observed. However, at lower $\dot\gamma(=10^{-6})$, $i_2/\dot\gamma^2$ is highly heterogeneous, exhibiting a visually clear long-range correlation of much steeper velocity gradients than the average gradient.
Although any siginificant anisotropic structures, such as slip or segregation layers and shear bands, are hardly observed near the boundary walls, notable anisotropies in the solvent dynamics are present in the inner region.
One such anisotropy is observed in $i_2$ around the particles, particularly evident in relatively less dense regimes, as shown in Figs. 2(a)-(c):
$i_2$ around each particle domain exhibits square-like patterns, which does not indicate structural anisotropy but results from the definition of $i_2$ [38]. Another observed anisotropy is the string-like patterns in $i_2$, extending horizontally along the $x$-axis and vertically along the $y$-axis, at higher $\dot\gamma$ in dense regimes, as seen in Figs. 2(e) and (f). There are numerous short segments of anisotropic string patterns, characterized by higher values of $i_2$ (in Figs. 2(e) and (f)), which are interpreted as localized manifestations of gliding motions along the shear and flow directions.

To supplement the observed distinction in the spatial patterns of $i_2$, we calculate the probability distribution of $i_2/(\eta_s\dot\gamma^2)$, $p[i_2/(\eta_s{\dot\gamma^2})]$, as shown in Fig. <ref>.
At $\phi=0.27$, except for at the two higher $\dot\gamma$ values, $p[i_2/(\eta_s{\dot\gamma^2})]$ collapses onto a single curve, which is responsible for the almost flat $\dot\gamma$ dependence of $\eta$.
In contrast, at $\phi=0.64$, with decreasing $\dot\gamma$, the functional form of $p[i_2/(\eta_s{\dot\gamma^2})]$ becomes broader, with its tail extending to larger values ($i_2/\dot\gamma^2\gg1$). That is, steeper velocity gradients than the applied shear rate $\dot\gamma$ are more enhanced as $\dot\gamma$ is lowered.

Structure factor $S(k)$ of $i_2$ against $ks_{\rm S}/2\pi$ at $\phi=0.27$ (a), $0.41$ (b), $0.54$ (c), and $0.64$ (d) plotted as solid curves. In our simulations, the available wave numbers are given by $k=2\pi n/L$ with $n$ being an integer.
For larger $\phi$, (c) and (d), $S(k)$ is fitted to the Ornstein-Zernike form, $S(0)/(1+\xi^2k^2)$ (broken curves), at smaller $k$.
This fitting determines the correlation length $\xi$, which decreases with increasing  $\dot{\gamma}$. The insets of (c) and (d) show $\xi$ scaled by the size of the small particle $s_{\rm S}$.

To further quantify the spatial correlation of the viscous dissipation in the solvent, we define the structure factor of the local velocity gradient as

\begin{eqnarray}
S( { {\mbox{\boldmath$k$}}} ) =  \dfrac{1}{\mathcal A}\langle |{\hat i}_{2}(\mbox{\boldmath$k$})|^2 \rangle.
\end{eqnarray}

Here, the Fourier transform of $i_2({\mbox{\boldmath$r$}})$ is defined by
${\hat i}_{2}(\mbox{\boldmath$k$})=\int d{\mbox{\boldmath$r$}} e^{-i{\mbox{\boldmath$k$}}\cdot{\mbox{\boldmath$r$}}}
Although a window function is not used here, the resultant structure factor for $k>2\times 2\pi/L$ is hardly altered by whether implementing it or not.
In the subsequent analysis, our focus is on capturing the trend of decreasing characteristic length scale with an increase in $\dot\gamma$. Despite observing some anisotropies depending on $\phi$ and $\dot\gamma$, replacing $S( { {\mbox{\boldmath$k$}}})$ with its angle average $S(k)$ may not significantly affect the results.
Figure <ref> displays $S(k)$ for various conditions.
As expected from the snapshots of $i_2$, $S(k)$ is significantly enhanced at smaller $k$ in the strongly shear-thinning regime, which confirms the spatially correlated viscous dissipation at larger $\phi$ and smaller $\dot\gamma$. Moreover, $S(k)$ at lower $k$ can be fit to the Ornstein-Zernike form

\begin{eqnarray}
S(k)=\dfrac{S(k\rightarrow 0)}{1+\xi^2 k^2},
\end{eqnarray}

from which we can determine the correlation length $\xi$. $\xi$ decreases as $\xi\sim \dot\gamma^{-z}$, with $z\cong 0.25$.
In recent simulations [39] for two-dimensional sheared frictionless particles without HIs, a similar exponent was reported for the correlation length extracted from the transverse velocity correlation functions.

§.§ Dynamic coupling between particle rearrangement and solvent dissipation

Coincidence of the large local dissipation $i_2$ and  relative motion of particles:
The left and right panels represent typical snapshots of $i_2$ (top), $I_2$ (middle), and $I_1$ (bottom)
for $\phi=0.64$ at $\dot\gamma=10^{-6}$ and $10^{-4}$, respectively.
In (e) and (f), points for which $|I_1|/\dot{\gamma}<1$ are omitted.

Let us investigate the link between the observed correlated solvent hydrodynamics and the particle motion.
To examine this, we introduce the following approach inspired by the finite element method (FEM) [40].
Delaunay tessellation (triangulation) is implemented for a sufficiently high-density suspension, with each particle position and velocity being identified with those of a corresponding vertex.
For the $e$-th triangle, the vertices are labeled as $i=1,2,3$, and the velocity ${\mbox{\boldmath$U$}}_i^{e} = (U_{x,i}^{e}, U_{y,i}^{e})$ is assigned to the $i$-th vertex.
We quantify the relative particle motion by evaluating the divergence ($I_1^e$) and the rotationally invariant shear gradient ($\sqrt{I_2^e/2}$) of the particle velocities as

\begin{eqnarray}
I_1^e = {\dfrac{\partial}{\partial x}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{x,i}^{e}}\biggr] + {\dfrac{\partial}{\partial y}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{y,i}^{e}}\biggr]
\end{eqnarray}

\begin{eqnarray}
I_2^e = ({D_{xx-yy}^e})^2 + ({D_{xy+yx}^e})^2,
\end{eqnarray}

\begin{eqnarray}
D_{xx-yy}^e &=& {\dfrac{\partial}{\partial x}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{x,i}^{e}}\biggr] - {\dfrac{\partial}{\partial y}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{y,i}^{e}}\biggr], \nonumber \\
D_{xy+yx}^e &=& {\dfrac{\partial}{\partial x}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{y,i}^{e}}\biggr] + {\dfrac{\partial}{\partial y}}\biggl[\sum_{i} N_i^{e}({\mbox{\boldmath$r$}}){U_{x,i}^{e}}\biggr] , \nonumber
\end{eqnarray}

and $N_i^e({\mbox{\boldmath$r$}})$ is known as the shape function in FEM. Please refer to the SI for details of the expression of $N_i^e({\mbox{\boldmath$r$}})$ and derivations of $I_1^e$ and $I_2^e$.

In Fig. <ref>, we plot $I_1$, $I_2$ (based on the particle velocity) and $i_2$ (based on the solvent velocity) at low and high $\dot\gamma$.
The regions with larger values of $I_2$ and $i_2$ are well overlapped, demonstrating that stronger solvent dissipation is associated with larger relative shearing particle motion.
More interesting is that at such high-shear-gradient locations, significant changes in the local particle density, whose characteristic time scale is faster than that of the applied shear flow ($1/\dot\gamma$), are observed.
At a lower $\dot\gamma$, such local density changes occur cooperatively, and their spatial patterns resemble those of the local shear gradients of the solvent and particle velocities.
Due to the incompressibility condition of the solvent component, which is completely satisfied in our simulations, the induced solvent flow has a transverse nature (pure shear or rotation). Such a transverse nature of the solvent dynamics should prevent enhancement of the longitudinal modes of the particle dynamics, whereas we here observe a significant emergence of density changes occurring faster than the time scale of the applied flow ($|I_1|/\dot\gamma \gg1$).
Note that areas of increases and decreases in particle-density changes are alternatively adjacent to each other with distances of several particle sizes, suppressing the density changes at larger length scales. These behaviors of $I_1$ and $I_2$ are confirmed by evaluating their structure factors shown in the SI.
Naively, particle rearrangement in suspensions causes simultaneous  flow of the surrounding solvent, and such a forcibly flowing solvent further pushes other particles, in which flow channels of solvents should be formed by creating gaps between particles. For larger $\dot\gamma$, flow channels with sufficiently wide gaps are created. In contrast, for smaller $\dot\gamma$, the shear stress is not sufficiently large, resulting in narrower gaps, by which flow channels or streamlines are extended to almost the same extent.
In this case, as observed in Fig. 5(e), significant positive and negative particle density changes are alternatively induced in space, indicating that to make flow channels extend over a longer length scale ($\sim\xi$) under narrower average gaps, gap rearrangement against the incompressibility condition is necessary.

§ DISCUSSION AND REMARKS

We note that shear thinning in non-Brownian suspensions is not solely attributed to HIs; notably, in the Stokes regime of hard-sphere suspensions with a Newtonian solvent, shear thinning does not occur [37, 41, 42].
A possible origin of shear thinning could be particle softness or repulsive interactions (effective softness).
In Ref. [7], numerical evidence shows that shear thinning becomes more significant with increasing repulsive-interaction range. This finding is consistent with the results of repulsion-range-control experiments [11]; especially, when particles can be considered almost hard spheres as the repulsion range is negligibly small, shear thinning was hardly observed [11].
In stationary states, a larger repulsion range effectively increases particle size, while under shear flow, an increase in the characteristic force magnitude may reduce the effective particle size [10, 11, 5, 9].
This reduction in the particle size can be anisotropic [44], as particle overlaps (gaps) are more enhanced than in stationary states along the compression (extension) axis.
In our system, HIs with the incompressibility condition should impose strong restrictions on the kinetic paths of particles compared to systems without HIs.
Such restrictions potentially prevent the particle configuration from undergoing relaxation, subsequently leading to an enhancement of particle overlaps (gaps). In this context, we infer that the involvement of HIs plays a crucial role in determining the quantitative aspects of shear-thinning behaviors.

The issues mentioned above will be investigated in detail elsewhere by developing comparative simulations with and without HIs. Presently, let us consider phenomenologically how the macroscopic rheology is related to the steady-state properties of solvent flows and the particle configuration as follows.
The characteristic magnitude of the local shear gradients is set as $\hat{\dot\gamma}$.
Then, from Eq. (<ref>), ${\bar \eta}\dot\gamma^2 \sim \eta_s \hat{\dot\gamma}^2$, and by making use of the relation ${\bar \eta} \sim \dot\gamma^{-\lambda}$, we obtain ${\hat{\dot\gamma}}\sim {\dot\gamma}^{1-\lambda/2}$.
With $h$ and $s$ being the average “gap” and particle sizes, respectively, we may assume that $\hat{\dot\gamma}$ scaled as ${\hat{\dot\gamma}}\sim \dot\gamma s/h$ [45], where $\dot\gamma s$ is the typical relative velocity between neighboring particles.
Consequently, we obtain $h\sim \dot\gamma^{\lambda/2}$, indicating that for larger $\dot\gamma$ associated with larger shear stress, the average gap size becomes wider.
We further assume that the average overlap between particles grows by nearly the same degree as the average gap size and thus that the particle stress almost linearly changes with $h$ [46].
Therefore, with $G$ being a constant with dimensions of the shear modulus, we may obtain ${\bar \eta} \dot\gamma\sim G h/s$.  Substituting ${\bar \eta}\sim \dot\gamma^{-\lambda}$ and $h\sim\dot\gamma^{\lambda/2}$ into ${\bar \eta} \dot\gamma\sim G h/s$ gives $\lambda\sim 2/3$. This value of $\lambda$ is close to the simulation results for smaller $\dot\gamma$ and larger $\phi$.

Next, let us set the structural relaxation period as ${\hat \tau}$:
during ${\hat\tau}$ the average particle configuration undergoes shear deformation with an average shear strain of $\dot\gamma {\hat \tau}$ on the order of $h/s$, and after this period, the configuration may relax.
Thus, ${\hat \tau}\sim h/\dot\gamma s\sim {\hat{\dot\gamma}^{-1}}$, which immediately means that the time scales of the structural relaxation and interstitial solvent flows are almost identical.
Moreover, the relation $\bar\eta \dot\gamma\sim G h/s$ can be rewritten as $\bar\eta\sim G h/\dot\gamma s \sim G{\hat{\dot\gamma}^{-1}}\sim G{\hat\tau}$. Namely, the steady-state viscosity is controlled by $h$ or $\hat\tau$, whose values depend on the details of the relaxation mechanism.
Here, we implicitly assume that the average size of gaps or overlaps, denoted as $h$, is determined through the probability distributions of $i_2$. In the SI, we describe the changes in microstructures with varying $\dot\gamma$ are described for dense case ($\phi=0.64$), but a direct evaluation of $h$ necessitates more careful analysis, including precisely determining the reference when measuring the particle gaps or overlaps.

In this study, we have reported the spatially correlated viscous dissipation in shear-thinning non-Brownian suspensions. Its correlation length $\xi$ decreases with increasing $\dot\gamma$ as $\xi\sim {\dot\gamma}^{-z}$, where $z\cong 0.25$.
This study has also indicated a close link with the thinning of the viscosity as $\eta\sim \dot\gamma^{-\lambda}\sim \xi^{\lambda/z}$.
We have further argued that the thinning of $\eta$ is linked to an increase (decrease) in the size of the gaps (overlaps) between particles, which, however, would not occur in the absence of the particle softness in the first place.
At larger $\dot\gamma$, the gap size is greater, by which, to the same extent, the flow channels are shorter and the local velocity gradients are weaker. Such tuning of the local velocity gradients reduces the macroscopic viscosity.
Then, we have demonstrated that the correlated solvent hydrodynamics are directly coupled with many-body particle motion.
Intriguingly, significant longitudinal particle motion is induced against the solvent incompressibility over the length scale of $\xi$, which may be associated with cooperative creation of gaps and flow channels.
Our preliminary simulation results show a significant difference in the observed viscosities with and without HIs in the present ranges of $\dot\gamma$ and $\phi$.  However, this difference seems less pronounced as $\dot\gamma$ decreases and $\phi$ increases.
More detailed results and discussions of how solvent hydrodynamics influence the macroscopic rheology will be presented elsewhere.

We thank Prof. R. Seto for valuable discussions, particularly regarding the significance of the repulsive-interaction range on the qualitative aspects of rheological behaviors.
This work was supported by KAKENHI (Grants No. 25000002 and 20H05619) and the special fund of Institute of Industrial Science, The University of Tokyo.

§ MICROSTRUCTURES UNDER SHEAR FLOW

$g_{\rm LL}({\mbox{\boldmath$r$}})$ (left) and its cross-sectional  view (right) for $\dot\gamma=10^{-6}$ (a), $10^{-5}$ (b), and $10^{-4}$ (c).
In the right panels the purple and dark green lines represent $g_{\rm LL}({\mbox{\boldmath$r$}})$ along the extension and compression axes, respectively.

Typical snapshots of the force-bonds at $\dot\gamma=10^{-6}$ (a)-(c) and $10^{-4}$ (d)-(f) at $\phi=0.64$ for various threshold values.

In this section,  we examine how the microstructures are modulated by varying the shear rate $\dot\gamma$
through investigations of the partial pair distribution function (PDF) defined for both the same and different species as   [47]

\begin{eqnarray}
g_{\mu\mu}({\mbox{\boldmath$r$}}) = \dfrac{\hat{\mathcal A}}{{\hat N}_\mu ({\hat N}_\mu-1)} {\sum_{i \in \mu}}'{\sum_{j(\ne i)\in \mu}}'\langle \delta({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r$}}_{i}^{(\mu)}+{\mbox{\boldmath$r$}}_{j}^{(\mu)})\rangle,  \nonumber \\ \label{gAA}
\end{eqnarray}

\begin{eqnarray}
g_{\rm SL}({\mbox{\boldmath$r$}}) = \langle \dfrac{\hat{\mathcal A}}{{\hat N}_{\rm S} {\hat N}_{\rm L}} {\sum_{i \in \rm S}}'{\sum_{j \in \rm L}}' \delta({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r$}}_{i}^{\rm (S)}+{\mbox{\boldmath$r$}}_{j}^{\rm (L)})\rangle,  \label{gAB}
\end{eqnarray}

respectively.  Here, ${\mbox{\boldmath$r$}}_{i}^{(\mu)}$ represents the position of the $i$-th $\mu$-species particle ($\mu=\rm L,S$).
In Eqs. (<ref>) and (<ref>), the prime denotes that the summantion is taken over the particles for which $-0.3H \le y_{i}^{(\mu)}\le 0.3H$.
Notice that the walls are located at $y=\pm 0.5H$.
Additionally, ${\hat N}_{\mu}$ is the average number of $\mu$-species particle in the region of $-0.3L\le y \le 0.3L$ and $\hat{\mathcal A}$ denotes the area of this region given by $\hat{\mathcal A}=0.6 {\mathcal A}$.

In the left panels of Fig. <ref>, we show $g_{\rm LL}({\mbox{\boldmath$r$}})$ at $\phi=0.64$ for various values of $\dot\gamma$ in the shear-thinning regime, demonstrating that a structural anisotropy becomes evident with increasing $\dot\gamma$.
Along the extension direction (${\hat x}={\hat y}$), the first peak of $g_{\rm LL}({\mbox{\boldmath$r$}})$ becomes smaller and more broad towards the outer region with an increase in $\dot\gamma$, indicating that dilution occurs.
On the other hand, along the compression axis, almost the opposite occurs.
Such behaviors are more clearly shown in the right panels of Fig. <ref>, showing the cross-sectional views of $g_{\rm LL}({\mbox{\boldmath$r$}})$ along the extension and compression axes. Note that essentially the same behaviors are observed for $g_{\rm SS}({\mbox{\boldmath$r$}})$ and $g_{\rm SL}({\mbox{\boldmath$r$}})$, which are not shown here.

These observations in the partial PDF are supplemented by measurements of the force bond.
We represent the interaction force magnitude between the $i$- and $j$-th particles as $|f_{ij}|=|{\partial U}/{\partial {{\mbox{\boldmath$r$}}}_{ij}}|$ as the magnitude of the interaction forces between $i$- and $j$-th particles.
In Figs. <ref>(a) and (d), we display typical snapshots of patterns of the force bonds, where the pairs of particles with non-zero interactions ($|f_{ij}>0|$) are described by connecting bonds, for $\dot\gamma=10^{-6}$ and $10^{-4}$, respectively, at $\phi=0.64$.
Here, the color of each bond indicates the scaled value of $|f_{ij}|$ by $2\pi s_{\rm S} P_p$, with $P_p$ being the particle pressure defined as [47]

\begin{eqnarray}
P_p = -\dfrac{1}{2\mathcal A}\sum_i\sum_{j>i} { {\mbox{\boldmath$r$}}}_{ij}\cdot \dfrac{\partial U}{\partial {{\mbox{\boldmath$r$}}}_{ij}}.
\end{eqnarray}

For the interaction potential used in this study, as given by Eq. (<ref>), although the potential is steep, the cutoff distance is slightly longer than usual.
Consequently, when we interpret particle pairs experiencing a non-zero interaction force as being in contact, most particles would be nearly in a perfect contact state for sufficiently large value of $\phi$.
Nevertheless, this interpretation is oversimplified.
It is crucial to recognize that the characteristic force magnitude varies with the shear rate. As the shear rate increases, so does the driving external force, leading to a greater characteristic force. Therefore, if the interaction force between particles is significantly smaller than this characteristic force, these pairs can be considered ineffective in contact.

At this stage, identifying an exact threshold value for judging the interaction force as relevant for contact or not is still difficult.
Instead, by adjusting the threshold value, bonds with force magnitudes below this threshold are not displayed, and subsequently, we assess how their overall patterns change with the shear rate.
In Figs. <ref> (b) and (c), at \(\dot{\gamma} = 10^{-6}\), we display the same snapshots as (a) but with different threshold values.
Furthermore, similar adjustments are made to (d), resulting in (e) and (f).
A noticeable difference between the cases of $\dot\gamma=10^{-6}$ and $10^{-4}$ is observed even for the threshold values of $|f_{ij}|/(2\pi s_{\rm S} P_p)=0.1$,
It is observed that more bonds (mainly along the extension axis) disappear at higher shear rates:
While larger forces are anisotropically exerted along the compression axis, forces along the extension axis are weaker, indicating an effective decrease in contact (or effective opening of gaps) along the elongation axis with an increase in the shear rate.

Before closing this section, it is noteworthy to mention the following point.
Although we here scale the interaction forces by $2\pi s_{\rm S} P_p$, the characteristic force magnitude may be alternatively determined by particle shear stress. Notably, a similar trend was found when the force is scaled by the shear stress.

§ QUANTIFYING RELATIVE PARTICLE MOTION

(a) Snapshot of $I_1$. Each circle indicates the position of a particle, the color indicates the particle size, and the arrow indicates the velocity. Schematic pictures of $I_1$ (b) and $I_2$ (c). $I_1$ and $I_2$ measure the degree of divergent and shearing motions of the triangular element, respectively. The blue arrows schematically illustrate typical accompanying flows of the surrounding solvent.

We explain the details of our analysis of relative particle motion. To quantify the relative motion, we define two quantities, $I_1$ and $I_2$, whose physical meanings are schematically shown in Fig. <ref>: $I_1$ directly measures the rate of change of the triangular area, while $I_2$ measures the degree of shear deformation rate of the triangle. These quantities are calculated for each element.
We perform Delaunay triangulation, in which the vertices are the positions of the particles ${\mbox{\boldmath$R$}}_i$.
The solid lines in the left panel of Fig. <ref> show the partitioning.
The quantities $I_1$ and $I_2$ are tied to each triangular element. In this way, we can regard $I_1$ and $I_2$ as spatial variables.

(Left) A first-order triangular element.  On each element, the approximate function ${\tilde f}(x,y)$ is described as a plane (a). The plane function ${\tilde f}(x,y)$ can be expressed as a product of shape functions $N_i^e(x,y)$ within the element and the values at the vertices $f_i^e$ (b-d).

Here, we define the two quantities $I_1$ and $I_2$ as invariants of a tensor.
The derivation follows the simple procedure of the conventional finite element method (FEM) [40]. We pick up a triangular element $e$, whose vertices are numbered as $1,2$, and $3$ (see Fig. <ref>).
We let $f^e_i$ denote the value of an arbitrary variable ${\tilde f}$ at the vertex $i$ ($i=1,2,3$) on the element $e$.
The vertex positions are denoted as $(x^e_1, y^e_1)$, $(x^e_2,y^e_2)$, and $(x^e_3,y^e_3)$.
Then, we formally express $f_i^e$ as

\begin{align}
\mqty(f^e_1\\f^e_2\\f^e_3) &= \mqty(1 & x^e_1 & y^e_1 \\ 1&x^e_2&y^e_2 \\ 1&x^e_3&y^e_3)\mqty(\alpha^e_1\\\alpha^e_2\\\alpha^e_3) =: S\mqty(\alpha^e_1\\\alpha^e_2\\\alpha^e_3),
\end{align}

where the parameters $\mqty(\alpha^e_1&\alpha^e_2&\alpha^e_3)$ are given by

\begin{align}
\mqty(\alpha^e_1\\\alpha^e_2\\\alpha^e_3)
= S^{-1}\mqty(f^e_1\\f^e_2\\f^e_3)
\end{align}

\begin{eqnarray}
S^{-1} &=&  \mqty( a^e_1 & a^e_2 & a^e_3 \\ b^e_1 & b^e_2 & b^e_3 \\ c^e_1 & c^e_2 & c^e_3 )  \nonumber \\
&:=&\frac{1}{2A_e} \mqty( x^e_2y^e_3 - x^e_3y^e_2& x^e_3y^e_1 - x^e_1y^e_3&x^e_1y^e_2 - x^e_2y^e_1\\
y^e_2 - y^e_3&y^e_3 - y^e_1&y^e_1 - y^e_2\\
x^e_3 - x^e_2&x^e_1 - x^e_3&x^e_3 - x^e_2), \nonumber\\
\end{eqnarray}

where $A_e=\det(S)/2$.
In the linear FEM, by defining the shape function as $N^e_i(x,y) = a^e_i + b^e_ix + c^e_iy\quad(i=1,2,3)$, we express the field variable ${\tilde f}$ on the element $e$ by interpolating the values assigned to the vertices as

\begin{align}
\tilde{f}(x,y) &= \sum_{i=1,2,3} N_i^e({{\mbox{\boldmath$x$}}})f_i^e \nonumber \\
&=\mqty(N^e_1(x,y) & N^e_2(x,y) & N^e_3(x,y))\mqty(f^e_1\\f^e_2\\f^e_3),   \label{f_tilde}
\end{align}

where ${{\mbox{\boldmath$x$}}}=(x,y)\in e$.

Let $(V_{xi}^e, V_{yi}^e)$ be the velocity of the $i$-th particle on the element $e$.
According to the above explained procedure, we can approximate the particle velocity `'field” at an arbitrary position ${{\mbox{\boldmath$x$}}}$ on the element $e$ as the linear interpolation of $(V_{xi}^e, V_{yi}^e)$:

\begin{eqnarray}
V_x^e({{\mbox{\boldmath$x$}}})&=& \sum_{i=1,2,3} N_i^e({{\mbox{\boldmath$x$}}})V_{xi}^e,  \\
V_y^e({{\mbox{\boldmath$x$}}})&=& \sum_{i=1,2,3} N_i^e({{\mbox{\boldmath$x$}}})V_{yi}^e.
\end{eqnarray}

Then, we can compute the deformation rate as

\begin{eqnarray} \mqty(D_{xx}^e & D_{xy}^e \\ D_{yx}^e & D_{yy}^e ) := \mqty(\frac{\partial}{\partial x}{{V}_x^e({{\mbox{\boldmath$x$}}})} &\frac{\partial}{\partial x}{{V}_y^e({{\mbox{\boldmath$x$}}})}\\ \frac{\partial}{\partial y}{{V}_x^e({{\mbox{\boldmath$x$}}})} &\frac{\partial}{\partial y}{{V}_y^e({{\mbox{\boldmath$x$}}})}),
\end{eqnarray}

which further defines

\begin{eqnarray}
I_1^e &:=& D_{xx}^e + D_{yy}^e,  \\
I_2^e &:=& (D_{xx}^e - D_{yy}^e)^2 + (D_{xy}^e + D_{yx}^e)^2.
\end{eqnarray}

$S_{I_1}(k)/k^2$ against $ks_{\rm S}/2\pi$ at $\phi=0.41$ (a), $0.54$ (b), $0.64$ (c), and $0.68$ (d).
At a lower $k$, enhancement of $S_{I_1}(k)$ is suppressed and behaves as $k^a$ ($a\sim 2$).

Structure factor $S_{I_2}(k)$ of $I_2$ against $ks_{\rm S}/2\pi$ at $\phi=0.41$ (a), $0.54$ (b), $0.64$ (c), and $0.68$ (d).
For larger $\phi$,  $S_{I_2}(k)$ at lower $k$ grows with decreasing $\dot\gamma$.
Such behaviors are qualitatively similar to those of the structure factor of $i_2$ shown in Fig. 4 in the main text.

To quantify the spatial correlation of the transverse and longitudinal particle motions,  we define the following structure factors:

\begin{eqnarray}
S_{I_1}( { {\mbox{\boldmath$k$}}} ) &=&  \dfrac{1}{\mathcal A}\langle |{\hat I}_{1}(\mbox{\boldmath$k$})|^2 \rangle,\\
S_{I_2}( { {\mbox{\boldmath$k$}}} ) &=&  \dfrac{1}{\mathcal A}\langle |{\hat I}_{2}(\mbox{\boldmath$k$})|^2 \rangle,
\end{eqnarray}

where ${\hat I}_{m}(\mbox{\boldmath$k$})=\sum_{e} A_e {\exp}({-i{\mbox{\boldmath$k$}}\cdot{\mbox{\boldmath$R$}_e^{G}}}) I_m^e$ ($m$=1,2) and ${\mbox{\boldmath$R$}_e^{G}}=(1/3)\sum_{i=1,2,3}{\mbox{\boldmath$R$}_i}$ is the geometrical center of the $e$-th triangular element.
Note that, although a window function is not used here, the resultant structure factor for $k>2\times 2\pi/L$ is hardly altered by whether implementing it or not.
In the continuum limit, it may be reduced to ${\hat I}_{m}(\mbox{\boldmath$k$})=\int d{\mbox{\boldmath$r$}} {\exp}({-i{\mbox{\boldmath$k$}}\cdot{\mbox{\boldmath$r$}}}) w({\mbox{\boldmath$r$}}) I_m({\mbox{\boldmath$r$}})$.
As denoted in the main text, anisotropy and localization near the walls of the particle motions are not remarkable for the present ranges of $\phi$ and $\dot\gamma$, and thus we may safely replace $S_{I_m}( { {\mbox{\boldmath$k$}}})$ by its angle average $S_{I_m}(k)$.

Figures <ref> and <ref> display $S_{I_1}(k)$ and $S_{I_2}(k)$, respectively.
As demonstrated in Fig. 5 in the main text, the regions with larger values of $I_2$ and $i_2$ almost correspond to each other, which indicates that stronger solvent dissipation is associated with larger relative shearing particle motion.
This correspondence is further supported by Fig. <ref>: the spatial correlation of $I_2$ and its intensity grow as $\phi$ increases and $\dot\gamma$ decreases. Such tendencies of $S_{I_2}(k)$ are qualitatively similar to those of $S(k)$, the structure factor of $i_2$, which is shown in Fig. 4 in the main text.
Unlike $S_{I_2}(k)$, Fig. <ref> shows that an enhancement of $S_{I_1}(k)$ at the lower $k$ is highly suppressed: as indicated in the main text, particle density increases and decreases occur alternatively in space at several particle length scales, which may suppress $I_1$ at larger scales.
Since the rearrangement dynamics of particles are inevitably accompanied by solvent flows surrounding them, it is necessary to form flow channels by enhancing the particle gaps. However, the incompressibility of the solvent components tends to prevent particle density changes. The obtained form of $S_{I_1}(k)$ may reflect a compromise of these contradicting requirements.

[1]
J. Mewis and N.J. Wagner,
Colloidal Suspension Rheology,
(Cambridge University Press, New York, 2012).

[2]
J. J. Stickel and R. L. Powell,
Fluid mechanics and rheology of dense suspensions,
Annu. Rev. Fluid Mech. 37, 129 (2005).

[3]
M. M. Denn and J. F. Morris,
Rheology of Non-Brownian Suspensions,
Annu. Rev. Chem. Biomol. Eng. 5, 203 (2014).

[4]
R.I. Tanner,
Review Article: Aspects of non-colloidal suspension rheology  featured,
Physics of Fluids 30, 101301 (2018).

[5]
C. Ness, R. Seto, and R. Mari,
The Physics of Dense Suspensions,
Annu. Rev. Condens. Matter Phys.  13, 97 (2022).

[6]
R. Seto, R. Mari, J. F. Morris, and M. M. Denn,
Discontinuous Shear Thickening of Frictional Hard-Sphere Suspensions,
Phys. Rev. Lett. 111, 218301 (2013).

[7]
R. Mari, R. Seto, J. F. Morris, and M. M. Denn,
Discontinuous shear thickening in Browniansuspensions by dynamic simulation,
Proc. Natl. Acad. Sci. U.S.A.  112, 15326 (2015).

[8]
M. Wyart and M. E. Cates,
Discontinuous shear thickening without inertia in dense non-brownian suspensions, Phys. Rev. Lett. 112, 098302 (2014).

[9]
G. Chatté, J. Comtet, A. Niguès, L. Bocquet, A. Siria, G. Ducouret, F. Lequeux, N. Lenoir, G. Overlez, and A. Colin, Shear thinning in non-Brownian suspensions,
Soft Matter 14, 865 (2018).

[10]
I. M. Krieger,
Rheology of monodisperse latices,
Adv. Colloid, Interface Sci. 3, 111 (1972).

[11]
B. J. Maranzano and N. J. Wagner,
The effects of interparticle interactions and particle size on reversible shear thickening: Hard-sphere colloidal dispersions,
J. Rheol. 45, 1205 (2001).

[12]
R. Mari, R. Seto, J. F. Morris, and M. M. Denn,
Shear thickening, frictionless and frictional rheologies in non-brownian
suspensions, J. Rheol. 58, 1693 (2014).

[13]
L. Lobry, E. Lemaire, F. Blanc, S. Gallier, and F. Peters,
Shear thinning in non-Brownian suspensions explained by variable friction between particles,
J. Fluid Mech. 860, 682 (2019).

[14]
A. Vázquez-Quesada, R. I. Tanner, and M. Ellero,
Shear thinning of noncolloidal suspensions,
Phys. Rev. Lett. 117, 108001 (2016).

[15]
A. Vázquez-Quesada, A. Mahmud, A. Dai, M. Ellero, and R. I. Tanner,
Investigating the causes of shear-thinning in non-colloidal suspensions: experiments and simulations,
J. Non Newtonian Fluid Mech 248, 1 (2017).

[16]
R. I. Tanner, C. Ness, A. Mahmud, S. Dai,  and J. Moon,
A bootstrap mechanism for non-colloidal suspension viscosity,
Rheol. Acta 57, 635 (2018).

[17]
M. Kroupa, M. Soos, and J. Kosek,
Slip on a particle surface as the possible origin of shear thinning in non-Brownian suspensions,
Phys. Chem. Chem. Phys. 19, 5979 (2017).

[18]
A. Vázquez-Quesada, A., P. Español, and M. Ellero,
Apparent slip mechanism between two spheres based on solvent rheology: theory and implication for the shear thinning of non-Brownian suspensions,
Physical Review Fluids, 3, 123302 (2018).

[19]
J. Comtet, G. Chatté, A. Niguès, L. Bocquet, A. Siria, and A. Colin,
Pairwise frictional profile between particles determines discontinuous shear thickening transition in non-colloidal suspensions,
Nat. Commun. 8, 15633 (2017).

[20]
E, Lemaire, F. Blanc, C. Claudet, S. Gallier, L. Lobry, and F. Peters,
Rheology of non-Brownian suspensions: a rough contact story,
Rheo. Acta, 62, 253 (2023).

[21]
A. Furukawa and H. Tanaka,
Key Role of Hydrodynamic Interactions in Colloidal Gelation,
Phys. Rev. Lett. 104,  245702 (2010).

[22]
S. Yamanaka, A. Furukawa, and H. Tanaka
Complex dynamical interplay between solid particles and flow in driven granular suspensions,
Phys. Rev. E 100, 012907 (2019).

[23]
A. V. Nguyen Le, Adrien Izzet, and G. Ovarlez, and A. Colin,
Solvents govern rheology and jamming of polymeric bead suspensions,
J. Colloid Interface Sci. 629, 438 (2023).

[24]
Y. Lin, P. Lin, Y. Wang, J. Chen, Z. He, T. P ahtz, N. Phan-Thien,
Shear thinning of non-Brownian suspensions and its variation at different ambient conditions, Physics of Fluids 35, 023301 (2023).

[25]
B. Bernu, J. P. Hansen, Y. Hiwatari, and G. Pastore, Soft-sphere model for the glass transition in binary alloys: Pair structure and self-diffusion,
Phys. Rev. A 36, 4891 (1987).

[26]
R. Yamamoto, and A. Onuki, Dynamics of highly supercooled liquids: Heterogeneity, rheology, and diffusion,
Phys. Rev. E  58, 3515 (1998).

[27]
Y. Nakayama and R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersions, Phys. Rev. E. 71, 036707 (2005).

[28]
J. J. Molina, K. Otomura, H. Shiba, H. Kobayashi, M. Sano, and R. Yamamoto,
Rheological evaluation of colloidal dispersions using the smoothed profile method: formulation and applications, J. Fluid Mech. 792, 590 (2016).

[29]
R. Yamamoto, J. J. Molina, and Y. Nakayama,  Smoothed profile method for direct numerical simulations of hydrodynamically interacting particles,
Soft Matter, 17, 4226 (2021).

[30]
L.D. Landau and E.M. Lifshitz,  Fluid Mechanics
(Pergamon Press, New York, 1959).

[31]
A.J.C. Ladd,
Short-time motion of colloidal particles - numerical-simulation via a fluctuating lattice-Boltzmann equation.
Phys. Rev. Lett. 70, 1339 (1993).

[32]
M. E. Cates, K. Stratford, R. Adhikari, P. Stainsell, J.-C. Desplat,  I. Pagonabarraga, and A.J. Wagner, Simulating colloid hydrodynamics with lattice Boltzmann methods.
J. Phys.: Condens. Matter 16, S3903 (2004).

[33]
H. Tanaka and T. Araki,
Simulation Method of Colloidal Suspensions with Hydrodynamic Interactions: Fluid Particle Dynamics,
Phys. Rev. Lett. 85, 1338 (2000).

[34]
A. Furukawa, M. Tateno, and H. Tanaka,
Physical foundation of the fluid particle dynamics method for colloid dynamics simulation, Soft Matter, 14, 3738 (2018).

[35]
P. Español and P. Warren,
Statistical Mechanics of Dissipative Particle Dynamics,
Europhys. Lett. 30, 191 (1995).

[36]
S. JI, R. Jiang, R.G. Winkler, and G. Gompper,
Mesoscale hydrodynamic modeling of a colloid in shear-thinning viscoelastic fluids under shear flow,
J. Chem. Phys. 135 134116 (2011).

[37]
I.M. Krieger, A Dimensional Approach to Colloid Rheology,
Trans. Soc. Rheol. 7, 101 (1963).

[38]
When an isolated circular particle is put under simple shear flow, $i_2$ around the particle exhibits four-fold symmetry, creating square-like patterns around the particle domain. Even in semi-dilute regions, the $i_2$ field is supposed to be approximated as a superposition of contributions from each particle, resulting in square-like (dark-colored) patterns around the particle regions. This effect is emphasized by the color tones selected for Fig. 2.

[39]
Peter Olsson and S. Teitel,
Dynamic length scales in athermal, shear-driven jamming of frictionless disks in two dimensions,
Phys. Rev. E 102, 042906 (2020).

[40]
T. J. R. Hughes,
The finite element method: linear static and dynamic finite element analysis
(Courier Corporation, 2012).

[41]
A. Lemaître, J.N. Roux, and F. Chevoir,
What do dry granular flows tell us about dense non-Brownian suspension rheology?,
Rheol. Acta 48, 925 (2009)

[42]
F. Boyer, E. Guazzelli, and O. Pouliquen,
Unifying Suspension and Granular Rheology,
Phys. Rev. Lett. 107, 188301 (2011).

[43]
Particle overlaps (gaps) are more enhanced than in stationary states along the compression (extension) axis, whereby the particle-size reduction may occur anisotropically [44].

[44]
A. Furukawa,
Quantification of the volume-fraction reduction of sheared fragile glass-forming liquids and its impact on rheology,
Phys. Rev. Research 5, 023181 (2023).

[45]
P. Mills and P. Snabre,
Apparent viscosity and particle pressure of a concentrated suspension of non-Brownian hard spheres near the jamming transition,
Eur. Phys. J. E 30, 309 (2009).

[46]
Increasing $\dot\gamma$ enhances particle overlap mainly along the compression axis ($\hat x=-{\hat y}$). At the same time, mainly along the elongation axis (${\hat x}={\hat y}$), the gap between particles is also increased to a similar degree.
Note, however, that our preliminary simulations show that the distortion of the particle configuration is not exactly symmetric along the compression and elongation axes and that the average interparticle distance becomes smaller with increasing $\dot\gamma$.

[47]
J-P. Hansen and I.R. McDonald, Theory of Simple Liquids (Academic Press, London and New York, 1976).

[40]
T. J. R. Hughes,
The finite element method: linear static and dynamic finite element analysis
(Courier Corporation, 2012).


# LLM2Loss: Leveraging Language Models for Explainable Model Diagnostics

Shervin Ardeshir
Netflix
<EMAIL_ADDRESS>

###### Abstract

Trained on a vast amount of data, Large Language models (LLMs) have achieved
unprecedented success and generalization in modeling fairly complex textual
inputs in the abstract space, making them powerful tools for zero-shot
learning. Such capability is extended to other modalities such as the visual
domain using cross-modal foundation models such as CLIP[13], and as a result,
semantically meaningful representation are extractable from visual inputs.

In this work, we leverage this capability and propose an approach that can
provide semantic insights into a model’s patterns of successes, failures, and
biases. Given a black box model, its training data, and task definition, we
first calculate its task-related loss for each data point. We then extract a
semantically meaningful representation for each training data point (such as
CLIP embeddings from its visual encoder) and train a lightweight model which
maps this semantically meaningful representation of a data point to its task
loss. We show that an ensemble of such lightweight models can be used to
generate insights on the performance of the black-box model, in terms of
identifying its patterns of failures and biases.111Disclaimer: we used
Microsoft Bing as an assistive tool in writing this manuscript.

##  1 Introduction

Figure 1: An overview of the proposed approach. Given a black-box model ($f$)
trained on a task ($f(x)\rightarrow y$), we aim to estimate the function
$h(g(.))$ which estimates its loss on task $f$, i.e. $l^{f}(f(x),y)$, such
that $g(x)$ is an embedding compatible with a language model. As a result, at
inference time, we can feed textual descriptions and derive insights into the
model’s performance on such data points.

Large language models learn to generate natural language texts based on a
large corpus of textual data. They have shown promise in modeling the semantic
space, which is the representation of meaning and concepts in natural
language. Using natural language prompts, which are inputs that guide the
model to complete a certain task or generate a certain output, large language
models can perform various natural language processing tasks without requiring
any task-specific training data or labels. This is also known as zero-shot
learning, where the model leverages its general knowledge and linguistic
skills to infer the desired output from the prompt. Built on-top-of, or in-
conjunction-with language models, vision-to-text models such as CLIP are able
to map visual space to the same semantic space as natural language by
objective functions such as cross-modal contrastive[13]. This means that such
models learn to align visual and language representations by maximizing their
similarity for pairs of images and texts that are semantically related and
minimizing their similarity for pairs that are not [13]. By doing so, vision-
to-text models can achieve a better understanding of both modalities and
perform various tasks such as image captioning, image retrieval, visual
question answering, etc. As a result, such language based models [4, 13] are
capable of producing highly structured embeddings spaces, which also are
inherently semantic, and thus, human interpretable/understandable. Here we aim
to leverage this characteristic to derive some insight from the performance of
a black-box machine learning model, in terms of understanding it’s potential
areas of strengths, failures, and biases post-hoc and in absence of any given
ground-truth semantic categories.

Given a black box model and its training data, we first calculate its loss for
each data point using a single forward pass. We then extract a semantically
meaningful, and Language-model-compatible representation for each training
datapoint (such as CLIP). We then train a lightweight diagnosis model which
maps the semantically meaningful representation of a datapoint to its task
loss. We show that such a model can provide semantic insight in terms of
identifying patterns of failures and biases in the performance of models.

##  2 Related Work

Language models have shown promise in capturing fairly complicates structures
within the semantic space [4, 18, 5]. On the other hand, understanding of
model robustness and estimating it across an abstract and contineous
embeddinng space has been explored in [2]. Here we aim to add the capabilities
of LLMs to that domain and allow for estimation of robustness in a way that is
semantically interpretable. This is a fairly new area, and relevant to this
line of work may be semantic adversarial approaches [8, 12, 17] which rely on
generative models in order to generate adversarial examples, and [6] where
interpretable controls are discovered for generative models. However, our work
explores a complementary direction, and simpler approach to understanding a
models failures and biases without the need for generating any artificial
input data. Another line of related work which would be relevant to ours is
that of [3, 14, 11] where the authors aim at using latent spaces to improve
model robustness.

Zeros shot capability of CLIP [13] embeddings have shown promise in many
applications such as manipulation of NERFs [19], avatars [7], shape generation
[16], and debiasing [3].

##  3 Approach

An overview of the approach is shown in Figure 1. We assume that we are given
a black-box model trained to map input $\mathbf{X}=\\{x_{i}\\}_{i=1,...,N}$ to
output $\mathbf{y}=\\{y_{i}\\}_{i=1,...,N}$ via function
$f:\mathbf{x}\rightarrow\mathbf{y}$ through minimization of loss function
$l^{f}(f(x),y)$. In our experiments, we set $f$ to be a face recognition model
trained on the CelebA dataset[10]. We also assume that the format of $x_{i}$s
is compatible with some Language compatible model $g$. By language compatible,
we mean that there is a foundation language-based model that can map this
input to a common space with language models. Given that the input data for
the face recognition model is visual, in our experiments we used CLIP visual
encoder as the language-compatible encoder $g$, and CLIP textual encoder as
$\tilde{g}$.

###  3.1 Training

Our goal is to train model $h$, such that it can predict model $f$’s loss for
datapoint $x_{i}$, i.e. $l^{f}(f(x_{i}),y_{i})$, solely from the content of
its input $x_{i}$. Given that loss values are often continuous, we experiment
with two variations of this mapping functions, namely, a regression model, and
a ranking model. We trained both regressor and ranker alternatives in an
ensemble setup to be able to have uncertainty on each of the generated
insights. We trained our ensemble models in a very simple setup, and by
training 5 independent runs with randomized initializations.

####  3.1.1 $h$ as a Loss Regressor

The regression model is trained using MSE loss and defined as:

$l_{regr}^{h}=MSE(h(g(x_{i})),l^{f}(f(x_{i}),y_{i})),$ (1)

in other words, $h$ estimates the absolute value of the task loss.
Intuitively, this would estimate how much error model $f$ is expected to have
for a given input.

####  3.1.2 $h$ as a Loss Ranker

An alternative is to train $h$ using a ranking loss. In this setup we sample
pairs of datapoints $(x_{i},x_{j})$, and train $h$ to be able to predict which
datapoint leads to a higher loss through $f$. We train this model using a
binary-crossenropy loss defined as:

$l_{rank}^{h}=BCE(\Delta h(g(x)),\Delta l^{f}(f(x),y)>0),$ (2)

in which:

$\Delta h(g(x))=h(g(x_{i}))-h(g(x_{j})),$ (3)

and:

$\Delta l^{f}(f(x),y)=l^{f}(f(x_{i}),y_{i})-l^{f}(f(x_{j}),y_{j}).$ (4)

Intuitively, this model predicts if model $f$ does better for datapoint
$x_{j}$ compared to datapoint $x_{i}$. In our setup, the face recognition
model was trained using a multi-class classification setup (predicting
identity as a class label from face image content), and therefore the per-
datapoint loss $l^{f}(f(x_{i}),y_{i})$ is the output of a categorical cross-
entropy loss function, and therefore a continuous value.

####  3.1.3 $h$ as a Classifier

Another aspect of trying to understand a model’s strengths and weaknesses is
to understand whether the model has "seen" a specific type of input during
training or not. In other words, is input $\mathbf{x}$ in-distribution when
compared to it’s training dataset. A datapoint being seen or unseen can have a
significant impact on a models predictions for this datapoint being
meaningful. More often than not, such predictions are not meaningful. As a
result, we also believe that the estimated values for the losses may or may
not be valid for a datapoint depending on such phenomenon. Thus, we experiment
with equipping our diagnostics to a secondary output which indicates if a
datapoint is in or out of distribution. We model such an output as a binary
classification problem: in-distribution (i.d.) vs out-of-distribution
(o.o.d.), or seen vs unseen. Given that each training datapoint results in a
512 dimesnional CLIP embedding, we generate random 512 dimensional vectors as
out-of-distribution areas of the embedding space. We then train a binary
classifier on in-distribution vs out-of-distribution embeddings. Our
assumption is that we only have access to the training datapoints of the task
of interest, and thus we can only use random negatives. Our experiments
(included in Section 4.1.1 and Figure 2) indicate that such a simple approach
achieves meaningful results in terms of allowing a user to understand what
type of datapoints have been seen by a model. Also, our experiments suggest
that the loss-regressor (explained in Section 3.1.1) already has this
capability and is strongly capable of OOD detection, which in fact suggests
that there may not be a need for a separate OOD detection component. However,
we still experimented with this variation to provide a good point of
comparison to the OOD detection capabilities of the loss-regressor and loss-
ranker models.

###  3.2 Inference

As shown at the bottom row of Figure 1, at inference time, we aim to provide
some insights on a given textual input $x_{text}$. Intuitively, given textual
input $x_{text}$, $-h(\tilde{g}(x_{text}))$, will estimate, how good the
model’s performance is on $x_{text}$ type datapoints. As an example, we can
feed textual inputs such as $x_{text}$="face image of a long-haired smiling
male wearing glasses" and get an estimated loss for that. Assuming
generalization, this is an estimate of how good the black-box model $f$ is on
its own task, on inputs that resemble $x_{text}$. Given that absolute values
of a loss are often hard to interpret, a better way would be to feed two
different textual inputs and compare their estimated performance. As an
example feeding "face of a smiling long-haired male", and "face of a smiling
long-haired female", and observe if the model’s performance has a bias towards
the former or the latter. In the OOD classifier variation of $h$, the same
analysis would tell us how familiar a model is with the inputs that resemble
$x_{text}$, in other words, if semantic concept $x_{text}$ is in-distribution
with respect to $f$’s training data.

##  4 Experiments

In our experiments, we used celebA dataset[10] and used a trained face-
recognition model as $f$. We use categorical cross-entropy loss of the face
recognition model as $l^{f}$, and the target for $h$, and experimented with
simple linear and (2-layer) non-linear MLPs for the architecture of $h$, for
both variations of regressor and ranker. Given that we used CLIP [13] encoders
as $g$ and $\tilde{g}$, $g(x_{i})$ is a 512-dimensional vector resulting from
the image encoder of CLIP. In inference mode, we use the text encoder of CLIP
as $\tilde{g}$. The main application of our proposed approach in realistic
scenarios would be better demonstrated using qualitative examples. Thus we
demonstrate a few examples of how we can use the trained model h to derive
insights on model bias and failures. We then aim to measure the correlation
between our generated insights to their true values based on attributes
available for the celebA [10] dataset.

###  4.1 Loss/Prediction Validity: Explainable/Semantic OOD Detection

As mentioned earlier, an important aspect of understanding a models strengths
and weaknesses is to understand on which datapoints a models predictions are
valid. In other words, whether a datapoint is in-distribution or out-of-
distribution for a model. The classifier variation of $h$ is explicitly
trained to do so. In Section 4.1.1 we visualize how insights about the
validity of predictions of a model could be derived using this variation pf
$h$, and in Section 4.1.2 we quantify the performance of different variations
of $h$ in terms of their capability in performing OOD detection in the
semantic space (more details in Section 4.1.2).

####  4.1.1 Qualitative

Here we visualize how insights about the validity of predictions of a model
could be derived using the classifier variation of $h$. Figure 2 shows the
output of prediction of $h_{classifier}$ on a set of textual inputs. Right
shows the output of the classifier post the Sigmoid activation, and left shows
the output prior to it, to signify the differences. As it can be observed,
textual inputs such as face, frontal-face, and frontal face headshot are
predicted to be more in-distribution compared to concepts such as cat, dog,
building. Given that model $f$ in our experiments is a face recognition model
trained on celebA it’s only been exposed to face images. Thus such textual
inputs are correctly predicted to be more in-distribution / familiar to the
model. Interestingly, even within different variations of textual input
containing the word face, the ones containing the phrase "frontal face" has
much higher scores than the one containing "side face". Knowing that celebA
primarily contains frontal face headshots, our anecdotal qualitative examples
seem to capture the seen vs unseen domains fairly meaningfully. We argue that
the reason why the post Sigmoid probabilities are all close to 1 is that in
our very simple training setup ($x_{tr}$ vs. $x_{rand}$), the negative samples
are randomly generated. As a result, the classifier can draw the boundary
between random vectors and a manifold containing all natural images, resulting
in all natural textual inputs being classified as seen/in-distribution. We
argue that this is something that can easily improved once out-of-distribution
samples come from natural images/inputs and not randomly generated. However,
our experiments indicate that even without that and with this extremely simple
way of randomly generated negatives, the model meaningfully learns that faces
are more familiar to the model compared to concepts such as cars,dog, etc. We
later quantify this capability of $h_{classifier}$ and compare it to
$h_{regressor}$ and $h_{ranker}$, in section 4.1.2.

Figure 2: Explainable OOD detection. Right: post-sigmoid. Left: pre-sigmoid.
It can be observed that phrases describing face data are correctly identified
as better descriptions of the data that the model has seen.

####  4.1.2 Quantitative

Table 1 shows the performance of our different variations of $h$ on semantic
OOD detection task. In this task, we feed semantic inputs describing the seen
dataset (such as "face", "portrait", "human face", etc) in addition to out-of-
distribution data. We then feed class labels of some out-of-distribution /
unseen datasets as $x_{text}$ to the inference pipeline
$h(\hat{g}(x_{text}))$. We then evaluate whether the in-distribution class
labels are retrieved correctly and with a higher score, i.e. measure the AUROC
of $h(\hat{g}(x_{text}))$. Rows 2,3, and 4 of Table 1 reflects that metric for
classes of CIFAR10 [9], and super-classes and classes of CIFAR100 [9]. It can
be observed that all models perform better than chance, and interestingly, the
loss-regressor variation (i.e. $h_{reg}$) performs better than the
classification variation (i.e $h_{class}$) which is directly trained to
perform this task. As a result, the loss regressor model in its ensemble form
performs very reliably not only in terms of estimating the loss (as it’s
directly trained on it), but also the validity of it.

OOD dataset/model | $h_{rank.}$ | $h^{ens-mean}_{rank.}$ | $h_{reg.}$ | $h^{ens-mean}_{reg.}$ | $h_{class.}$ | $h^{ens-mean}_{reg.}$
---|---|---|---|---|---|---
random | $0.670\pm 0.217$ | 0.686 | $1.0\pm 0.0$ | 1.0 | $1.0\pm 0.0$ | 1.0
CIFAR10-classes | $0.576\pm 0.049$ | 0.62 | $0.919\pm 0.047$ | 1.0 | $0.848\pm 0.009$ | 0.840
CIFAR100-super-classes | $0.555\pm 0.058$ | 0.578 | $0.951\pm 0.036$ | 0.989 | $0.903\pm 0.012$ | 0.915
CIFAR100-classes | $0.700\pm 0.101$ | 0.722 | $0.951\pm 0.016$ | 0.978 | $0.878\pm 0.005$ | 0.88
Table 1: Loss-validity, i.e. explainable OOD Detection. In-distribution
dataset is always celebA. Left column specified classes from what out-of-
distribution dataset were used in this evaluation.

###  4.2 Loss Estimation:

The regressor and ranker models aim directly at estimating the absolute value
and ranking of the model’s loss given input data. We acknowledge that it would
be hard to derive any insight from an absolute itself, however, comparing
performance of a model across two (or more) semantic concepts could be eye
opening in terms of model’s expected areas of strength and weakness. In this
section, we first go over qualitative examples of how insights could be
extracted by performing such comparisons in Section. 4.2.1. We then provide
quantitative results in terms of how successful each of our three variations
are in terms of achieving such ranking in Section 4.2.2.

####  4.2.1 Qualitative

Here we show examples of how our approach at inference mode could be used to
generate explanations on the potential failure areas of the model, and
potentially uncover bias in the model’s performance across pairs of input
types. We visualize two examples of such insights in Figure 3. On both
examples we are visualizing $h(\tilde{g}(x_{1}))-h(\tilde{g}(x_{2}))$ where
$x_{1}$ and $x_{2}$ are two text inputs across which we would like to estimate
model’s performance disparity. Right shows the results of feeding the two
following textual inputs $x_{1}$=" face of an old male", and $x_{2}$=" face of
a young female". The figure shows the quadrants of probabilities that loss is
higher for $x_{1}$ compared to $x_{2}$. It can be observed that the model
thinks model $f$ would do better on inputs resembling "face of a young female"
compared to inputs resembling "face of an old male", as the loss is predicted
to be lower for "face of a young female" compared to "face of an old male".
Given that we trained ensembles of models, each individual model will have a
point-estimate prediction on the estimated-loss difference, leading to a
distribution of estimations from the ensemble. This allows us to measure the
likelihood of $h(\tilde{g}(x_{1}))-h(\tilde{g}(x_{2}))$ to be above (or below)
zero, indicating that the model is biased towards $x_{2}$ compared to $x_{1}$
(or vice-versa). The bars show the quadrants of the distribution and as it can
be observed the whole distribution is above zero indicating a fairly confident
estimation. On the other hand, the example on the left-hand side indicates
that for textual inputs $x_{1}$=" smiling long-haired person wearing glasses",
and $x_{2}$="non-smiling long-haired person wearing glasses" the prediction is
less certain. In the following section, we aim at quantifying this capability.

Figure 3: Two examples of pairs of textual inputs fed to the model for bias
estimation. Right is the results of feeding the pair $x_{1}$="face of an old
male", and $x_{2}$="face of a young female". And left is the results of
feeding $x_{1}$="smiling long-haired person wearing glasses", and
$x_{2}$="non-smiling long-haired person wearing glasses". It can be observed
that the model has a more confident prediction about the bias across
$x_{1}$="face of an old male" vs $x_{2}$="face of a young female". The results
on the right shows that the model confidently thinks that $f$ has a lower loss
for $x_{2}$, and therefore it’s performance is biased towards it.

####  4.2.2 Quantitative

Ranking datapoints based on their loss: In this experiment we pass each of the
images in celebA test dataset through $h(g(x))$ to get its estimated loss. We
then pass the same images through $f$ and measure their actual losses
$l(f(x),y)$. We then evaluate how well the two sets correlate. We evaluate
this for each type of $h$ (namely ranker, regressor, and classifier), and
report the results in Table 2. Given that we trained 5 models for each setup,
we report the results in ensemble mode (using mean and median aggregation of
the scores). We report the performance of the models individually in the
supplementary material. It can be observed that all models positively
correlate, However, when it comes to the ranking correlation metrics, the
ranking model strongly outperforms the other two variations, which is expected
as it is trained for this objective. Interestingly, the classifier seems to
consistently perform better than chance even though it’s objective does not
have any notion of absolute or relative value of loss for a datapoint. This
suggests the relativity of datapoints being in-distribution. In other words,
the datapoints that are more in-distribution (and thus probably over
represented), tend to have lower losses. This is consistent with the
observation in Section 4.1.2 where we observed the ranker and regressor having
great performance in terms of OOD detection.

| Pearson | Spearman | Kendall Tau
---|---|---|---
$h^{ens}_{class}$ | (0.01463, 0.086) | (0.13181, 3.850e-54) | (0.08958, 1.499e-53)
$h^{ens}_{reg}$ | (0.06851, 9.97e-16) | (0.27099, 3.75e-229) | (0.18626, 3.76e-225)
$h^{ens}_{rank}$ | (0.13288, 5.33e-55) | (0.56047, 0.0) | (0.39860, 0.0)
Table 2: Evaluating performance of the three variations of $h$ in terms of
ranking test data-points based on their estimated loss i.e. $h(g(x)$. The
correlation is evaluated between the predicted losses and the actual losses
$l^{f}(f(x),y)$ calculated based on model $f$ and ground-truth label $y$.

Ranking categories based on model’s predicted bias: Given that celebA[10]
dataset provides textual attributes for each of the data points, we measure
the ground-truth bias of model f across those textual attributes. To do so we
measure AUROC of the face recognition model for each binary of a semantic
attribute. As an example we measure mAUROC for all images with attribute
"Young=True", and for all images with attribute "Young=False". We then measure
the bias of the model across the two by measuring the Rawlsian [15] max min of
the two values. These values for each attribute are reported in the first
column of Table 5 under GT-True AUROC, GT-False AUROC and GT-RMM (Ground-truth
Rawlsian max-min). The higher the Rawlsian Max-Min the higher the disparity of
performance across the two categories. We then predict the model’s bias using
$h(\hat{g}(x_{text}))$. To do so we calculate the Rawlsian max-min between the
following two values: $h(\hat{g}(x_{text}="Young"))$,
$h(\hat{g}(x_{text}="NotYoung"))$. We simply infer on the attribute and itself
appended to "Not". We report this estimated Rawlsian Max-Min in column CLIP-
RMM of Table 5. Finally, we measure how our predicted bias CLIP-RMM is
correlated to the ground-truth bias GT-RMM. In Table 3 we report the
correlations in terms of Pearson, Spearman, and Kendall-tau ranking
correlation. It can be observed that for all three, the results positively
correlate with statistically significant p-values. Similar to the per
datapoint loss-based ranking task, the ranker variant (i.e. $h_{rank}$)
expectedly performs better in terms of predicting model’s bias across
categories. The regressor variant (i.e. $h_{reg}$) also seems to perform
meaningfully, while the classifier (i.e. $h_{class}$) variant, as expected
does not perform meaningfully, and seems to have a statistically non-
significant correlation with high p-values.

##  5 Limitations

Given that As all of our estimations depend on zero-shot capabilities of a
language model, namely CLIP[13], disproportionate reliability across different
areas of the embedding space could lead to certain semantic concepts being
better represented. This shortcoming may be to some extent being modeled by
modeling reliability across latent embedding spaces along the lines of [1],
however, the approach to integrate such notion is non-trivial. We conduct a
simple experiment to validate the existence of such effect. We report per-
attribute AUROC of clips zero shot accuracy on the test data and report it in
Table 5. This is simply measured by performing zero-shot accuracy on an image
and passing two textual inputs such as "Young" and "Not Young" and comparing
that to the ground-truth value of attribute "Young" for that datapoint. We
then measure how this zero-shot accuracy correlates with model robustness
across each attribute. We report the results in Table 4. As it can be observed
there is a statistically significant correlation. This suggests that on
attributes where CLIP prior is stronger, the face-recognition model happens to
be more robust. We hypothesize that this may be due to high correlation
between the biases of the two models, however, it would be hard to verify this
hypothesis.

- | Kendall-tau | Pearson | Spearman
---|---|---|---
$h_{reg}$ | (0.19, 3e-04) | (0.27, 2e-04) | (0.24, 5e-03)
$h_{class}$ | (0.064, 0.560) | (0.034, 0.834) | (0.072, 0.656)
$h_{rank}$ | (0.253, 0.021) | (0.324, 0.040) | (0.366, 0.019)
Table 3: Table of correlations between the predicted, and the actual biases. Correlation values are reported as tuples with the format: (correlation-coefficient, p-value).  Kendall-tau | Pearson | Spearman
---|---|---
(0.2743, 0.0126) | (0.291, 0.0679) | (0.3482, 0.0276)
Table 4: Bias Correlation. The correlation between the AUROC of a zero-shot clip based classification for a category, compared to the robustness of a face-recognition model across that category. attribute | GT-True AUROC | GT-False AUROC | GT-RMM | CLIP-RMM | 0-shot AUROC
---|---|---|---|---|---
Male | 0.9969 | 0.9912 | 0.0057 | 0.0040 | 0.8932
Young | 0.9932 | 0.9939 | 0.0006 | 0.0030 | 0.7501
Oval Face | 0.9960 | 0.9923 | 0.0037 | 0.0003 | 0.4775
Pale Skin | 0.9921 | 0.9934 | 0.0013 | 0.0010 | 0.7606
Big Lips | 0.9918 | 0.9941 | 0.0023 | 0.0016 | 0.5836
Big Nose | 0.9961 | 0.9926 | 0.0034 | 0.0003 | 0.5114
Pointy Nose | 0.9932 | 0.9934 | 0.0002 | 0.0009 | 0.6190
Narrow Eyes | 0.9912 | 0.9937 | 0.0025 | 0.0029 | 0.5345
5 o Clock Shadow | 0.9977 | 0.9929 | 0.0047 | 0.0058 | 0.3840
Arched Eyebrows | 0.9920 | 0.9939 | 0.0019 | 0.0020 | 0.6638
Attractive | 0.9926 | 0.9941 | 0.0014 | 0.0007 | 0.6934
Bags Under Eyes | 0.9951 | 0.9929 | 0.0021 | 0.0001 | 0.4412
Bald | 0.9957 | 0.9933 | 0.0024 | 0.0049 | 0.5656
Bangs | 0.9931 | 0.9934 | 0.0003 | 0.0008 | 0.8513
Black Hair | 0.9962 | 0.9923 | 0.0038 | 0.0011 | 0.8187
Blond Hair | 0.9898 | 0.9939 | 0.0041 | 0.0039 | 0.8847
Blurry | 0.9884 | 0.9936 | 0.0053 | 0.0013 | 0.6912
Brown Hair | 0.9926 | 0.9936 | 0.0010 | 0.0023 | 0.7756
Bushy Eyebrows | 0.9970 | 0.9928 | 0.0041 | 0.0024 | 0.7114
Chubby | 0.9965 | 0.9932 | 0.0032 | 0.0006 | 0.5286
Double Chin | 0.9967 | 0.9932 | 0.0034 | 0.0006 | 0.3657
Eyeglasses | 0.9936 | 0.9934 | 0.0002 | 0.0041 | 0.8945
Goatee | 0.9973 | 0.9932 | 0.0041 | 0.0022 | 0.5699
Gray Hair | 0.9972 | 0.9933 | 0.0039 | 0.0030 | 0.7494
Heavy Makeup | 0.9931 | 0.9936 | 0.0005 | 0.0015 | 0.4951
High Cheekbones | 0.9944 | 0.9923 | 0.0021 | 0.0004 | 0.5354
Mouth Slightly Open | 0.9939 | 0.9929 | 0.0009 | 0.0003 | 0.6058
Mustache | 0.9972 | 0.9932 | 0.0040 | 0.0021 | 0.5196
No Beard | 0.9927 | 0.9974 | 0.0046 | 0.0049 | 0.2825
Receding Hairline | 0.9955 | 0.9932 | 0.0023 | 0.0030 | 0.5250
Rosy Cheeks | 0.9941 | 0.9933 | 0.0007 | 3.8824 | 0.5282
Sideburns | 0.9976 | 0.9932 | 0.0043 | 0.0002 | 0.4113
Smiling | 0.9943 | 0.9924 | 0.0018 | 0.0010 | 0.9432
Straight Hair | 0.9938 | 0.9933 | 0.0004 | 0.0031 | 0.6136
Wavy Hair | 0.9924 | 0.9939 | 0.0015 | 0.0034 | 0.7287
Wearing Earrings | 0.9937 | 0.9933 | 0.0004 | 0.0009 | 0.7150
Wearing Hat | 0.9937 | 0.9934 | 0.0003 | 0.0012 | 0.8060
Wearing Lipstick | 0.9921 | 0.9948 | 0.0026 | 0.0013 | 0.8018
Wearing Necklace | 0.9908 | 0.9938 | 0.0029 | 0.0004 | 0.5787
Wearing Necktie | 0.9973 | 0.9931 | 0.0041 | 0.0030 | 0.3425
Table 5: per-attribute model performance AUROC, true performance-disparity
(GT-RMM) calculated using ground-truth CelebA[10] attributes, zero-shot
estimated disparity (CLIP-RMM), alongside accuracy of it’s zero-shot
classification (0-shot AUROC).

##  6 Discussion and Future Work

Given that this work is in the preliminary stages, there are many aspects of
it that could be improved and further explored. As an example, we believe that
even though the ranking model seems to lead to better quantitative estimations
(according to Table 3), the regressor model has other benefits that may have
not been manifested in this experiment. As an example, we believe that since
the regression model estimates absolute values of single datapoint loss
functions, we can use error bars on the single datapoint loss estimations in
conjunction with the pairwise ones to reason about the certainty of the
relative predictions. We aim to experiment with different definitions of
probability to take that into account. We also believe that not only more
extensive experimentation should be conducted across other types of modalities
and tasks, but also more extensive experimentation could be done within the
confines of the celebA dataset. Given that language models have shown to be
better at modeling more specific textual prompts, we suspect that our approach
could show its impact more, in more detailed group descriptions (the ones
shown in Figure 3). We aim to explore this in the future.

##  7 Conclusion

We explored a simple approach to leverage the strengths of language models in
zero-shot learning, to provide semantic insights into a trained model’s
patterns of strengths, weaknesses and biases. Our experiments suggest some
promising preliminary results in providing insights in terms a models
predicted familiarity and predicted performance on an arbitrary semantic
concept. We hope this encourages further investigation in this line of
research.

## References

  * [1] Shervin Ardeshir and Navid Azizan.  Embedding reliability: On the predictability of downstream performance.  In NeurIPS ML Safety Workshop.
  * [2] Shervin Ardeshir, Cristina Segalin, and Nathan Kallus.  Estimating structural disparities for face models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 10358–10367, June 2022.
  * [3] Hugo Berg, Siobhan Mackenzie Hall, Yash Bhalgat, Wonsuk Yang, Hannah Rose Kirk, Aleksandar Shtedritski, and Max Bain.  A prompt array keeps the bias away: Debiasing vision-language models with adversarial learning.  arXiv preprint arXiv:2203.11933, 2022.
  * [4] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  Advances in neural information processing systems, 33:1877–1901, 2020.
  * [5] Aakanksha Chowdhery, Sharan Narang, Jacob Devlin, Maarten Bosma, Gaurav Mishra, Adam Roberts, Paul Barham, Hyung Won Chung, Charles Sutton, Sebastian Gehrmann, et al.  Palm: Scaling language modeling with pathways.  arXiv preprint arXiv:2204.02311, 2022.
  * [6] Erik Härkönen, Aaron Hertzmann, Jaakko Lehtinen, and Sylvain Paris.  Ganspace: Discovering interpretable gan controls.  Advances in Neural Information Processing Systems, 33:9841–9850, 2020.
  * [7] Fangzhou Hong, Mingyuan Zhang, Liang Pan, Zhongang Cai, Lei Yang, and Ziwei Liu.  Avatarclip: Zero-shot text-driven generation and animation of 3d avatars.  arXiv preprint arXiv:2205.08535, 2022.
  * [8] Ameya Joshi, Amitangshu Mukherjee, Soumik Sarkar, and Chinmay Hegde.  Semantic adversarial attacks: Parametric transformations that fool deep classifiers.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 4773–4783, 2019.
  * [9] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [10] Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep learning face attributes in the wild.  In Proceedings of International Conference on Computer Vision (ICCV), December 2015.
  * [11] Qi Qi and Shervin Ardeshir.  Improving identity-robustness for face models.  arXiv preprint arXiv:2304.03838, 2023.
  * [12] Haonan Qiu, Chaowei Xiao, Lei Yang, Xinchen Yan, Honglak Lee, and Bo Li.  Semanticadv: Generating adversarial examples via attribute-conditioned image editing.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XIV 16, pages 19–37. Springer, 2020.
  * [13] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, et al.  Learning transferable visual models from natural language supervision.  In International conference on machine learning, pages 8748–8763. PMLR, 2021.
  * [14] Vikram V Ramaswamy, Sunnie SY Kim, and Olga Russakovsky.  Fair attribute classification through latent space de-biasing.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 9301–9310, 2021.
  * [15] John Rawls.  Justice as fairness: A restatement.  Harvard University Press, 2001.
  * [16] Aditya Sanghi, Hang Chu, Joseph G Lambourne, Ye Wang, Chin-Yi Cheng, Marco Fumero, and Kamal Rahimi Malekshan.  Clip-forge: Towards zero-shot text-to-shape generation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 18603–18613, 2022.
  * [17] Axel Sauer and Andreas Geiger.  Counterfactual generative networks.  arXiv preprint arXiv:2101.06046, 2021.
  * [18] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al.  Llama: Open and efficient foundation language models.  arXiv preprint arXiv:2302.13971, 2023.
  * [19] Can Wang, Menglei Chai, Mingming He, Dongdong Chen, and Jing Liao.  Clip-nerf: Text-and-image driven manipulation of neural radiance fields.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 3835–3844, 2022.