A two-fluid model for the formation of clusters close to a continuous or almost continuous transition

Experiments have shown that spatial heterogeneities can arise when the glass transition in polymers as well as in a number of low molecular weight compounds is approached by lowering the temperature. This formation of “clusters” has been detected predominantly by small angle light scattering and ultrasmall angle x-ray scattering from the central peak on length scales up to about 200 nm and by mechanical measurements including, in particular, piezorheometry for length scales up to several microns. Here we use a macroscopic two-fluid model to study the formation of clusters observed by the various experimental techniques. As additional macroscopic variables, when compared to simple fluids, we use a transient strain field to incorporate transient positional order, along with the velocity difference and a relaxing concentration field for the two subsystems. We show that an external homogeneous shear, as it is applied in piezorheometry, can lead to the onset of spatial pattern formation. To address the issue of additional spectral weight under the central peak we investigate the coupling to all macroscopic variables. We find that there are additional static as well as dissipative contributions from both, transient positional order, as well as from concentration variations due to cluster formation, and additional reversible couplings from the velocity difference. We also briefly discuss the influence of transient orientational order. Finally, we point out that our description is more general, and could be applied above continuous or almost continuous transitions


Introduction
The appearance of heterogeneities on many length and time scales as the glass transition is approached from higher temperatures has been of interest for a number for years. In parallel it has also become clear that the formation of spatio-temporal variations is also playing a role in connection with liquid-liquid transitions in onecomponent systems -compare, for example, a recent paper by Tanaka's group (Takae and Tanaka 2020). Experiments revealing spatio-temporal heterogeneities have come in several groups. One class are light scattering experiments for the central peak (Forster 1975;Berne and Pecora 1976) -also called the Rayleigh line - (Fischer 1993;Kanaya et al. 1994;Patkowski et al. 2001a;Patkowski et al. 2001b;Fischer et al. 2002) including also photon correlation spectroscopy (Walkenhorst et al. 1998;Fischer et al. 2002). Closely related are ultrasmall angle x-ray scattering experiments on the central peak (Patkowski et al. 2000;Fischer et al. 2002). In all these experiments one observes spatial heterogeneities (clusters) on length scales from ∼ 10nm → 200nm. In addition, one finds additional spectral weight under the central peak compared to what one would expect from a simple fluid (Forster 1975;Berne and Pecora 1976). We would also like to refer to the light scattering experiments on the ultraslow relaxation of hydrogen-bonded dynamic clusters in glassforming aqueous glucose solutions (Sidebottom 2007). In passing we mention that large-scale spatial heterogeneities have also been observed by light scattering in boron oxide glasses (Bokov 2016). A second fairly large class of observations of clusters of large spatial extent up to about ∼ 15μm) are mechanical experiments, in particular using piezorheometry (Collin and Martinoty 2003;Pozo et al. 2009). While the former have been done on linear polymer melts of rather low molecular weight, the latter were done in the isotropic phase of a liquid crystalline sidechain polymer. This study (Pozo et al. 2009) also complements earlier studies in the isotropic phase of sidechain liquid crystalline polymers above the nematic-isotropic phase transition (Gallani et al. 1994;Martinoty et al. 1999). In all cases when piezorheometry was used, a gel-like behavior is observed below a critical sample thickness of about 20μm indicating the maximum cluster size to be of that order of magnitude. Since low frequency ac shear of small strain amplitudes ∼ 10 −4 ) is applied, it was concluded (Collin and Martinoty 2003) that density fluctuations of large spatial extent and long lifetimes are present indicating transient positional order. Other mechanical experiments in the same spirit of detecting cluster formation have been done using low frequency anelastic spectroscopy (Wu and Zhu 2007;Shang et al. 2009). Another technique that has been used to investigate cluster formation is dielectric spectroscopy (Fischer 1993;Fischer et al. 2002;Kaminski et al. 2010). There is also evidence for the formation of dynamic heterogeneities from multi-dimensional NMR experiments indicating length scales of interest for cluster formation of about 3nm (Tracht et al. 1998). To sum up there is a large body of experimental evidence indicating cluster formation on many length and time scales, in particular as the glass transition is approached from above, mainly in polymers but also in low molecular weight materials. More recently it has also become clear in the field of biological physics that spatial heterogeneities in crowding cellular fluids on length scales of about 100nm play an important role (Stiehl and Weiss 2016;Collins et al. 2019;Donth and Weiss 2019) thus forming a bridge to the nonliving systems discussed above.
In the field of complex fluids two-fluid behavior is quite frequent in systems such as polymers and microemulsions. It arises, for example, for multiphase flows (Drew and Passman 1998), viscoelasticity of concentrated emulsions (Hebraud et al. 2000) and flow of colloidal suspensions (Lhuillier 2001). As examples for flows we mention the rheological behavior of polymer solutions and blends (Onuki 1989;Saito et al. 2001) and of polymer migration and phase separation under flow (Sun et al. 1999;Araki and Tanaka 2001).
We also emphasize that our description is more general, and could be applied above continuous or almost continuous transitions for which spatial heterogeneities (clusters) can grow to a size for which macroscopic dynamics enters the picture, that is above phase transitions, which are of second or weakly first order. We do not cover critical fluctuations in the renormalization group sense and their nonlinear properties.
The macroscopic description of two-fluid systems of any type brings along two additional macroscopic variables, namely the velocity difference between the velocities of the two subsystems and the concentration of one subsystem, φ, in addition to the total density ρ, which is always conserved, and the mean velocity, v i . Whether the concentration density φ is conserved or a relaxing quantity depends on the system of interest. For example, for mixtures of immiscible liquids or for microemulsions it will be a conserved quantity, while for one-component systems, which have a tendency towards cluster formation such as smectic clusters in a nematic liquid crystal (Brand and Pleiner 2021a) or singlecomponent substances showing a liquid-liquid transition (Takae and Tanaka 2020), φ will be a relaxing quantity. For the pioneering papers for the macroscopic dynamics of twofluid systems discussing fluid mixtures, two-fluid nematics and elastomer-fluid mixtures we refer to the publications by Pleiner and Harden (Pleiner and Harden 2003;. This work has since been generalized to several other systems including two-fluid effects in magnetorheological fluids (MRFs) , to immiscible compound materials in solids and gels (Pleiner et al. 2021) as well as to a nematic liquid crystal with smectic clusters to address the issue of the breakdown of flow alignment (Brand and Pleiner 2021a). This approach has also been used to describe bio-inspired materials (Pleiner et al. 2013;. There is one important aspect for the macroscopic dynamics of normal fluid two-fluid systems to be kept in mind: the velocity difference is always a macroscopic variable, which relaxes on a long, but finite time scale. This can be traced back to the fact that in normal fluid systems there is only the barycentric velocity, which is associated with a truly conserved variable, the density of linear momentum. Therefore, all the normal fluid systems discussed so far must be distinguished from superfluids for which one has two velocities associated with truly hydrodynamic variables. The additional truly hydrodynamic superfluid velocity is associated with a spontaneously broken continuous symmetry, namely gauge invariance (Khalatnikov 1965;Hohenberg and Martin 1965;Forster 1975). Hydrodynamic equations for superfluids have been derived in a number of papers; compare, for example, for superfluid 4 He (Khalatnikov 1965;Hohenberg and Martin 1965;Forster 1975), for the various superfluid phases of 3 He in the bulk and in aerogels (Graham 1974;Graham and Pleiner 1975;Liu 1976;Brand et al. 1979;Liu 1979;Brand and Pleiner 1982; and for superfluid neutron star matter (Brand and Pleiner 1981b).
To derive macroscopic equations for a fluid with clusters we use the systematic approach of linear irreversible thermodynamics (de Groot and Mazur 1962;Martin et al. 1972;Pleiner and Brand 1996). We discuss this subject in detail in the "Two-fluid model for a fluid with clusters" section. First, we clarify the macroscopic variables for such a fluid and then we present the basic hydrodynamic equations including the macroscopic variables for a twofluid description of a fluid with clusters: the relative velocity w i and the concentration of clusters. Since it became clear that also transient positional order and transient orientational order play an important role in the interpretation of the experiments, these macroscopic variables are incorporated as well into the basic equations.
In the "Response to an external homogeneous shear" section we investigate in detail the consequences for experiments on the mechanical response, in particular for piezorheometry. From our model it turns out that the ground state with constant shear S can become unstable via a stationary or an oscillatory instability.
In the "Influence of two-fluid effects on the central peak" section the focus is on light scattering experiments, but we also address the results obtained from photon correlation spectroscopy and from ultrasmall angle x-ray scattering. In this section we also discuss the issue of transient orientational order.
In the "Summary" section we give the conclusions and a perspective for future investigations in this field, where spatial heterogeneities on many different length and time scales play an important role.
In Appendix we list the complete set of dynamic macroscopic equations for the two-fluid model including transient positional order: the dynamic equations, and the dissipative (associated with entropy productions) as well as the reversible (zero entropy production) currents.

Variables
The hydrodynamics of a simple fluid is described by the momentum density g i , the mass density ρ, and the total energy density ε representing the local conservation laws of the fluid.
As the fluid is cooled down, the situation frequently changes as one approaches, for example, the glass transition, from above: clusters (also denoted as density variations and/or spatial modulations of the density) arise upon cooling. In the following we will investigate a new two-fluid model for a fluid with such clusters.
On the macroscopic level we describe the system as a homogeneous mixture of a solvent part -the original simple fluid -with density ρ s , and parts composed of clusters with a different density ρ c , due to the different packing of the molecules. As it is known experimentally the clusters are in addition associated with a relaxing strain field, U ij introduced in the spirit of polymer dynamics (Brand et al. 1990). The nonlinear Eulerian strain tensor used here has been discussed previously (Temmen et al. 2000;Pleiner et al. 2000;Temmen et al. 2001). It fits into the GENERIC framework (Grmela 2002), is suitable to describe general viscoelastic phenomena ) and describes polymer rheology to a large extent (Müller et al. 2016a;. Throughout the present paper we consider systems that are on average spatially isotropic to keep the presentation as simple as possible. We refer to the recent paper on the breakdown of flow alignment (Brand and Pleiner 2021a) on how one can incorporate anisotropy effects as needed.
In contrast to binary mixtures, the mass densities are not conserved individually but are allowed to exchange by mutual relaxation -in addition to a diffusive mass transport according to two different velocities. The solvent and the cluster mass density, ρ s and ρ c , respectively, add up to the total density ρ = ρ c +ρ s . Similarly, the two velocities give rise to a cluster momentum density, g c i = ρ c v c i , and to a solvent one, g s i = ρ s v s i that add up to the total momentum For details of these two-fluid aspects we refer to Pleiner and Harden (2003).
As additional variables compared to the one-fluid description, we therefore take the relative velocity, w i = v s i − v c i , between the cluster and the solvent velocities, the relaxing strain field U ij as well as the mass concentration of the solvent, φ = ρ s /ρ.
We assume that deviations of the densities relax towards their ground state valueṡ according to the appropriate conjugates μ s and μ c , the two relative chemical potentials. The dots in Eq.
(1) denote all the other contributions. Since the total density ρ, being a conserved variable, must not relax, we have to require the condition μ s /τ s + μ c /τ c = 0. The first law of thermodynamics relates changes of the variables to changes of the energy density ε by the Gibbs relation (Martin et al. 1972;Callen 1985;Pleiner and Brand 1996).
The entropy density σ represents the thermal degree of freedom of the system. The appropriate thermodynamic conjugates are the temperature T , the chemical potential μ, the osmotic pressure Π, the mean velocity v i = g i /ρ, the stress tensor Φ ij , and m i , the conjugate field to w i . The chemical potential is the partial derivative of the energy density with respect to the density (Martin et al. 1972;Reichl 1980;Pleiner and Brand 1996). We have split the strain field and the elastic stress tensor into traces, U kk , Φ kk , and traceless parts, which is common in isotropic elasticity. In that caseΦ ij dŨ ij + 1 3 Φ kk dU ll enters the Gibbs relation. Rotational invariance of the Gibbs relation (2) leads to the conditioñ From the requirement that the energy of the system is a first order Eulerian form of all extensive variables, one gets for the pressure p ≡ −(∂/∂V ) ε dV = −(∂/∂V )E the Gibbs-Duhem relation

Basic hydrodynamics
Hydrodynamics allows to describe the statics and the dynamics of a system in two separate steps. The statics is given by the relation between the conjugate fields and the variables, and the dynamics relates the time derivatives of the variables to the phenomenological currents, which themselves are expressed by the conjugates or their gradients. The static behavior of our macroscopic system is conveniently described by the energy functional from which the conjugate fields follow by taking partial derivatives according to the Gibbs relation, Eq.
(2). The explicit form of the conjugates is listed below. The remaining relations between velocities and momenta v i = g i ρ and are not really static, but nevertheless follow from the energy density, in particular from the kinetic energy density ε kin = (1/2ρ c )[g c ] 2 + (1/2ρ s )[g s ] 2 = (1/2)g 2 + (α/2)w 2 . The w i -dependence of the chemical potential and the osmotic pressure are due to the ρand φ-dependence of α.
From the energy functional, Eq. (5), amended by the kinetic energy, ε kin = (1/2)g 2 + (α/2)w 2 , the static relations are found by partial derivation and read There is a total of six static susceptibilities from the binary mixture fluid (C V , α φ , α ρ , κ φ , κ π , κ μ ), two elastic Hooke-like moduli from the two elastic media (longitudinal c l and transverse c tr ), and three general susceptibilities describing the cross-coupling between the fluid and the elastic degrees of freedom (α u , κ u , and κ ρ ).
The full nonlinear dynamic equations will be given and discussed in the Appendix, Eqs. (A.1)-(A.7). In order to define the phenomenological currents, we list here simplified dynamic equations, where all nonlinear transport and convection contributions are omitteḋ with 2A ij = ∇ i v j + ∇ j v i and R the energy dissipation function. The dynamic equation forε follows from Eqs. (12)-18 making use of the Gibbs relation (2). Note that Eqs. (13) and (14) can equivalently be written aṡ where I D φ contains, among others, the relaxations shown in Eq. (1).
We have split the phenomenological currents into reversible and dissipative parts, denoted by superscripts R and D, respectively. The second law of thermodynamics requires with the equal sign (> sign) for * = R ( * = D). We note that the important role of the behavior of macroscopic variables under time reversal has been discussed in some detail before (Brand et al. 2014;2018).

Ground state under external shear flow
To analyze the response of a fluid with clusters under an external shear flow we assume a linear velocity profile where v z is now the barycentric/mean velocity, because this is the quantity which is controlled in a classical shear flow experiment. Deviations of S will be mostly neglected, and in addition we set v x = 0 and consider a strictly two-dimensional situation. Throughout this section we focus on the limit of small shear rates to address piezorheometry, a technique using very small shear rates (Collin and Martinoty 2003;Pozo et al. 2009). We discuss a linear stability analysis in this limit. Naturally, and in general, in particular for large shear rates, a linear profile can no longer be expected to apply.
To analyze the piezoelectric experiments of Collin and Martinoty (2003) we use as macroscopic variables the concentration variation, δφ, the relative velocity, w i and the component U yz of the strain tensor, which is directly generated by the external shear flow, S. The diagonal components of the strain tensor, U yy and U zz , are directly generated by the shear flow only in quadratic order and will be neglected in the following.
Starting with the Eqs. (14), (16) and (17), the linearized equations of the macroscopic variables considered here take the forṁ In Eqs. 23-26 ζ Π denotes the relaxation rate of Π, the thermodynamic conjugate of the concentration φ, ξ the relaxation rate of m i , the thermodynamic conjugate of w i , and ζ tr the relaxation rate of the mechanical stresscompare also Eq. (A.9) for the introduction of these three relaxation rates in the framework of the dissipation function R. All the static susceptibilities in Eqs. (23)-(26) have been introduced in Eqs. (5)-(11).
It is easy to see that there is an exact homogeneous solution of Eqs. (23)-(26) whilė leads, for a time-periodic shear rate S = S 0 exp(iΩt), to This result is obtained as a homogeneous solution of the linearized macroscopic equations for time-periodic shear and shows that only the product of the relaxation rate of the mechanical stress and the transverse elastic constant enters. We also note the intuitive result that the component U yz is directly proportional to the amplitude S 0 of the shear rate. In the following we focus on time-independent shear.

Linear instability of the ground state with constant shear rate S = S 0
We will now investigate for which values of the external driving force, a static shear S 0 , the solutions Eqs. (27)-(29) become unstable. The goal of a linear stability analysis is then to find out for which values of the shear S 0 the spatially homogeneous solution (Eq. (27)-(29)) becomes unstable to spatially inhomogeneous solutions signaling the onset of spatially heterogeneous cluster formation. At the onset of such an instability the concentration variations, δφ, the relative velocity, w i , and the variations of the strain field, U yz , will start to grow from zero. Frequently one encounters two possible types of instabilities: stationary instabilities, which set-in without an explicit time-dependence, and oscillatory instabilities, which are associated with a finite frequency at onset. Below we will discuss both options in detail.
Along the lines of a linear stability analysis for the onset of a hydrodynamic instability (Chandrasekhar 1961), we start with the plane wave ansatz We look for which values of ω and S we find a solution with χ ≡ 0, where (δφ, w y , w z ,Û yz ) do not decay to zero anymore. As for the ground state, we do not consider the longitudinal strains.
Here we will assume that k y is not fixed externally, while k z has a fixed value to accommodate the boundary conditions at top and bottom plate of the shear setup, e.g., k z = π/L with L the thickness of the layer. The critical value of k y is chosen such that the threshold value for χ = 0 has its minimum possible value.

Stationary instability
We first consider the case of a stationary instability with ω ≡ 0. The solvability condition for Eqs. (23)-(26) for the deviations from the zeroth order solution Eq. (29) reads with ζ 1 = 1 2 ρ 2 ζ Π , ζ 2 = ρ 2 s ζ tr ,ζ = ζ 1 + ζ 2 , and S = 2αS. The critical wavevector k crs y , for which S becomes minimum, has to fulfil In Eq. 32 the superscript crs refers to critical and stationary. What one is looking for in a linear stability analysis is the value of S for which an instability first arises as a function of k. EliminatingS from the coupled Eqs. (31) and (32) a condition for the critical wavevector can be found with the abbreviations f (k z ) = k 2 z − ξζ and g(k z ) = k 4 z + 2k 2 z ξζ + 2ζ 1 ζ 2 ξ 2 leading to the critical values Equations (34) and (35) contain a number of special cases. For example, looking at the limit ξ → 0corresponding to long relaxation times of w i -we obtain k crs y → k z , corresponding to circular rolls, and S crs → 0. Therefore the picture emerges that small relaxation rates of the relative velocity, w i leads to a reduction of the threshold value for the shear rate necessary to trigger a stationary instability.
A similar analysis can be done for the opposite case, where the relaxation rates of δφ and U yz , associated with ζ Π and ζ tr , respectively, are much smaller than the relaxation rate of w associated with ξ . In the limit ζ 1 ∼ ζ 2 → 0 we again get circular rolls (k crs y → k z ), however at a finite threshold value S crs → 2αξ . Nevertheless, increasing ζ 1 or ζ 2 , the threshold S crs increases (as well as k crs y ). It is also instructive to compare this result with the one previously obtained for the breakdown of flow alignment in nematic liquid crystals (Brand and Pleiner 2021a). In this case one has the same macroscopic variables except for the strain U yz . Inspecting Eqs. (67) and (68) of Brand and Pleiner (2021a), we find that the two results (k crs y → k z and S crs → 2αξ ) coincide provided phase diffusion is neglected.
These are important special cases, since they show that there are larger thresholds for larger relaxation rates. As an intuitive picture we thus find small stationary thresholds for small values of the relaxation rates, i.e., for slow temporal changes (or large lifetimes), as one expects for large clusters. The opposite case, large relaxation rates corresponding to fast dynamics (short lifetimes) lead to high threshold values for small clusters.
The picture just outlined can be compared qualitatively with the data presented for the correlation time and the correlation length experimentally for various polymers and low molecular weight compounds (Fischer 1993;Kanaya et al. 1994;Walkenhorst et al. 1998;Patkowski et al. 2000;Patkowski et al. 2001a;Patkowski et al. 2001b;Fischer et al. 2002). Typically either the correlation length or the correlation time were measured as a function of temperature. And the data found compare well with the intuitive picture in the sense that the correlation length/average size of the clusters increases as the temperature decreases and the correlation time/ lifetime of the clusters associated with the ultraslow mode also increases as the temperature decreases. What does not seem to be available in the literature so far is a plot of the correlation time versus the correlation length for one chosen material. In this case a more quantitative comparison with the analysis presented here would be possible.
For the intuitive picture outlined there are a number of similarities to the situation near the onset of stationary thermal convection in a simple fluid, where the threshold for convective onset, the critical temperature difference, T c , scales with the product of the viscous contributions, the kinematic viscosity, ν, and the thermal diffusivity, κ (Chandrasekhar 1961). This is also quite intuitive: for small viscosities one needs a smaller temperature difference to get thermal convection started. This is the case for water and liquid metals. For large viscosities as for silicon oil, glycerol or for honey etc. the required temperature difference to start thermal convection is much higher.

Oscillatory instability
Next we investigate the conditions for an oscillatory instability ω = 0. Now the solvability condition for Eqs. (23)-(26) for the deviations from the zeroth order solution Eq. (29) takes the form of a fourth order polynomial in ω where A = r 1 + r 2 + 2r 3 (37) where we use as abbreviations the inverse relaxation times of, respectively, the concentration φ, the transverse strain U yz and the relative velocity components, r 1 = ζ Π /κ φ , r 2 = ζ tr c tr , and r 3 = αξ , as well as the wave velocities (squared) related to concentration and strain, v 2 1 = α/ρ 2 κ φ and v 2 2 = αc tr /2ρ 2 s . Note that D is related toK defined in Eq. (31) by D = v 2 1 v 2 2K making sure that in the static limit, ω = 0, Eq. (36) reduces to Eq. (31).
Separating the real and imaginary part of Eq. (36) one can eliminate ω 2 from Eqs. (41) and (42) resulting in the condition To evaluate k cro y we have to take the derivative of Eq. (43) with respect to k y and equate it to zero with Multiplying Eq. (44) with −k y and adding it to Eq. (43) one can eliminate the terms ∼ S 2 giving rise to an equation linear in S. The resulting expression for S can then be inserted into Eqs. (43) or (44) to obtain a closed fourth order polynomial equation for (k cro y ) 2 . However this is a very complicated expression that is difficult to interpret. Therefore, we will employ some approximations.
First, we assume that fluctuations of the solvent concentration φ relax on a very long time scale and that the appropriate wave velocity is small which is obtained for a small static susceptibility 1/κ φ .
In that case D ≡ 0, and Eqs. (43) and (44) for the critical quantities reduce to C = AB and C = AB , which are easily solved by where the last relation follows from Eq. (42), reading here ω 2 = B. For D = 0 there is a second solution of Eq. (43), C = 0, which however leads to ω 2 = 0 and is therefore a stationary solution already discussed in the "Stationary instability" section. Similarly, assuming that the elastic modulus c tr (instead of 1/κ φ ) is small, meaning r 2 → 0 and v 2 2 → 0, one gets the solutions Eqs. (51) and (50) with the subscripts 2 replaced by 1 Second, we choose a weaker approximation, r 1 = r 2 ≡ R = r 3 and v 2 1 = v 2 2 ≡ v 2 that preserves most of the general structures of the full solution, but reduces the number of possible different combinations of r 1,2,3 and v 1,2 . Using the general procedure described above we find We note that the critical quantities diverge for R = r 3 or r 1 = r 2 = r 3 , meaning there is no oscillatory instability in that case. This remains true, even if v 2 1 = v 2 2 . Also in the special case v 2 1 = v 2 2 , there is no singularity, if the three relaxation rates are not equal.
The positivity of ω 2 cro requires either R > r 3 or, for r 3 > R, an implicit condition r 3 > r c 3 with r c 3 large enough to reduce the negative contribution in Eq. (54) sufficiently. Thus, there is no oscillatory instability for R < r 3 < r c 3 . The critical shear rate diverges, S cro → +∞ for R − r 3 → 0 + , while it is always negative and finite for r c 3 − R ≥ 0. We emphasize, however, that the fact that in a linear stability analysis the threshold for an instability diverges or that a line of instabilities ends, is well known from pattern forming instabilities. As examples we refer to the onset of thermal convection in miscible binary fluid mixtures (Hurle and Jakeman 1971) and in nematic liquid crystals (Lekkerkerker 1977).
The critical wave vector is obtained from showing explicitly that (k cro y ) 2 is always positive. It also contains the singularity at R = r 3 . We observe that for the approximate case discussed here the equation to determine the critical wavevector is linear in (k cro y ) 2 , while in the general case combining Eq. (43) with Eq. (44) results in a quartic polynomial for (k cro y ) 2 -as already pointed out above.

One-fluid description
As already discussed in the introduction, extensive light and ultrasmall angle x-ray scattering studies have shown that the spectral weight under the central peak is larger than expected on the basis of the hydrodynamics of simple fluids (Fischer 1993;Kanaya et al. 1994;Walkenhorst et al. 1998;Patkowski et al. 2000;Patkowski et al. 2001a;Patkowski et al. 2001b;Fischer et al. 2002). These experiments stimulated work on the macroscopic dynamics of cluster formation above the glass transition (Brand and Kawasaki 2003).
Although the scattering amplitude S(k, Ω) reflects density fluctuations, the thermal degree of freedom is involved due to the thermal expansion effect. Being a diffusive, non-propagating excitation, heat conduction gives rise to the central peak. In a simple fluid it is related to the dispersion relation where K is the thermal conductivity, ρ 0 the total density in equilibrium and C p the specific heat at constant pressure in the general case (Forster 1975;Berne and Pecora 1976). In our notation it follows from the entropy diffusion where κ = K/T 0 with the thermal conductivity K from Eq. (57) and taking into account the static coupling between entropy and density described by the thermal expansion coefficient α ρ , Eq. (7), resulting in ω p = ω V C V /C p . The frequency ω p describes the width of the central peak and governs its height (Forster 1975;Reichl 1980) It was shown (Brand and Kawasaki 2003) that transient positional order (as well as transient orientational order) associated with the formation of clusters contribute to macroscopic dynamics and give rise to additional weight under the central peak. Furthermore, it was pointed out (Brand and Kawasaki 2003) that so-called isotropic biaxial nematic order (Mermin 1979) of the cubic or icosahedral type (Liu 1981) also contributes to the central peak.

Two-fluid description
In the present section we evaluate how two-fluid effects of the nature studied in the "Two-fluid model for a fluid with clusters" section and in the appendices can contribute additional spectral weight under the central peak. When compared to previous work (Brand and Kawasaki 2003) we have as additional macroscopic variables the concentration φ and the relative velocity w i . First, we evaluate how static and dissipative contributions associated with two-fluid effects can affect the central peak.
As for the static contributions related to concentration variations δφ we have in the energy In addition to the diagonal susceptibility associated with concentration variations, 1/κ φ , we have cross-coupling terms to the scalar variables of a simple fluid, density variations, δρ, and entropy variations, δσ . Furthermore there is a static cross-coupling term to the diagonal strains of the transient network. Due to the static cross-coupling terms all terms in Eq. (60) contribute to the central peak.
For the dissipation we obtain from the dissipation function given in Eq. (A.9) the contributions where Π is the thermodynamic conjugate of φ. In Eq. (61) the contribution ∼ ζ Π is associated with the relaxation of the concentration and the contribution ∼ D with diffusion. They will typically both be of interest here since the clusters observed via light scattering go in size up to ∼ 2000Å and are of large, but finite spatial extent. This effect can also be seen immediately from the fourier transform being ∼ (ζ Π + Dk 2 ) for these two diagonal contributions. We see from Eq. (61) that concentration gradients couple dissipatively to temperature gradients and thus to the central peak via heat diffusion as well as to both, diagonal strains and strain gradients.
As for the second macroscopic variable associated with two-fluid effects, the velocity difference w i , the coupling to the central peak is a bit more subtle in nature.
From Eqs. (8) and (9) we obtain for the two relevant thermodynamic conjugate fields Π and the chemical potential μ From inspection of Eqs. (62) and (63) we see that both, Π and μ, pick up contributions quadratic in the velocity fields. These in turn will give rise only to corrections to the central peak, which are of higher order than those studied for δφ in this section so far.
For the dissipative contributions involving the relative velocity we have from Eq. (A.9) with m i the thermodynamic conjugate of the relative velocity w i . The only cross-coupling of w i is to A ij = 1 2 (∇ i v j + ∇ j v i ). This will lead for the central peak to contributions from the relative velocity w i , which are again of higher order compared to those of the concentration variations.
Therefore, neither in the statics nor in the dissipative dynamics, does the velocity difference w i couple to the variables relevant for the central peak in lowest order. However, there is a reversible dynamic coupling to the concentration variable and therefore indirectly to the entropy, Eqs. (A.4), (A.6), (A.13) and (A.14), Taking the strong damping limit, ω αξ , Eq. (66) simply leads to a replacement D → D+1/(ρ 2 0 ξ) in Eq. (65) thus adding to the spectral weight under the central peak.
In the scattering amplitude there are, apart from the central peak, also the two Brillouin peaks, centered at finite frequencies, that reflect propagating sound modes. Amending the one-fluid description by a Maxwell-type relaxation of the viscosity, the Brillouin peaks are changed (Boon and Yip 1980). Since in our two-fluid description the relaxing strain as well as the relative velocity couple dynamically to the momentum density, the Brillouin peaks will be affected as well. This is, however, beyond the scope of the present paper.

Orientational order
So far we have concentrated throughout this manuscript on two-fluid effects for the strain field as an additional macroscopic variable. In parallel to previous work (Brand and Kawasaki 2003) we briefly outline possible effects on the central peak associated with transient orientational order and with truly biaxial nematic order of optically isotropic symmetry.
For transient orientational order in an isotropic phase the analysis proceeds in parallel to that of transient positional order by replacing the strain as variable with Q ij , the traceless symmetric order parameter characterizing nematic orientational order (de Gennes 1971;1975;Pleiner et al. 2002). There is one linear static coupling term to δφ, manifest in the free energy as For the thermodynamic force ij associated with Q ij that gives rise to where . . . indicate the expression previously derived in (Brand and Kawasaki 2003). There are two diagonal terms in the energy complementing the cross-coupling given in Eq. (67). One is the analog of Frank's free energy for Q ij , which reads (de Gennes 1971;1975) and the other one takes the form K φφ (∇φ) 2 . Thermostatic stability then requires for ζ φ in the one constant approximation for ε Q the inequality This inequality shows that ζ φ is bounded from above by the two diagonal coefficients, but its sign is not fixed: it can be positive or negative. The static coupling term given in Eq. (67) is a higher order gradient term of the type discussed in Pleiner and Brand (1980). In this reference we discussed -in addition to other effects -static cross-coupling terms between gradients of the nematic director field on the one hand and gradients of the density and the entropy density on the other. We should also point out that, of course, the term 1 2 κ φ (δφ) 2 exists in the energy density (compare Eq. (5)), which is of lowest order in k.
In contrast to the case of transient positional order, cf. the contribution ∼ 1/κ u in Eq. (5), the static coupling in Eq. (68), is not gradient-free. Nevertheless, it can contribute significantly to the central peak, since from light scattering we know that length scales of ∼ 10 → 100 nm play an important role.
In addition, there are dissipative linear cross-couplings, which read in the dissipation function which can also contribute to the spectral weight under the central peak as a higher gradient effect. We will briefly comment on the influence of more general types of orientational order on the central peak. First, for permanent biaxial nematic order of lower symmetry, such as tetragonal or orthorhombic, optical anisotropy will arise, which is not compatible with the light scattering observations of the central peak. On the other hand, there is the so-called isotropic biaxial nematic order (Mermin 1979) of cubic or icosahedral type (Liu 1981). As has been shown in simulations (Tomida and Egami 1995), the latter case can occur above the glass transition, at least locally. The order parameter for cubic biaxial nematics is a fourth rank tensor (Nelson and Toner 1981), while for icosahedral biaxial nematics it is a sixth rank tensor Steinhardt et al. (1981a, b).
As variables for those biaxial nematic phases one can use three Eulerian angles, Θ i (Brand and Pleiner 1981a;Liu 1981). Then, the coupling to concentration variations proceeds via the density or entropy density dependence of the Frank elastic constants (Brand and Kawasaki 2003).
In closing this section we point out that there is another class of systems for which transient positional order as well as transient orientational order play an important role for the understanding of the macroscopic properties, namely the sponge or L 3 phase in lyotropic liquid crystals close to the phase transition to the isotropic liquid phase (Pleiner and Brand 1991;.

Summary
Using the macroscopic dynamics approach we have studied in this paper a two-fluid model for cluster formation. We have focused predominantly on the description of experimental results above the glass transition. The experimental results analyzed cover cluster sizes between ∼ 10nm and ∼ 20μm. They have been mainly obtained by optical techniques such as light scattering from the central (Rayleigh) peak and photon correlation spectroscopy as well as ultrasmall angle x-ray scattering for length scales up to ∼ 2000 A. The second large group of experiments used two types of mechanical experiments, namely piezorheometry and anelastic mechanical measurements. Length scales for clusters of up to ∼ 20μm are covered with frequencies in the range from 10 −2 H z to 10 3 H z. For the light scattering experiments we find that the additional macroscopic variable, concentration, adds additional weight to the central peak due to all its cross-coupling terms to total density, temperature fluctuations and strains. We also obtain the result that the relative velocity contributes to light scattering from the central peak mainly to higher (quadratic) order. As for light scattering we find that both, transient positional order as well as transient orientational order, contribute additional extra weight to the central peak. To describe the results of the mechanical experiments we have studied the influence of a constant shear. We find that both, a stationary instability and oscillatory instability, for the onset of cluster formation are possible. So far the experimental results, in particular on piezorheometry, suggest so far the onset as a stationary instability in the low frequency limit. It might be worthwhile, however, to check experimentally whether low frequency oscillations can be found as well.
In the main body of the paper we have also pointed out in the "Stationary instability" section some similarities between our analysis of the influence of clusters with the breakdown of flow alignment. As a perspective we mention that it appears worthwhile to investigate in the future the analog of stress-oscillations in shear start-up and of backflow upon cessation of shear for a two-fluid system of the type considered here. This could be done along the lines of Müller et al. (2016a, b), where these processes have been studied in the framework of nonlinear transient elasticity for the one-fluid case. Such an extensive study is, however, beyond the scope of the present paper.
As already mentioned in the introduction, the issue of fluid-fluid transitions (Takae and Tanaka 2020) is clearly also an interesting subject for the macroscopic two-fluid description presented in the present paper. It will therefore be most interesting to see whether some of the predictions made here re. additional spectral weight under the central peak and the formation of large clusters as monitored by piezorheometry can be verified. So far we are not aware of such experimental investigations for these systems. Another class of systems of high interest in the present context of two-fluid macroscopic dynamics is of the nature studied recently by the group of Carsten Tschierske (Reppe et al. 2020). These systems, showing various isotropic phases as well as cubic liquid crystalline phases of several types, could provide a handle to study the influence of chirality on macroscopic two-fluid effects.
For the two-fluid effects associated with optically isotropic orientational order we have concentrated here on the influence of long range biaxial nematic order and quadrupolar transient orientational order to address this issue. We note, however, that tetrahedral or octupolar order of short or long range might also be of interest in the context of two-fluid effects. That tetrahedral order could arise for liquid crystals has been pointed out first by Fel (Fel 1995). Tetrahedral order is especially interesting since it breaks parity symmetry. In the following phase transitions (Radzihovsky and Lubensky 2001) as well the macroscopic dynamics of isotropic single-component tetrahedral liquid crystals have been investigated Brand and Pleiner 2010;2017), for a recent review we refer to Pleiner and Brand (2016a). As for the question whether twofluid effects related to tetrahedral phases are important, we are not aware of any experimental literature, neither above the glass transition nor for tetrahedral systems not showing a glass transition. We note, however, that simulations on densely packed tetrahedra by Sharon Glotzer's group (Haji-Akbari et al. 2009) indicate that there can be ordered tetrahedral regions embedded into disordered tetrahedra. Therefore we have addressed the issue of two-fluid effects for tetrahedral/octupolar phases with different types of order quite recently (Brand and Pleiner 2021b).
Finally we mention an open problem re. the glass transition for systems with clusters. While the mode coupling theory for the glass transition of spatially homogeneous phases composed of small molecules has been worked out in detail -compare the review by Goetze and Sjögren (1992), it will be interesting to see whether mode coupling calculations in the same spirit can be carried out for spatially heterogeneous phases with clusters on many different length scales.

A.1 Dynamic equations
The full dynamical equations for the strain, the concentration φ and the other macroscopic degrees of freedom areε with 2A ij = ∇ i v j + ∇ j v i , and R the energy dissipation function. These equations follow from Pleiner and Harden (2003), Pleiner andHarden (2004), andPleiner et al. (2020) and contain, apart from the reversible (superscript R) and irreversible, dissipative (superscript D) phenomenological currents, also transport and convection whenever appropriate. The latter are reversible and, indeed, all transport contributions (including the isotropic pressure) add up to zero entropy production, R = 0. This is automatically achieved by using the mean velocity for all those transport contributions (as in a one-fluid description) in Eqs. (A.1)-(A.7). In a two-fluid description; however, there exist (for most of the dynamic equations) reversible currents, cf. Section A.3, that modify the transport and convection contributions. The resulting actual transport terms lead to zero entropy production, since the contributions from the reversible currents individually fulfil R = 0 due to suitable counter terms. The contributions from the reversible currents are nonuniversal, but come with phenomenological (reversible) transport parameters reflecting the fact that in a two-fluid description there is no general law that fixes the transport velocities (Pleiner and Harden 2003). On the other hand, the procedure described above gives the possibility, by choosing the reversible transport parameters appropriately, to adjust the general result to some desired simplified models, guaranteeing zero entropy production of transport and convection terms. At the end of Section A.3 we show how this procedure works for the present case.
Here is also a good location to discuss the symmetry properties of the stress tensor and the closely related question of the conservation law for angular momentum. These issues have been of interest for a long time. Fields of interest ranged from field theory (Belinfante 1939) over fluids, solids and liquid crystals (Martin et al. 1972) to the superfluid A-phase of 3 He (Pleiner 1977). In these references, in particular in Appendix A of Martin et al. (1972), it has been shown that angular momentum conservation for intrinsic as well as extrinsic angular momenta can be obtained by a stress tensor, which is the sum of a symmetric part and the total divergence of an antisymmetric tensor of third rank.
The concentrations of the clusters (and the solvent) are treated as relaxing quantities, rather than as conserved ones, which is manifest by the form of the current I D φ in Eq. (A.4). The energy conservation law, Eq. (A.1), is redundant due to the Gibbs relation, Eq. (2).
In the whole set of dynamic equations the mean velocity v i has been chosen as the convective velocity for all variables. This ensures zero entropy production of the convective derivatives. Due to various material dependent contributions in the reversible currents (see below), the actual convective velocities can be different from v i and can be specific for the different variables.
For the phenomenological parts of the currents the second law of thermodynamics requires with the equal sign (> sign) for * = R ( * = D).

A.2 Dissipative currents
The dissipative parts of the currents introduced above can be deduced from a potential, the dissipation function R, that reads in bilinear approximation where the rank-4 tensors have the same form as the elastic tensor, e.g., ζ ij kl = ζ l δ ij δ kl + 1 2 ζ tr (δ ik δ jl +δ il δ jk − 2 3 δ ij δ kl ), containing two independent terms denoting transverse and longitudinal contributions. The requirement of positivity of the entropy production associated with the second law of thermodynamics leads to a number of positivity requirements and bounds for the dissipative coefficients entering Eq. (A.9). All the diagonal coefficients must be positive leading to ζ Π , ζ l , ζ tr , ξ, κ, D, ξ l , ξ tr , ν l , ν tr , ν (w) l and ν (w) tr all being positive. In addition, all cross-coupling terms associated with the non-diagonal coefficients are bounded from above by products of diagonal coefficients to guarantee positivity of the entropy production leading to the inequalities (D (T ) ) 2 < κD, ζ 2 φ < ζ Π ζ l , The general expression for strain diffusion is ξ ij klmn (∇ m Φ ij )(∇ n Φ kl ), containing five coefficients (Mason 1958): We have approximated this by using ξ ij kl = δ mn ξ ij klmn in order to simplify the currents. This approximation seems to be justified, since there is already strain relaxation (∼ ζ ij kl ).
From Eq. (A.9) the following dissipative currents are obtained, with ∇ i ∇ j = ∇ i ∇ j − 1 3 δ ij ∇ 2 k . Diffusion and thermodiffusion is written in the usual way with D = ρ d and D (T ) = α d (T ) . For the relaxation of φ we can similarly write ζ Π = α/(ρτ φ ) where τ φ = ρ 2 (τ c + τ s ). The relative velocity, w i , always relaxes, since it is not related to any broken symmetry, nor to a conservation law.
The same reasoning applies to the variables φ and U ij . All other variables are conserved and show diffusional behavior.

A.3 Reversible currents
For the reversible parts of the currents in Eqs. (A.1)-(A.7) we find, requiring zero entropy production j (σ )R i = β 1 m i , (A.17) σ R ij = 2β 2 m i w j + β 8 δ ij Π (A.18) where we have kept nonlinear contributions when they contribute to the transport and convection of variables. For the reversible coefficients in Eqs. (A.17)-(A.21) either sign is possible and no bounds can be given in general from the thermodynamic requirements associated with the second law of thermodynamics. These reversible currents can be used to tune the transport and convective velocities of the different variables, compare Pleiner and Harden (2003, Sec. 6). In particular, the β 1 term adds to the transport of the entropy density, Eq. (A.2), but since we assume that the mean velocity, v i , is the appropriate transport velocity, we put β 1 = 0. In our model the transport velocity ofρ c andġ c i is v c i , while forρ s andġ s i it is v s i . To achieve this we have to take γ = 0 = β 8 , cf. Eqs. (19) and (20), and β 2 = β 4 = 1/2 and β 3 = 1/ρ s − 1/ρ c Pleiner and Harden (2003). Similarly, we assume that translation, transport, and convection of Φ ij , Eq. (A.7), is done by the velocity v c i , which requires β 5 = −β 6 = −β 7 = ρ −1 φ −1 .