Quantification of the volume-fraction reduction of sheared fragile glass-forming liquids and its impact on rheology

This study determines the volume-fraction reduction of sheared fragile glass-forming liquids. We consider a group of hypothetical systems that consist of particles with anisotropic particle-size modulations yet have almost the same average particle configuration as actual systems under shear flow. Our molecular dynamics (MD) simulations demonstrate that one specific hypothetical system can reproduce the relaxation dynamics of an actual sheared system, and we identify the shear-flow effect on the particle size with anisotropic size-modulation of this specific system. Then, based on the determination of the particle size and the resultant volume fraction, we rationalize how slight decreases in the volume fraction significantly reduce the viscosity snf provide a nonlinear constitutive equation. Notably, the obtained rheological predictions, including the crossover shear rate from Newtonian to non-Newtonian behavior, can be expressed only in terms of experimental observables, showing a good agreement with the MD simulation results. Our perspective on the volume fraction under shear flow may provide new insights into the conventional concept of free-volume.


I. INTRODUCTION
Shear thinning is one of the most ubiquitous non-Newtonian flow behaviors in glassy materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].When an imposed shear rate γ is smaller than the crossover value γc , the shear viscosity η and the structural relaxation time τα under a given flow are the same as those at equilibrium ( γ = 0), η (eq) and τ (eq) α , respectively.In contrast, when γ > γc , η and τα decrease significantly as γ increases.This nonlinear flow response usually causes more complex phenomena, such as shear banding and fracture, drastically altering the mechanical properties.Thus, understanding and controlling the shear thinning properties are of particular importance in the design of processing of glassy materials; however, there is still no general consensus regarding the underlying mechanism of shear thinning.
For many fragile glass-formers under equilibrium conditions, it is known that the density n and the temperature T are not independent parameters, whereas their combined variable determines physical states, i.e., the so-called (power law) density scaling [37][38][39][41][42][43][44][45].The density scaling is naively interpreted as follows: increasing (decreasing) the temperature increases (decreases) thermal fluctuations and overlap between neighboring particles, which results in an effective decrease (increase) in the particle size and the volume fraction.In the sense of the density scaling, determining the effective volume fraction by setting the temperature and the (number or mass) density and then evaluating how the effective volume fraction changes by varying these parameters are fundamental.Now we may ask how such a physical picture for equilibrium liquids is modified for nonequilibrium liquids under an external flow field.In an applied flow field, the interaction potential remains unchanged, while the average structural configuration is anisotropically distorted according to the flow symmetry.This distortion subsequently modifies the overlap properties from those at equilibrium which can influence determining the effective particle size.In other words, the shear rate γ, in addition to the temperature T and the density n, may serve as an extra parameter to control the effective particle size.
The present study addresses these issues with the aid of molecular dynamics (MD) simulations of model fragile glass-formers.Attempts to quantitatively estimate the free volume or the effective volume fraction by directly analyzing actual sheared systems have not yet succeeded.Instead, in this study, we consider a group of hypothetical liquid systems with different particle sizes but with the same two-body pair correlation function as an actual sheared system.We find one such hypothetical system that can reproduce the relaxation dynamics of the actual sheared system.Then, we determine the effective particle size of the actual system by identifying it with the particle size of the specific hypothetical system.The detailed approach is as follows.(i) During the structural relaxation period τα , the particle configurations, on average, are preserved.Thus, in a flow field with shear rate γ the average structure undergoes shear deformation with a strain of γ = γ τα .The two-body pair correlation func-tion describes the extension (compression) of the particle configurations along the extension (compression) axis by γ/2 relative to the equilibrium state.However, we notice that the two-body pair correlation function can be interpreted differently.As shown in Sec.II, the particle configuration under shear flow can be reproduced by hypothetical anisotropic distortion operations applied to particles: we can prepare a group of hypothetical systems that share two-body pair correlation functions that are almost identical to that of the actual sheared system.(ii) The hypothetical distortion operation reproducing the steady structure under shear flow is described by two parameters: the degree of shear distortion and the particle overlap.There is arbitrariness in setting these parameters.Nevertheless, as shown in Sec.III B, a specific distortion operation can even reproduce the relaxation dynamics of the actual sheared system, whereby such arbitrariness can be removed.By identifying this specific operation with the actual shear-flow effect controlling the effective volume fraction, we provide a prescription for quantitatively determining the volume-fraction reduction in sheared fragile glass-forming liquids.The effective particle size along the extension axis of an externally applied flow field is the invariant reference, while in other directions it anisotropically decreases.(iii) Furthermore, in Sec.III C, by incorporating the volume fraction under shear flow, determined above, into the Doolittle equation of the structural relaxation time, we give a nonlinear constitutive equation.This constitutive equation can explain how a slight decrease (increase) in the volume fraction (free volume) significantly reduces the viscosity with high predictability.Although we previously derived similar rheological predictions based on more heuristic arguments in Refs.[32,33], here, we provide a detailed numerical investigation of the physical origin of the shearinduced reduction of the volume fraction or enhancement of the free volume.

II. THEORETICAL BACKGROUND
In this section, we first provide a detailed explanation of the theoretical background of our study, whose validity is examined in the following section using MD simulations.A binary particle system is employed in our MD simulations to prevent crystallization, whereas, in this section, a monodisperse system is assumed for the simplicity of the expressions.The generality of the discussion presented below is not lost under this assumption.
We assume that the constituent particles interact via the following inverse power law (IPL) potential: where r is the distance between two particles and ζ is an exponent that is sufficiently larger than 1.Hereafter, a system whose constituent particles interact via U (A) (r) is referred to as an A-system.

A. The effective volume fraction
In typical simulation studies of liquids, the thermal energy scale is comparable to ϵ.For such a case, when ζ ≫ 1, neighboring particles are strongly prevented from getting closer than a distance σ to each other.Therefore, in most literature, σ is conventionally set to the (soft) core or the particle size.In this setting, the particle volume and the volume fraction are simply given as and respectively.Here, N is the total number of particles, V is the system volume, and the particle number density is denoted as n(= N/V ).
Although the above setting of the volume fraction is simple, the nature of particle packings intrinsically depends on the temperature.That is, as stated in the introduction, increasing (decreasing) temperature increases (decreases) thermal fluctuations, and therefore, even at a fixed n, the overlap between neighboring particles is enhanced (reduced), resulting in an effective decrease (increase) in the particle size and the volume fraction.Such a competing effect between repulsive interparticle interactions and thermal fluctuations is incorporated in the effective volume fraction as follows.In the present system with the IPL potential Eq. ( 1), the physical state is characterized by a scaling variable [46][47][48], where the temperature T is measured in units of the Boltzmann constant.A similar scaling variable can be defined in binary mixtures with the additive IPL potentials [57] and even in more general fragile systems that can be mapped onto those of effective IPL systems [40].In Eq. ( 4), the factor (ϵ/T ) 3/ζ characterizes the degree of the particle overlap, and therefore, we interpret ϕ as the effective volume fraction at equilibrium.However, whether the somewhat simplified volume fraction φ or the effective volume fraction ϕ is used does not matter in a practical sense.When n is varied at a fixed T , φ and ϕ are essentially the same.On the contrary, when T is varied at a fixed n, the physical states are usually characterized not by ϕ but by T .Under an applied flow field, the particle structures and the resultant particle-overlaps are anisotropically modulated according to the given flow symmetry.We would like to naively ask whether such a flow-induced anisotropy affects the above scaling.More specifically, we want to know whether the effective volume fraction is changed by varying the shear rate γ even at fixed n and T .If so, we may further raise a question about how ϕ is changed and its impact on the relaxation dynamics.In the following, we provide a detailed explanation of our perspective on these issues and the physical background.

B. Particle configurations under stationary shear flow
First, let us consider the average particle configurations under shear flow.In this study, a simple shear flow with the following mean velocity profile is assumed, where the x axis is along the direction of the mean flow, the y axis is along the mean velocity gradient, and x is the unit vector along the x axis.For a fixed T condition, φ and ϕ are essentially the same at γ = 0. Therefore, we hereafter reset the reference size and the volume fraction at γ = 0 to be σ and ϕ 0 = nπσ 3 /6, respectively.Then, we examine how they effectively vary as γ changes.Throughout this study, the average configurations are considered in terms of the two-body pair correlation function.Under the shear flow of Eq. ( 5), according to the flow symmetry, g s (r; ϕ 0 ) is generally expressed as [50,51] s (r; ϕ 0 ) + xŷg (1)  s (r; where x = x/r = sin ψ cos θ and ŷ = y/r = sin ψ sin θ, with ψ and θ representing the polar and azimuthal angles, respectively.Here, g s (r; ϕ 0 ) represents the isotropic part and xŷg (1) s (r; ϕ 0 ) is the leading order deviation from g (0) s (r; ϕ 0 ).It is known that g (1) s (r; ϕ 0 ) is responsible for the non-zero average shear stress, Σ xy , which is given as [50,51] Σ xy = 1 2 n 2 drx 2 ŷ2 r dU (A) dr g (1)  s (r; ϕ 0 ), The deviatoric part g s (r; ϕ 0 ) is approximately given by [52][53][54] g (1)  s (r; ϕ where c g is a numerical constant of order unity.As shown in the next section, this approximate form of g s with an appropriate value of c g reproduces the simulation results well.Equation ( 8) can be understood as a consequence of balancing advection and relaxation in the steady state [52] as follows By taking the leading order term of γ, we obtain Eq. ( 8).Substituting Eq. ( 8) into Eq.(6) yields where λ = c g γ τα (λ ≪ 1).Equation (10) indicates that the average particle configuration under the shear flow is approximately given by distorting the reference configuration g (0) s (r; ϕ 0 ), i.e., by elongating and compressing by (1 + λ/2) and (1 − λ/2) along xŷ = 1/2 and xŷ = −1/2, respectively.
C. Reproduction of the sheared configuration by distorting an arbitrary reference Next, we note that Eq. ( 10) can be expressed differently.Namely, within the present leading order approximation in λ(≪ 1), we can formally rewrite Eq. (10) as where we set b ∈ [−1/2, 1/2].Note that Eq. ( 11) includes Eq. ( 10) when b = 0.This formal re-expression of g s (r; ϕ 0 ) indicates that an identical sheared particle configuration can be obtained by applying a hypothetical distortion to an isotropic reference configuration described by g s (r; ϕ 0 )].A more detailed explanation of this reexpression is presented in Appendix A.
By denoting the operation of Eq. ( 12) as D λ,b , we may express the present hypothetical operation as These situations are schematically illustrated in Fig. 1.Note that, as observers, we generally consider the operation D λ,0 to correspond to the actual occurring distortion.

D. Reproduction of the sheared configuration by distorting constituent particles
In Eq. ( 11), we consider the system to be anisotropically distorted while the constituent particles remain undistorted.However, we may interpret Eq. ( 11) differently: particles are anisotropically distorted, while the system remains undistorted.That is, with σ being the reference size, the particle size is anisotropically modulated as which is assumed to be set by the following interaction potential (1 ) The particle volume and the volume fraction are given as and respectively.Hereafter, a system where constituent particles interact via U (B) (r) is referred to as a B-system.By denoting the operation of Eq. ( 14) as Dλ,b , we may express the present hypothetical operation, Eq. ( 14), as The situations for three different values of b are schematically displayed in Fig. 2. In the next section, we demonstrate that the two-body pair correlation function of the B-system obtained by MD simulations for various values of b and λ nearly exactly reproduces the stationary sheared structure.
Assuming that the virial theorem holds in the present B-system, the compressibility factor p/nT is expressed as A more detailed derivation is presented in Appendix B. In Eq. ( 19), p formally represents the equilibrium pressure of the A-system at a volume fraction of (1 − 3λb)ϕ 0 and a number density of n.The pressure is smaller for larger b at a fixed number density n.Furthermore, the pressure component along arbitrary directions does not depend on the direction and, thus, is isotropic.This resultant isotropy reflects the absence of distortion of the system itself.Based on the view presented in this subsection, Eq. ( 11) implies that the configurations of differently modulated particles [by Eq. ( 14)] share almost identical two-body correlation functions [55].This can be naively understood as follows.At a fixed system size, enhancing the particle size (decreasing b) reduces fluctuations but increases the pressure.These two competing pressure and fluctuation effects balance, and the resultant particle configuration remains almost unchanged.

E. Tuning the overlap by tuning the potential
In the B-system, the equation of motion of the i-th particle's dynamics may be simply given by, where R i is the position of the i-th particle and R ij = R i − R j .Equation ( 20) is formally rewritten as where Equations ( 20) and ( 21) correspond to the views described in Secs.II D and C, respectively.That is, regarding the equation of motion, the difference in these views corresponds to the difference in whether F i is regarded as the intrinsic or extrinsic force.Note that neither Eqs. ( 20) nor ( 21) describe the actual particle dynamics under shear flow.Still, these dynamics may reproduce the stationary structural configurations of the sheared system, as shown in the next section.In Secs.II B-D, we discussed that the "average" shear-flow effects on the particle configurations can be reproduced by applying hypothetical distortions to the system or the particles even without shear flow.However, in reproducing the steady structure, g s (r; ϕ 0 ) ∼ = g (0) there is an arbitrariness in choosing b.In the next section, we use MD simulations to demonstrate that only a specific operation with b = 1/2 can appropriately reproduce the actual relaxation dynamics, removing such arbitrariness.

III. NUMERICAL RESULTS
In this section, following the arguments presented in the previous section, we perform MD simulations to demonstrate that the structure and the dynamics of the sheared system can be mapped onto a system in which the constituent particles interact via the anisotropically modulated potential U (B) with a specific value of b(= 1/2).
A. Simulations of the sheared A-system with isotropic potentials In our simulations, we employ a binary mixture of large (L) and small (S) particles interacting via the (soft core) IPL potentials given by [56][57][58] where µ, ν = L, S and r is the distance between two particles.σ µν = (σ µ +σ ν )/2, where σ µ is conveniently set to the size of the µ species particle in the reference state; that is, similar to the setting of the particle size in the monodisperse case (Sec.IIA), we also set the reference particle size when γ = 0 to be σ µ (µ = L, S).Under this setting, the reference particle volume and the volume fraction are given as and respectively.The mass and size ratios are m L /m S = 2 and σ L /σ S = 1.2, respectively.The units for length and time are σ S and (m S σ 2 S /ϵ) 1/2 , respectively.The total number of particles is N = N L + N S = 8000 and N L /N S = 1.The temperature T is measured in units of ϵ/k B .The fixed particle number density and the linear dimension of the system are N/V = 0.8/σ 3 S and L = 21.54,respectively.In this simulation, under simple shear flow, Eq. ( 5), the equations of motion are solved using Lee-Edwards periodic boundary conditions with a Gaussian thermostat [59].
In Fig. 3, we show the γ − η curves for the present model, which exhibit shear-thinning behavior.Crossovers from Newtonian to non-Newtonian flow behavior at γτ (eq) > 1 with τ (eq) being the equilibrium relaxation time are observed in many soft matter systems: γτ (eq) > 1 indicates the dominance of advective effects over equilibrium structural relaxation mechanisms (the so-called constitutive instability) in flows.Therefore, similar crossovers might be expected to occur in glass-forming liquids.However, as shown in Fig. 3, by focusing on the average degree of the shear "distortion" at the crossover, we find that shear thinning starts when γ is several orders of magnitude smaller than 1/τ (eq) α , which indicates quite a small average structural distortion (λ ≪ 10 −2 ).This is also the case in most experiments [12][13][14] and simulations [6-8, 10, 11, 31-34, 36] of supercooled liquids, where the onset of shear thinning occurs at approximately γτ (eq) α ∼ 10 −2 ∼ 10 −3 .This large time-scale separation may exclude the possibility of the usual constitutive instability [11,[31][32][33], and, therefore, is an important characteristic of rheological features observed near the crossover from Newtonian to non-Newtonian flow behaviors.In the standard simulation of supercooled liquids, the Lindemann length is typically approximately 0.1 times the particle size; therefore, the amplitude of strain fluctuations at the particle scale is much larger than the average shear distortion γ τα .In such a situation, the shear-flow effects should be much weaker than the thermal fluctuation effects and are regarded as a small perturbation to the structures.However, as shown in Sec.III C, they should have a strong impact on the dynamics., indicating γcτ α (eq) ≪ 1.In the inset, we plot 1/τ (eq) α , γc, and the theoretically predicted crossover shear rate ϕ0(∂τ (eq) α /∂ϕ0) −1 (Eq.( 43) derived below in Sec.IIIC) against 1/T .We find that γc quantitatively corresponds with ϕ0(∂τ Similar to Eq. ( 6) for the monodisperse case, under the shear flow of Eq. ( 5), the pair correlation functions g s,µν (r; ϕ 0 ) are expressed as [51] g s,µν (r; ϕ 0 ) = g (0)  s,µν (r; ϕ 0 ) + xŷg (1)  s,µν (r; where g s,µν (r; ϕ 0 ) represents the isotropic part and xŷg (1) s,µν (r; ϕ 0 ) is the leading order deviation from g (0) s,µν (r; ϕ 0 ).Similar to Eqs. ( 8) and (10) presented in Sec.II B, the deviatoric part g s,µν (r) is approximately given as [52][53][54] g (1)  s,µν (r; and we obtain Here, c g is a numerical constant of order unity and λ = c g γ τα (λ ≪ 1).As shown in Fig. 4(b), for the present model system, this approximate form of g s,LL (r) (b) for various shear rates at T = 0.285 (left) and 0.267 (right).To leading order in γ, the pair correlation function gs,LL(r; ϕ0) is expressed as gs,LL(r; ϕ0)=g s,LL (r; ϕ0).In (a), g s,LL (r; ϕ0) collapses onto a single curve.As shown in (b), for particle pairs in the first shell, g s,LL (r) is well approximated by g s,LL (r; ϕ0), where λ = cg γ τα(≪ 1), and we set cg = 0.65 and 0.45 for T = 0.285 and 0.267, respectively.Almost the same results are obtained for gs,SS(r; ϕ0) and gs,SL(r; ϕ0) as those for gs,LL(r; ϕ0) using the same value of cg.
being a constant of order unity [60], nearly reproduces the simulation results.

B. Simulations of the unsheared B-system with anisotropic potentials
As discussed in Secs.II C and D, there are two ways to interpret Eq. ( 29).One is that the system is anisotropically distorted while the particles remain undistorted.In this subsection, we show that the structure and the relaxation dynamics of the system obtained by the hypothetical operation Dλ,1/2 can reproduce those of the actual sheared system.For this purpose, let us consider the following interaction potentials for binary mixtures with the effective particle size of the µ-species being anisotropically modulated as Here, λ and b are the parameters controlling the degree of distortion and the size of the particles, respectively, as in the monodisperse case discussed in Sec.II D. The particle volume and the volume fraction are given as and respectively.The other settings are the same as those of the A-system presented in Sec.III A. We simulate the present model, where constituent particles interact via U (B) µν (r ij ) without shear flow, using velocity Verlet algorithms in the NVE ensemble [59].
Before proceeding, we note the following.Since the offdiagonal components of the stress tensor are not symmetric due to the asymmetric form of U (B) , the net torque is not exactly zero.However, the particle configurations are distorted so that the resultant local torques are sufficiently suppressed.Therefore, there are no strange rotational motions in the B-system.Furthermore, in both the sheared A-and unsheared B-systems, as long as λ is sufficiently small, the dynamics are almost isotropic.For the sheared A-system, the deviatoric particle motions, in which the contribution from convective transport by the average shear flow is subtracted, show minimal marked anisotropy at the two-body correlator level [61].However, although some anisotropies emerge from longer-time behaviors and are captured in the dynamic heterogeneity or shear bandings [11], their roles in the rheological properties remain poorly understood.
As discussed for the monodisperse case in Sec.II, Fig. 5 shows that the particle configurations of the B-system for different b(∈ [−1/2, 1/2]) have almost identical twobody correlators that are approximately described as For the sheared A-system, as long as λ = c g γ τα ≪ 1, g s,µν (r; ϕ 0 ) ∼ = g (eq) µν (r; ϕ 0 ).Therefore, from Eqs. ( 28) and ( 35), we deduce from which we may conclude that the unsheared Bsystem can approximately reproduce the average particle configurations of the sheared A-system.

Dynamics
In Fig. 6, the structural relaxation time of the Bsystem, which hereafter is denoted as τα (λ; b, ϕ 0 , T ), is plotted against λ for several values of b.In the present study, the structural relaxation time is defined as the relaxation time of the shear-stress autocorrelation function.For more details, please refer to Appendix C. At the same λ, despite almost the same particle configurations for different b (shown in Fig. 5), the behaviors of the structural relaxation times are quite different.When b < 0, τα (λ; b, ϕ 0 , T ) increases with an increase in λ, reflecting the increase of the volume fraction (and the resultant pressure).However, when b > 0, the opposite result occurs.Remarkably, when b = 0, the relaxation time remains unchanged, which suggests that a small anisotropy (λ ≪ 1) without any volume changes does not affect the structural relaxation.In Fig. 6, we also show the relaxation time of the sheared A-system τα ( γ; ϕ 0 , T ) against the degree of the average shear distortion λ = c g γ τα .We find that τα ( γ; ϕ 0 , T ) and τα (λ; b = 1/2, ϕ 0 , T ) nearly coincide with each other at various temperatures.
as schematically shown in Fig. 7, the average interparticle potential becomes minimal along the extension axis.Then, setting this direction to be the reference, in other directions, the potential energy is "lifted" up due to the shear flow, making extra particle overlaps in addition to overlaps due to thermal fluctuations.By subtracting the extra overlap regions, Eq. ( 37) with b = 1/2 describes the effective particle size of the actual sheared system, and the effective volume fraction is given by ϕ 0 (1−dλ/2) at a fixed T .Here, λ = c g γ τα is determined by the distortion of the two-body pair correlation function.
This shift of the potential energy reference under the shear flow does not alter the observables, such as the average energy, pressure, and shear stress.Although our simulations certainly support the present speculation regarding the reduction of the volume fraction induced by the shear flow, the detailed investigation based on first principles are required to provide further evidence.
3. Doolittle equation: Mapping onto the equilibrium system with the reduced volume fraction Note again that for λ ≪ 1, anisotropy is hardly noticeable in the dynamics at the two-body correlator level [61].We expect that the shear-flow effect is incorporated only through the reduction of the volume fraction and, thus, that the dynamics of the sheared system can be mapped onto the dynamics of the equilibrium system.
In Fig. 6, we also plot the Doolittle equation [62,63] with the reduced volume fraction Here, the parameters of Γ, ϕ c , and τ D 0 generally depend on the temperature and are separately determined at equilibrium using MD simulations.As shown in Fig. 8(a), the Doolittle equation approximates the volumefraction dependence of the structural relaxation time at equilibrium well.In Fig. 6, we find that τ D α (ϕ s ; T ) with Eq. ( 39) quantitatively reproduces the effect of τα ( γ; ϕ 0 , T ) and τα (λ; b = 1/2, ϕ 0 , T ) with the same λ, further supporting our perspective.
Because of a very steep volume-fraction dependence of the relaxation time as described in Eq. (38) for supercooled states, even an infinitesimal reduction of the volume fraction causes significant acceleration of the dynamics.

C. Nonlinear constitutive equations
As discussed in Sec.IIB, based on the results shown in Fig. 6, we suppose that the relaxation time under shear flow τα is mapped onto the equilibrium τ α as τα ( γ; ϕ 0 , T ) = τ (eq) α (λ; ϕ 0 , T ) = τ (eq) where S )/2 from Eq. ( 25), and the shearflow effect is taken into account through the reduced volume fraction ϕ s .This equation is essentially nonlinear in τα .If we know the functional form of τ (eq) α (for example, the Doolittle equation), we can solve Eq. ( 40) in terms of τα .Equation ( 40) can be regarded as the nonlinear constitutive equation and describes the rheological curves of the present model well.Since the volume fraction dependence of the shear modulus G is much weaker than that of τ α , the viscosity may be taken to be η for (ϕ 0 − ϕ s )/ϕ 0 ≪ 1.Note that the T -dependence of G is also much weaker than that of τ (eq) α .For γ τα ≪ 1, by expanding Eq. ( 40) in γ τα , we obtain τα ( γ; ϕ 0 , T ) ∼ = τ (eq) where a = 3c g /2 ( ∼ = 1 in the present system).Therefore, the crossover shear rate from Newtonian to non-Newtonian behavior γc is given by γc This crossover shear rate can be much smaller than 1/τ (eq) α ( γc τ (eq) α ≪ 1) near the glass transition point, indicating that the usual constitutive instability does not trigger the onset of the shear-thinning.The reference volume fraction ϕ 0 and the number density n are linearly related to each other (ϕ 0 ∝ n), and therefore, Eq. ( 43) is rewritten as Equation ( 43) is expressed only in terms of experimental observables, and should thus be useful in the process design of glassy materials.Similar predictions for rheological behaviors were obtained in Ref. [32] by more heuristic arguments; here, we rationalize the possible mechanism of the shear-induced reduction of the volume fraction.The constitutive equation (40) can be quantitatively approximated by the Doolittle equation ( 38) with the reduced volume fraction ϕ s in Eq. ( 39) or equivalently with the enhanced "free volume" ∝ (ϕ c − ϕ s ).Although our perspective on the enhancement of the free volume under shear flow is not exactly the same as the conventional perspective [18], we expect that the present study will provide new insights into the physical substance of the free volume.

IV. CONCLUDING REMARKS
In this study, we have discussed how the particle size and volume fraction of fragile liquids under shear flow are determined.Based on this determination, we derive a nonlinear rheological constitutive equation, which quantitatively describes the shear thinning behavior of fragile supercooled liquids.
In a shear flow with a shear rate γ, the particle structures relax with a time scale of τα , resulting in an average distortion on the order of γ τα .The extent to which neighboring particles can overlap is involved in determining the volume fraction; such particle overlap should be controlled by the degree of shear distortion (∝ γ τα ) in addition to the strength of thermal fluctuations (∝ T ).Under the simple shear flow of Eq. ( 5), dilution and densification occur along the extension and compression axes, respectively, with the same magnitude, resulting in the absence of a system volume change.In this situation, the overall number density n, which is uniquely determined, is invariant.However, unlike the number density, the volume fraction ϕ decreases as γ increases.In a system where constituent particles interact via simple short-range repulsive potentials, the average interparticle interactions become minimal along the extension axis (xŷ = 1/2).By setting the reference direction measuring the potential energy to this extension axis, the potential energy can be considered to be "lifted" up in other directions.By identifying this lift effect with the shear flow effect inducing extra particle overlaps (in addition to the overlaps due to thermal fluctuations), we obtain Eq. (37).Note that if the reference direction is set to xŷ = 0, which we as observers usually consider, the particle volume does not vary with the degrees of the shear distortion λ.However, this may not be true for particles.Please also refer to the discussion presented at the end of Sec.IIB3.
Because significant anisotropy is hardly observed in the dynamics at the two-body correlator level [61], we may consider the shear-flow effect to be incorporated primarily by a slight reduction of the volume fraction.Therefore, assuming that the sheared dynamics can be mapped onto the equilibrium dynamics, we obtained the nonlinear constitutive equation, i.e., Eq. ( 40).As clearly shown in Fig. 6, this constitutive equation quantitatively describes the shear-induced acceleration in the relaxation dynamics of fragile supercooled liquids.
Finally, we note the following: (i) Our preliminary simulations for another standard fragile model liquid, the Kob-Andersen model [65], also reproduce almost the same results as those obtained in the present paper.Although the model employed in this study assumes pairwise IPL potentials, the interactions of other systems are generally more complicated.However, it has been established that the dynamics of a wide class of fragile glass-formers can be reproduced by the corresponding IPL systems [40,45], for which the effective particle size and volume fraction are simply defined.
(ii) In Refs.[29,30], shear thinning in glassy liquids was discussed in the context of the shear-induced density inhomogeneity, inspired by the shear-induced phase separation established in polymeric systems [66][67][68].The critical shear rate for the onset of the inhomogeneous flow is given as γcr = (∂η (eq) /∂p) −1 T , where p is the pressure and is identified with the shear rate describing the crossover from Newtonian to non-Newtonian behavior.This prediction seems to agree with the experimental results for the supercooled melts of Zr-based bulk metallic glass-formers [14].However, as was pointed out in Ref. [14], the flow is still homogeneous when γ ∼ γcr , and inhomogeneous flow occurs when γ ≫ γcr , which contradicts the results of Refs.[29,30].Note that (∂η (eq) /∂p) T ∼ = G(∂τ (eq) α /∂p) T = G(∂n/∂p)(∂τ (eq) α /∂n) = GK T n(∂τ (eq) α /∂n) [32].Here, K T is the isothermal compressibility and GK T ∼ = 0.3, which was estimated from experimental results [14,64].Therefore, γc and γcr are comparable to each other since γc = 0.3 γcr , and the agreement between the predicted γcr and the experimental results may instead indicate the validity of the present mechanism of the shear-induced reduction of the volume fraction, which is not related to the density inhomogeneity.We will discuss which mechanism is selected under actual experimental situations in more detail elsewhere.
(iii) In this study, the considered systems are in (supercooled) liquid states, where thermal fluctuations exert important effects.However, thermal effects are irrelevant in amorphous solid states.In amorphous states, close links between the shear distortion of microscopic configurations and nonlinear rheological properties have been intensively studied in Refs.[69,70].At this stage, it is unclear how our approach for liquid states can be related to amorphous rheology.
(iv) In a hard core system, particle overlaps never occur.However, in this case, what becomes anisotropic is the collision frequency: namely, the collision frequency is enhanced along the compression axis while it is reduced along the extension axis.If we set the extension axis as the reference direction for measuring the collision frequencies, such anisotropies in collisions may be analogous to anisotropic particle overlaps in a soft core system.We will examine points (i)-(iv) in future work.This work was supported by KAKENHI (Grant No. 26103507, No. 25000002, and No. 20508139) and the JSPS Core-to-Core Program "International research network for nonequilibrium dynamics of soft matter".
In the last line r ′ is replaced by σs.Equation (B4) corresponds to the equilibrium pressure of the A-system with the reduced particle size (1 − λb)σ at a volume fraction of (1 − 3λb)ϕ 0 and a number density of n.

FIG. 2 :
FIG. 2: (Color online) A schematic showing the reproduction of the sheared particle configuration by anisotropically distorting constituent particles.(a) Anisotropically distorted constituent particles { Dλ,b : σ → σ[1 + λ(xŷ − b)]}.(b) Particle configurations obtained by Dλ,b .When b = 1/2, Dλ,1/2 anisotropically reduces the particle size, while when b = −1/2, Dλ,−1/2 anisotropically enhances the particle size.When b = 0, Dλ,0 represents particle distortion without changing the particle volume.Changing the particle size, which is controlled by changing b, varies the fluctuation and pressure effects.These competing fluctuation and pressure effects balance at a fixed number density n, resulting in an identical two-body pair correlation function even for different values of b.

FIG. 7 :
FIG.7:(Color online) (a) A schematic of the shift of the potential energy reference.For the sheared nonequilibrium system, the average interparticle potential energy E(θ) that a particle experiences is minimal (E0) along the extension axis (θ = π/4).Here, Eeq is the average value at equilibrium.In the nonequilibrium sheared system, our simulations indicate that E(θ) − E0 is considered to be due to extra overlaps due to the shear flow.(b) A schematic of the anisotropically modulated effective particle size.The dashed line represents the size at equilibrium.

FIG. 8 :
FIG. 8: (Color online) (a) The structural relaxation time of the equilibrium A-system τ (eq) α (ϕ0, T ) for several T , which can be fitted to the Doolittle equation (in terms of ϕ0) τ D α (ϕ0, T ) = τ D 0 exp[Γϕ0/(ϕc −ϕ0)], represented as the dashed curves.(b) A schematic for the acceleration of structural relaxation caused by the shear-induced reduction of the volume fraction.As discussed in the main text, due to a small anisotropy at the two-body correlator level, we expect that the dynamics of the sheared system can be mapped onto the equilibrium dynamics by incorporating the shear-flow effect only by reducing the volume fraction: ϕ0 → ϕs = (1 − 3cg γ τα/2)ϕ0.Close to the glass transition temperature, the volume-fraction dependence of τα becomes much steeper; thus, even a very small decrease in the volume fraction significantly accelerates the relaxation dynamics.