Phase separation and nucleation in mixtures of particles with different temperatures

The difference of activity in colloidal particles is sufficient to drive the phase separation in active and passive (or less active) mixtures having even only excluded volume interactions. In this report, we study the phase separation kinetics of such systems in the diffusive limit. This model is described by a mixture of diffusing particles coupled to different thermostats, thus having a non-equilibrium nature due to temperature differences. However, we show that indeed the system recovers equilibrium thermodynamics in the dilute limit. We obtain phase diagrams showing asymmetry in concentrations due to activity differences. By using a more general approach, we show the equivalence of phase separation kinetics with the well known Cahn-Hilliard theory. On the other hand, higher order expansions indicate the emergence of non-equilibrium effects leading to breakdown of the equilibrium analogy as well as posing a problem to define a scalar temperature. We lay out the general theory in terms of accessible parameters which we demonstrate by several applications. In this simple formalism, we capture interesting scaling laws for interfacial properties, droplet growth dynamics, and phase segregation conditions some of which equally observed in simulations.


I. INTRODUCTION
Many-body systems self-organize and exhibit macroscopic collective behavior, which is amenable to coarsegrained dynamics [1]. The principles of equilibrium thermodynamics which provide the essential tools for studying equilibrium self-assembly are not adapted to the new phenomena emerging in active systems [2]. While active systems are non-equilibrium in their nature, with each constituent having its own energy budget and objective, it is quite fascinating how these systems can display distinct characteristics which are sometimes non-trivially related to the known basic physical principles. A typical example is that of active phase separation.
In living organisms, activity-driven phase separation plays an important role by promoting self-organization and increasing the efficiency of biological functions inside cells [3,4], which are due to operate in a very crowded and noisy environment. The constant use of energy in unequal amounts by individual cellular sub-components reflects their varieties in dynamical and chemical activities [5] along with their innate physiological differences. In turn, these activity differences may enhance phase separation by creating effective attractions between alike components. Remarkably, similar characteristics and governing principles are observed at various length scales on diverse biological systems and they are now also guiding synthetic model systems [6][7][8][9][10].
It is now well known that activity differences in colloidal particles is sufficient to drive phase separation between active and passive (or less active) hard spheres having only excluded volume interactions. These systems can be considered to have particles with two different temperatures which mimics the effect of activity differences [11][12][13]. In principle, the constituent particles exchange energy, the hotter ones providing extra energy for the colder ones. This, in turn, indicates a non-equilibrium behavior [14]. Similarly, in polymeric systems, tiny activity differences in active/passive polymer mixtures enhance phase segregation [15,16]. If the system is diffusive, the effective temperature can simply be deduced from the effective diffusivity. For instance, at room temperature T , bacteria may reach an effective temperature T A ∼ 10 2 − 10 3 T in translational motion by chemotaxis [17]. These realizations are not limited to biological or polymeric systems, but also seen in plasma physics [18] and in the description of thermal phases in the interstellar medium [19,20]. Due to having different heating/cooling processes, ionized hydrogen reaches a temperature ∼ 10 2 times larger than atomic hydrogen distributed in the interstellar medium [21]. In many other contexts, the concept of effective temperature [22] proves to be quite useful in understanding large-scale phenomena that originate from microscopic motility or even chemical differences between the building components [23][24][25][26]. Considering the growing interest in these systems, it is important to study and establish a theory of phase separation starting from the microscopic dynamics.
In this work, we extend the theory of mixtures of particles with different temperatures introduced in Ref. [11] to inhomogeneous mixtures and we study the phase separation kinetics of solutions composed of two diffusing species, which are coupled to different thermostats. Starting from the microscopic model (Sec.II), we derive and present the theory in the dilute limit (Secs.II-IV) and obtain the scaling laws for phase behaviors. Even though the system has inherently non-equilibrium properties due to activity differences, it obeys an effective equilibrium thermodynamics in the dilute limit including the interfacial contributions, which generalizes the Cahn-Hilliard theory of solutions [27] to two temperature mixtures. In Section V, we demonstrate example applications of the theory in diverse systems. In Section VI, we consider the higher order corrections in concentrations which break the equilibrium analogy and raise the difficulty to define a scalar temperature. For purely hard-core interactions, we present higher order corrections in density to the theory by considering depletion interactions between spheres of different activities, and obtain results, which do not seem to be easily obtainable otherwise.

II. MICROSCOPIC MODEL
We start by introducing the microscopic model and we derive an effective thermodynamic description, which is valid when the solutions that we study are very dilute. We follow the lines of Ref. [11] and study a solution of mixed particles satisfying overdamped Langevin equations, in contact with reservoirs at different temperatures where U is the overall interaction potential between the particles. The position of particle m is denoted by x m , its friction coefficient by ζ m , its temperature by T m and ξ m (t) is a standard zero mean, unit variance, Gaussian white noise. Here, we consider two different species of particles each being in contact with a thermostat at temperature T A or T B . This dynamics presented as a Langevin equation for each particle in the system, can be reformulated as a Fokker-Planck equation for a multiparticle probability distribution P ({r}) where {r} is a vector whose components are the positions of all the particles. The Fokker-Planck equation is written in terms of the fluxes J m for each particle: Distinguishing the two species of particles A and B and integrating over all coordinates except for one, we obtain the single particle distributions p α (r 1 ) where now α, β = A or B [11]: ∂p α (r 1 ) ∂t = Tα ζα ∇ 2 r1 p α (r 1 ) + 1 ζα ∂ r1 p α (r 1 ) βf αβ , and N β is the number of β-type particles. We may write the two particle densities in terms of single particle distributions and a pair distribution function p αβ 2 (r 1 , r 2 ) = p α (r 1 )p β (r 2 )g αβ (r 1 − r 2 ). Accordingly, g αβ (r 1 − r 2 ) is the pair distribution function and in the long time limit, it reaches a steady-state value, i.e., g αβ (r 1 − r 2 ) → g ss αβ (r 1 − r 2 ). The knowledge of the steady state pair distribution function is sufficient to close this set of equations. In general, g ss αβ (r 1 − r 2 ) depends on the particle concentrations and can be expanded in powers of the concentrations. In the dilute limit, when considering only pair interactions, the solution of Eq. (2) imposes: where T αβ are the mobility-weighted average temperatures and are defined as T AA ≡ T A , T BB ≡ T B , and In order to determine the forcesf αβ , we rewrite the two-particle distribution function as p αβ 2 (r 1 , r) = p α (r 1 )p β (r 1 + r) exp[−u αβ (r)/T αβ ] by defining r 2 = r 1 + r. This change of variables inside the integrals yields Assuming that the concentrations vary slowly over a length scale which is much larger than the range of the pairwise interactions, we expand p α (r 1 + r) ≈ p α (r 1 ) + r · ∇p α (r 1 ) + 1 2 (r · ∇) 2 p α (r 1 ). While inserting in the integrand, the first term gives the uniform or homogeneous contribution to the force, the second term vanishes, and the third term provides an inhomogeneous contribution to the force.

A. Effective thermodynamic identities
We introduce the concentrations c α (x) = N α p α (x), and obtain closed equations for the concentrations: We have defined here the total mean forcef α = βf αβ acting on a particle of species α due to all the other particles. In this particular case, this total mean force is the gradient of a potential and we can write the conservation equation for the concentrations in the Cahn-Hilliard form: This equation defines the functions µ α as non-equilibrium analogs of chemical potentials.
We decompose the non-equilibrium chemical potentials as sums of a homogeneous part, which depends only on the concentration and a non-homogeneous part, which depends on the concentration gradients.
The quantities B αβ = (1 − e −u αβ (r)/T αβ )dr are the effective excluded volumes or second virial coefficients and and Λ αβ = 1 The non-equilibrium chemical potentials can themselves be calculated as the functional derivatives of an effective non-equilibrium free energy µ α = δF /δc α which is the functional derivative of the total free energy F [c A , c B ] = f dr with respect to the concentration c α (r). The reconstruction of free energy from the chemical potentials gives the free energy per unit volume which is given by where L αβ = −T αβ Λ αβ is negative. The free energy f 0 for uniform concentrations has already been derived in Ref. [11]. It has a Flory-Huggins form with differences in interactions dictated by the two different temperatures and the effective excluded volumes B αβ 's (second virial coefficients). The contrast in temperatures further enhances the tendency toward demixing that is inherent to the Flory-Huggins free energy. For instance, in the case where all the B α are identical (B A = B AB = B B ), (which corresponds in particular to equal-sized hard spheres), solely the difference in temperatures drives asymmetrically weighted interactions between particles that are responsible for phase separation. A crucial remark is that the friction coefficients ζ α or the mobilities, which are the inverse of the friction coefficients, only enter through the effective cross temperature T AB which becomes indistinct for T A = T B . Hence, a difference in diffusivities D α ∝ T α /ζ α at the same temperature T A = T B but ζ A = ζ B has no influence on the thermodynamics. This is expected a priori since in this case the system is at thermal equilibrium. It is remarkable that the concept of effective free energy can be extended to inhomogeneous solutions of particles with two different temperatures. In dilute limit, at lowest order in the concentration gradients, this shows compatibility with the Landau-Ginzburg theory [28]. A major difference in this case is that the steady-state is maintained at the expense of extra power input [29].

B. Internal stress
In order to derive the local stress tensor σ ij in the solution, we first calculate it from the free energy as could be done in an equilibrium system and then give a mechanical derivation based on the Irving-Kirkwood formulation of the stress, which leads to the same results.
a. Effective thermodynamic construction: In a deformation of the volume V + δV , the total work is obtained by integrating the work across each surface element. If we call dS j the surface element and u i the infinitesimal displacement along respectively i− and jdirection in Cartesian coordinates, the work associated to the deformation of the volume is On the other hand, if there is an effective free energy, the work done by such a displacement is δW = δF . Hence, we may obtain the stress tensor σ ij from the surface contribution to the change in free energy δF which can be written as: The first term on the right-hand side gives the surface integral f (u · dS) while the second one can be converted to a surface integral by expanding δf in terms of of the changes in the concentrations δc α induced by the deformation, i.e., We study here the stress in a steady-state, and as discussed in the previous sections, the chemical potentials µ A and µ B are constant throughout the volume. We can therefore eliminate the changes in concentrations δc α by using the conservation of the total numbers of particles A and B during the deformation V δc α + δV c α = 0. The total change in the free energy can then be written as a surface integral Finally, the variation of the concentration on the surface is given by δc α = −u·∇c α for both species. The resulting stress tensor reads: Accordingly, the pressure p can be deduced directly from the diagonal component of the stress, pδ ii = −σ ii which gives in three dimensions: The first term contains the locally uniform pressure p 0 = µ 0 A c A + µ 0 B c B − f 0 that is given by the standard Gibbs-Duhem equation and the gradient terms including the contributions of µ ∇ A , µ ∇ B and f ∇ determine the interfacial contributions.
b. Irving-Kirkwood method: An alternative, more general method to calculate the stress tensor without any reference to the equilibrium thermodynamics, has been proposed by Irving and Kirkwood [30], starting from the mechanical virial equation [31]. The stress tensor σ (v) ij is given by: where σ K ij is the stress in the absence of interactions (for an ideal gas), while the second part is the contribution to the stress due to interparticle potentials that we name σ u ij . The interaction part of the stress can be decomposed into a sum over the particle species α and β and the average in Eq. (17) can be calculated using the two-particle probability distribution p αβ 2 (r 1 , r). As in the previous paragraphs, we expand the two-particle densities around r 1 , and make a change of variables to calculate the average values. This leads to the following stress, The rest of the calculation is straightforward. We give the full result of this calculation in Appendix B. It turns out that the stress tensor calculated by the Irving-Kirkwood method is different than the stress calculated from the free energy; it can be written as σ ij is the stress obtained in the previous paragraph from the free energy Eq. (15). This shows that the stress is not defined in a unique way [32]. However, it conserves all the properties of σ (f ) ij for our analysis, and strictly does not alter the force-balance since ij , could just be obtained from a free energy perturbation by adding surface terms to the free energy given by Eq. (10a).
This result validates that in the limit of low concentrations, we can still use the thermodynamic approach to calculate the stress inside the solution. In the more general case where there is no effective free energy, one would need to rely on the Irving-Kirkwood description. A final note is that an alternative formulation of Irving-Kirkwood method can be achieved by using microscopic force-balance [33,34] in the time evolution equations Eqns (6), thus summing up all mean internal forces.

III. PHASE LINES AND THE CRITICAL POINT
We now use the effective thermodynamic description of the solution to calculate the phase diagram of a solution of particles at two different temperatures.

A. Non-dimensionalization of effective thermodynamic properties
Let us introduce first the volume fractions φ α = c α B α /ǫ α where ǫ α ≡ B α /v α is the conversion factor to molecular volume v α . The volume factions are well defined only if the total volume fraction is smaller than one. They must then satisfy φ A + φ B ≤ 1. We also define . Accordingly, we set the dimensionless free energy densitŷ B v A f while the total free energyF = f dr. As a result, we obtain the non-dimensionalized chemical potentialsμ α = δF /δφ α : where we ignored the density-independent terms. Note thatμ B has an extra scaling factor α v in order to conserve all the functional properties to construct the thermodynamic functions of the previous section. Similarly for pressure we havep = T −1 B v A p. This completes our transformation to dimensionless functionals X[c A (r), c B (r)] →X[φ A (r), φ B (r)]. As in the previous case, we can separate these into locally uniform and interfacial part, i.e.,X =X 0 +X ∇ .

B. Two phase coexistence
At zero-flux steady state for single particle concentrations, we should have uniform chemical potentials and pressure. For a mixed state this would suggest to have uniform concentrations (a single phase). However, if there are any two phases coexisting, they should satisfy the following conditions at their interface: where a and b denote the two coexisting phases. Together, these suggest no net particle exchange and forcebalance at phase boundary while concentrations continuously vary from one phase to the other. We obtain the coexistence curve (or binodal line) in our phase diagrams by numerically solving the above conditions.

C. Spinodal line
The stability of the uniform state φ A (r) = φ 0 A and φ B (r) = φ 0 B can be analyzed by linearizing ∂φ A /∂t and ∂φ B /∂t around the uniform state by introducing and φ 0 B denote the uniform states. As a result, we obtain in Fourier space the equation of the relaxation of a perturbation of wave vector q, δφ(q) = (δφ A (q), δφ B (q)) with the relaxation matrix: where κ −1 p is the inverse of the compressibility matrix obtained from the linearization of the chemical potentials and we have neglected terms of order q 4 . Note that in κ −1 p and for the remainder, we no longer write the superscript zero of φ α 's. The instability occurs when at least one of the eigenvalues of Γ becomes negative. Thus, it is enough to determine the limit where the determinant |κ −1 p | is negative. The non-dimensionalized spinodal line equation is obtained as for vanishing q. The instability occurs when ν s < 0. The spinodal line would be symmetric if ǫ A = ǫ B . It is clear that a larger contrast of activity, i.e., a larger α T makes the unstable region of the phase diagram larger.

D. Critical point
The critical point is calculated by finding the point where the two phases in equilibrium are identical [35]. This is the point along spinodal line where the fluctuations are maximum. Hence, we search for the point along the spinodal line where the gradient of the spinodal line in the volume fraction parameter space ∇ν s = φ A ∂ φA ν s + φ B ∂ φB ν s is aligned with the eigenvector e 0 of the inverse compressibility matrix, corresponding to the eigenvalue ǫ = 0. Accordingly, the volume fractions at the critical point (φ * A , φ * B ) satisfy: We observe that even at ǫ A = ǫ B where the spinodal line is symmetric, the location of the critical point can be shifted along the spinodal line by controlling the ratio α T /α v . One easy way to see that is to use conjugate of Eq.(25) and symmetrize these two forms for φ A and φ B . Accordingly, choosing α T > α v moves the critical point toward the B-rich part of the phase diagram while setting α T < α v moves it toward the A-rich side. However, if ǫ A = ǫ B , it seems more complicated to get a sense on such symmetrization and one should follow Eq.(C2). We will investigate further the properties of such asymmetric phase diagrams in Section V in order to get a hint on the structure of coexisting phases that can be either both liquid-like phases or a solid-like and a gas-like phase.

E. Condition for existence of phase separation
In the previous paragraphs, we outlined the calculation of phase diagrams using the effective free energy obtained in the limit of low concentrations, but imposing no constraints on the volume fractions and assuming that both phases remain fluid. If the volume fractions are high enough in one of the phases, this phase cannot be fluid and is solid. There is in this case equilibrium between a liquid or gaseous (very dilute) phase and a solid phase for which our concentration expansion is not accurate but approximate. Still, it qualitatively indicates a crystalline phase. Simulations show that this phase exhibits hexagonal packing in 2-dimensions [12], while both facecentered cubic and hexagonal-close packed structures in 3-dimensions [36]. On the other hand, for hard-sphere interactions and considering that the two phases remain fluid, we improve our approximation in Section VI by adding one order in concentration.
A general constraint on the volume fractions is that the total volume fraction is not space filling and that it is smaller than the critical concentration for space filling φ max (for a random packing of identical spheres φ max ≃ 0.64). Here, for generality, we stick with φ max = 1 following the general literature. This choice, together with low density approximation is sufficient to observe qualitative tendency to phase-separate while changing interaction and activity parameters.
A common approach to study the conditions for phase separation [37] is to impose this space filling condition for the spinodal line given by Eq. (24). This conjecture accepts the emergence of an instability region in the phase diagram as the sufficient condition for phase separation. An alternative more restrictive view would be that the critical point should exist inside a physical phase diagram. In this case, the valid condition is the existence of the critical point inside the physical regime. These two approaches are identical if the critical point is at the tip of the spinodal line while the latter condition delays the onset of the coexistence region. We evaluate both conditions for the cases we consider (Fig. 3).
profiles are in normal density coordinates with ψ(z) and η(z) where we see that η(z) ≪ ψ(z). The interface width is 2ℓ and away from the interface the densities take constant values at two coexisting phases.

IV. PHASE SEPARATION KINETICS
Knowing the coexisting phases and the effective thermodynamic description, we can develop the theory of phase separation kinetics for two temperature mixtures. We start here by determining the surface tension between two phases at equilibrium and then discuss phase separation kinetics. Note that the interface between the two phases is stable only if the surface tension is positive.

A. Surface tension
We consider a mixture with two phases at equilibrium and with a flat interface between the two phases. In this geometry, the two concentrations or the volume fractions vary only along one direction, say the z-direction so that ∂ y c α = ∂ x c α = 0. The stress is isotropic in the bulk of each phase but becomes anisotropic close to the interface. The interfacial tension between the two phases can be calculated from the stress distribution in the solution where the integration is from one phase (phase-a) to the other (phase-b). Using our calculated value of the stress tensor σ ij , Eq. (15), we see that the isotropic component of the stress proportional to δ ij cancels out and that only the non-diagonal components of the stress contribute to the surface tension: they vanish for σ xx but do not vanish for σ zz . In dimensionless form, In order to determine the surface tension γ from Eqns (26) and (27), we then need to determine the concentration profiles along z-direction. We first introduce the boundary conditions in the two phases in equilibrium: where the values with daggers are the constant values of the chemical potentials and the pressure at equilibrium. The concentration profiles can then be calculated from the coupled differential equations: where the tilted free energy is given by: While this is the generic form, the same treatment more specifically implies that ∆f = 2∆f 0 . The tilted free energy is the difference between the local free energy along the concentration profiles and the energy obtained from the so-called common tangent construction. Since the common tangent construction gives the minimal possible free energy, the tilted free energy is always positive ∆f > 0. We may therefore write an alternative form of the interfacial tension : The surface tension γ is therefore always positive and the interface between the two phases is stable. The set of equations (28) can be solved numerically by linearization of the two equations around one boundary (say phase a) and integrating up to the other boundary (phase b) using a shooting method. Alternatively, an analytical approximation can be obtained by considering the system close to the critical point as done in the next paragraph. This analytical approximation is in excellent agreement with the numerical results.
Surface tension near the critical point: In the vicinity of the critical point, we show in Appendix C how the effective thermodynamics can be expressed as function of a single order parameter ψ which is a linear combination volume fractions relative to the critical point. In addition, the other normal coordinate η gives the normal distance from the critical point, and a phase separation occurs when a solution exists η a ≈ η b > 0. Each value of η a defines the two coexisting phases ψ a and ψ b .The transformed coordinates are illustrated in Fig.1, where at first the concentration profiles φ A (z) and φ B (z) are obtained by solving Eq. (28) numerically. The effective tilted free energy density is given in (C11) as a function of the order parameter only, which is obtained by transforming (29) to normal coordinates. Minimization with respect to ψ results The solution of this equation gives the concentration profile where the interface width is ℓ = . The surface tension can then be calculated by integration of Eq. (27) where we keep only the terms involving the gradient of ψγ In order to look at the scaling variation with parameters of the mixture such as α T or α v , one must express the dimensionless quantities that we used as functions of these parameters. Taking as an example equalsized hard spheres where α v = 1, and hence α ζ = 1, ǫ A = ǫ B = β B = 8, the volume fractions at the critical point are given by φ * when α T ≫ 1. Considering a particle mixture of given volume fractions of particles (which could be named as the laboratory conditions) φ 0 A and φ 0 B which satisfy where v 0 is the molecular volume for finite volume fractions such that φ 0 0γ . The surface tension therefore increases as a power law of the ratio between the temperatures of the two types of particles α T (see Appendix D for details).

B. Kinetics of droplet growth
We consider a mixture quenched in the two phase region which is therefore supersaturated. We focus on a spherical droplet of phase-b with radius R, growing inside the background phase which has a composition close to phase-a. The volume fractions outside the droplet are different than the concentrations in the a phase due to the supersaturation, i.e., φ A → φ a A +δφ A and φ B → φ a B +δφ B . Note that the concentrations inside the droplet are also slightly different from the concentrations of the b phase but we can ignore this difference here.
This is a multi-scale problem of equations with two well separated length scales: the width of the interface ℓ is much smaller than the droplet size R. At the scale of the small length ℓ, the system is still close to equilibrium with chemical potentialsμ ′ A andμ ′ B which are the chemical potentials slightly shifted from the phase equilibrium values and calculated outside the droplet on its surface. The steady-state equations, which give the particle concentration profiles are In order to solve the so-called inner problem at the length scale ℓ, we choose a length δ such that ℓ ≪ δ ≪ R. We multiply the first equation by ∂φ A /∂r and the second one by ∂φ B /∂r, add them up and then integrate across the interface over a region of size δ such that in both phases, ∂φ A /∂r = ∂φ B /∂r = 0 away from the interface. As a result, we obtain the Gibbs-Thomson relation [38] : where for each species, ∆φ ab α = φ b α − φ a α and δμ α =μ † α − µ ′ α . Note that outside the droplet, as discussed below, the volume fractions and the chemical potentials vary over the large length scale R and do not change over the length δ. As in the previous section, we now transform the volume fractions to normal coordinates and obtain: The non-critical variable η is identical in the two phases so that ∆η ab = 0. The small variation of the chemical potential δμ ψ can be obtained from Eq.(C7) as δμ ψ ≃ 2k ψ ψ 2 a δψ using Eq.(C10) where δψ = ψ(R) − ψ a is the small shift of the order parameter on the surface of the droplet from its equilibrium value in phase-a. Then, rearranging Eq.(38), we have: We now study the dynamics of the growing droplet by studying the outer problem, calculating the order parameter profile of the droplet material. As the non-critical variable η has the same value in the two phases and as we find numerically that its variation is very small, we will assume here that there is no flux associated to this variable. The problem then has a single conserved order parameter ψ.
Outside the droplet, the order parameter ψ follows a diffusion equation ∂ψ/∂t = D ψ ∇ 2 ψ. We discuss the value of the effective diffusion constant D ψ in Appendix E. The boundary conditions for this diffusion equation are the value δψ(R) given by the Eq. (39) and the value at infinity ψ ∞ which measures the supersaturation. The solution of the diffusion equation is The growth of the droplet is due to the radial flux where the supersaturation is defined as ∆ = δψ(∞)/∆ψ ab where d 0 =γ/(k ψ ψ 2 a ∆ψ 2 ab ) is a length of the order of the interfacial width ℓ. This gives the critical nucleation radius of the droplet R c = d 0 /∆. Droplets smaller than d 0 collapse whereas droplets larger than d 0 grow.
At the early stages of the phase separation just after the quench, there are few droplets and the droplets that are larger than the critical radius grow as R ∼ t 1/2 . At long times, the value of the supersaturation decreases with time and a much more detailed analysis is required, which has been made by Lifshitz and Slyozov [39]. The supersaturation decreases as ∆ ∼ d 0 /R and the average droplet radius increases as R ∼ t 1/3 . Plugging in the value of D ψ obtained in Appendix E gives the scaling of droplet growth with time, R ∼ (r G t) 1/3 in which the growth rate of mean droplet volume r G ∼ (φ 0 A α T ) 1/2 v 1/3 0 T B /ζ that is linearly proportional to the geometric mean of T A and T B .
Here, we obtain these power laws from our two temperature model recovering the properties of equilibrium systems in the dilute limit. On the other hand, at higher order density expansions even though the solution parameters (∆, d 0 , D ψ , r G ) change, there is no reason to expect a different power law behavior than R ∼ t 1/3 as long as the droplet material is transported by diffusion. Similar examples include one component active fluids which respect the 1/3 law [40,41].
V. SOME APPLICATIONS OF THE THEORY As discussed earlier, the intra-and inter-species interactions can be controlled either by modifying the interaction potentials between particles or by changing the activity difference (the temperature ratio). In the first case one alters the virial coefficients while in the latter case, one changes both the entropic part of the effective free energy by changing the temperatures T α but also the weight of the interactions through the cross temperatures T αβ as seen in Eq. (10b). In general, for short-range interactions, the second virial coefficients can be written as B αβ ≈ b αβ + a αβ /T αβ , where b αβ is the effective pair excluded volume and a αβ is the additional interaction part. Moreover, b αβ can be approximated by its hard-core value. As an example, for spherical particles, αβ and a αβ ≈ 4π ∞ d αβ u αβ (r)r 2 dr where d αβ is the distance between the centers of the particles at contact for an αβ pair. If the additional interaction is purely attractive, a αβ < 0 is negative whereas it is positive if the additional interaction is repulsive. For a mixture of particles with given hard core sizes, this part of the virial coefficient can be tuned by chemical modifications. In addition, the contrast in activity which can be adjusted by the energy input per constituent, provides an extra handle [6,42] for controlling the phase separation. In this section, we show some examples of applications of the theory, in various systems.

A. Colloidal hard spheres with different temperatures
Pure hard spheres interact with each other only through excluded volume interactions, which do not allow them to interpenetrate: a αβ = 0 and B αβ = b αβ . The dramatic influence of activity contrast towards demixing is clearly observed in mixtures of hard spheres with equal sizes and hence equal mobilities [12]. In Fig.2a we show the phase diagram for equal-sized hard spheres α v = 1, α ζ = 1, then ǫ A = ǫ B = β B = 8 with temperature ratio α T = 20. The volume fractions of the cold and hot particles are not equal at the critical point where the volume fraction of the cold particles is larger, indicating asymmetric phases. As exemplified in the phase diagram, starting from a mixture in the unstable region, the system phase separates into a solid-like close-packed B particles (phase-β, which is not quantitatively well described by our low density approximation) surrounded by an A gas (phase-α) as illustrated by top left inset.
Using our formalism, we can also investigate the phase behavior of hard-sphere mixtures with different sizes. As we mentioned earlier, exclusively for active systems where T A = T B , do the mobilities (or the friction coefficients ζ α ) come into play via the effective cross temperature T AB when ζ A = ζ B . For a given size ratio α v , β B = (1 + α 1/3 v ) 3 and using Stoke's law α ζ = α 1/3 v . In Fig.2b we plot the phase diagram for mixtures of hard- Moreover, in a) we mark two coexisting phases a and b (purple squares) which are strongly asymmetric. Starting from the initial well-mixed setup, the system phase separates into a solid-like close-packed B particles (phase-b) surrounded by an A gas (phase-a) as illustrated in top left inset. In b) the phase diagram appears to be more symmetric indicating both liquid-like coexisting phases though shifted slightly towards A-rich side. This results from having αT /αv = 20/27 1, the ratio which controls the symmetry of two phases (see Section III D) spheres where the hot particles A are larger with a volume ratio α v = 27 and a temperature ratio α T = 20. The evolution of the phase diagram upon increasing the size ratio can be appreciated by comparison to equal-size hard-spheres at the same temperature ratio (Fig.2a-b). A more symmetrical phase diagram is predicted by increasing the volume ratio to reach α v ≈ α T and we expect coexistence between two liquid-like phases (Fig.2b). A further increase of size ratio shifts the phase diagram towards the A-rich side as expected for mixtures of hard spheres with different sizes at the same temperature. Another way to observe the effect of tuning both size and temperature ratios is to evaluate the conditions required for the existence of a phase separation given in Section III E as a function of temperature and size ratios. We show this phase diagram in Fig. 3a which displays reenterances in both size ratio and temperature ratio axes.
Moreover, in dilute mixtures of active swimmers and passive particles, the run and tumble mechanism of the swimmers with propulsion speed V A and reorientation time τ r dictated by the time step between two tumbling events, can be considered at long times as entirely diffusive. Thus, these systems can be described by the present theory. If the reorientation time of swimmers is much smaller than the mean collision time, i.e., τ r ≪ τ c , then (T A − T B ) ≈ P cs τ r /d where T B is the background temperature, d is the space dimension and P A cs = V 2 A ζ A is the mean power required to drive chemotaxis [43] for an active particle. Accordingly, keeping τ r high will serve lower dissipation rate to maintain the same radial diffusivity. Now if the A and B constituents have the same volume v 0 , at later stages of clustering, Eq.(E5) suggests that the growth rate of mean cluster size r G ∼ (φ A P cs τ r ) 1/2 v This can be generalized similarly when there are two different types of swimmers and so on.

B. Interacting particles with different temperatures
Another interesting scenario appears when the interactions between particles enhance mixing while the activity contrast opposes and boosts demixing. To illustrate an example, consider a system of particles with equal strength of interactions which are repulsive for identical particles and attractive between different particles. This could be achieved for example by mixing hot and cold particles of opposite net electric charges in a medium with a finite screening length. Similar other systems can be prepared by engineeering the chemical interactions. The interaction part of the virial coefficients can in this case be written as a AA = a BB = a 0 and a AB = −a 0 . We may further simplify the problem by considering spherical particles of equal sizes such that b αβ = b 0 = 8v 0 where v 0 is the molecular volume of the particles. We can evaluate the parameter range where a phase separation occurs. In Fig. 3b, we show the phase diagram in terms of the scaled interaction parameter a 0 /(T A v 0 ), and the temperature ratio α T = T A /T B for T A > T B . When a 0 > 0, the molecular interactions only promote mixing while at sufficiently large temperature ratios, the system can still reach phase separation. Interestingly, the phase diagram shows reentrance in both the lateral and vertical directions. In b) we use the same condition (in purple) for the existence of demixing for equal size hard spheres with short-range interactions of equal magnitude but opposing behavior (attractive vs. repulsive) for intra-and inter-species. The inset shows the favored interactions in both regimes. If a0 > 0, the molecular interactions only promote mixing while at sufficiently large temperature ratios, the system can still reach phase separation.

C. Active-passive polymer blends
Biopolymers play a key role in intra-cellular or intranuclear organization in many instances [44], while they often interact with active proteins which confer them an active character. In our context, such activity in biopolymers has been shown to enhance spatial segregation and maintain compaction. This is the case for example displayed in the structure and compaction of DNA inside the cell nucleus [45,46]. Similarly to colloidal particles, the active forces induce an effective temperature higher than the ambient one [47]. To give an example within our theory, we consider here a mixture of poly-A and poly-B chains in solution with equal lengths but with different temperatures (or activities) T A > T B and having only excluded-volume interactions. For dilute solutions, one can still expect an effective thermodynamic behavior with an effective free energy given by a Flory-Huggins theory. In the spirit of the Flory-Huggins meanfield theory, we suppose here that the interaction part of the free energy does not depend on the connectivity between the monomers and that we can employ the results obtained for colloidal particles. Therefore, in order to describe the polymers, we use the chemical potentials obtained in the section II A by making the transformation µ α → µ id α /N + Φ α , where N is the number of monomers in each chain. Here, we consider the size interactions to be identical such that ǫ A = ǫ B = β B = ǫ 0 . By using the phase separation conditions described in Section III (both approaches agree in this case when N ≫ 1), we observe that segregation requires (T A − T B )/T B = α T − 1 > (4ǫ −1/2 0 )N −1/2 . This result agrees with extensive simulations of active-passive polymer mixtures [15] where the same scaling law is observed, though in terms of effective temperatures.

A. General form
Up to this point, we have studied the general phase behavior of a suspension of mixed particles with different temperatures in the dilute limit. In this limit, the theory takes into account only two-particle correlations, the system has an effective thermodynamic description and the phase behavior can be obtained from the direct analog of the equilibrium construction of phase separation, despite the existence of non-equilibrium aspects such as the violation of detailed balance that are observed at the microscopic level [11]. A natural extension of this approach is then to calculate higher order corrections in concentration and see whether the effective thermodynamic description is preserved.
In order to answer this question, we must first obtain the steady-state pair distribution function g ss αβ (r 1 −r 2 ) at next orders in particle densities c α . A general strategy would be to start from the Fokker-Planck equation (2) for the multi-particle probability distribution P , and integrate up to the desired order in densities. One can then solve the remaining coupled equations to obtain the pair distribution functions g ss αβ (r 1 − r 2 ). This approach leads to the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [31]. For our problem, a nonequilibrium analog of this hierarchy is detailed in Ref. [11]. In the equilibrium case, the solutions are found by choosing a closure relation and simplified by setting the fluxes along each direction to zero. By contrast, in a non-equilibrium case, for T A = T B , there might exist non-vanishing fluxes associated to dissipation in the system [48]. This complication makes it difficult to obtain a systematic expansion at higher orders in densities.
A solvable example of Fokker-Planck equation (or the Langevin dynamics) at higher order has been given for pairwise harmonic potential interactions. In this case, a steady-state solution exists [49] for m ∂ m J m = 0 while the fluxes J m = 0 when T A = T B . As a result, it is not possible to formulate a solution in Boltzmann form with a scalar temperature.
Other classical approaches such as the Kirkwood superposition approximation would also fail. We refer the reader to the probabilistic interpretation of the Kirkwood closure in Ref. [50]. Nevertheless, in the following section, we demonstrate an alternative approach based on the calculation of depletion forces, to obtain third-order density corrections for pure hard-sphere interactions.

B. Hard spheres
In the case of mixtures with only hard-sphere interactions, an appropriate approach to expand at least to the next order, considers the depletion interaction between two particles due to a third particle [51]. This method has been used repeatedly in colloid science [52] and recovers exactly the third virial coefficients in hard-sphere mixtures with different radii. The full details of the calculation are given in Appendix F. Here as an example, we briefly sketch the results obtained for equal-sized hardspheres with different temperatures. The method is also applicable to particles of differing size ratios as shown in the Appendix F. The resulting pressure is given by the standard virial expansion: (42) where B and C are the second and third virial coefficients respectively, which are identical for all types of pairs and triplets of particles and α, β, γ = A or B.
In equilibrium systems, there is a direct route from pressure to chemical potentials using the Gibbs-Duhem where the chemical potentials are given by µ 0 α = ∂f 0 /∂c α . Our analysis for systems with two temperatures has shown that this remains applicable in the dilute limit approximation. However, without knowing explicitly the effective free energy, this approach is not founded. Thus, one should always start from the dynamic equations for the concentrations. If we insist in rewriting the interaction part of the chemical potentials defined in Eq. (8) as a series expansion in densities Φ α = Φ (1) α + Φ (2) α , we only obtain the value of the chemical potential gradient at second order in the densities: where C is the third virial coefficient between hard spheres of the same size (see Appendix F), which is equal for all types of interactions. At this order, ∇Φ A because the mixture contains particles of equal sizes. The total chemical potentials µ α are then obtained from Eqs. (9), (43). However, the gradient ∇Φ (2) α given by Eq. (43) is non-integrable to a potential function due to the mismatch between the cross terms. This incompatibility is related to the non-equilibrium character of solutions of particles at different temperatures and hence the mismatch vanishes for an equilibrium solution, when T A = T B . As a result, the non-equilibrium additional terms break down the routes to construct an effective thermodynamic theory. A similar breakdown has been previously observed in one component active fluids with density-dependent motilities [40] where the interfacial contributions lead to terms not integrable to a free energy. That issue has then been addressed in Ref. [53], by introducing a functional transformation using an alternative scalar order parameter. By contrast, here we are not able to define even the local chemical potentials at third order in the dynamical equations (6), (7). This could as well indicate the existence of bubbly phases [54] in which two steady state phases exist only locally and are separated by interfaces. A more detailed study would require to extend the analysis in gradient terms presented in the previous sections. We plan to address this question in a further study.
Since the term ∇Φ A is not integrable to a chemical potential form, it is not straightforward to determine all uniformly conserved quantities at the zero-flux steadystate. At phase equilibrium, we only find two conserved quantities since ∇p 0 = 0 and ∇(µ 0 A −µ 0 B ) = 0. The latter condition exists specifically for equal-sized hard-spheres, because ∇Φ A . A third condition however, can be obtained by linearization of one of the concentration fluxes around the critical point. Then, the calculated phase diagrams are shown in Fig.4 in comparison to the phase diagram obtained from the previous expansion using second order virial coefficients (dilute limit) for the same set of parameters α v = 1, α T = 40. Curiously, the contribution from the third order terms delays the onset of demixing.

VII. CONCLUSION
To summarize, in this work, we have outlined the general framework to study the characteristics of phase equilibria in active mixtures of diffusive type where the temperatures of the constituents are different. We obtain the phase diagrams and phase growth properties using two methods in parallel which cover both the equilibrium thermodynamic description through functional analysis and a more general approach considering steady-state solutions of the Fokker-Planck equations. For each observable and method, we compare and contrast the nonequilibrium T A = T B scenarios to the equilibrium ones T A = T B . This allows us to reexplore landmark methods and concepts developed to establish the foundational principles of equilibrium solution theory. Resemblances appear in steady-state solutions while this state is only maintained with a net power dissipation when T A = T B . Even though these systems carry a non-equilibrium profile, it turns out that they recover a direct analog of the equilibrium construction in the dilute limit. Thus, we were able to construct a Cahn-Hilliard theory in dilute limit now generalized to two temperature mixtures.
Yet, this analog of the equilibrium picture provides a rich palette to explore various aspects of active/passive (or less active) mixtures. As we demonstrate by examples, the theory has broad applications in diverse systems at different length scales. In this simplified formalism, we capture interesting scaling laws for interfacial properties, droplet growth dynamics, and for the phase segregation condition, some of which being equally observed in simulations. Our results also suggests a means to calibrate the composition of coexisting phases (liquid-liquid vs. gaslike and solid-like) by controlling the ratio α T /α v which could motivate experimental applications.
Higher order corrections (though not available for the general case) break down the direct analogy with the equilibrium construction in the case of pure hard-core interactions. However, the results do not indicate significant qualitative differences with the equilibrium behavior, and hence in general, the qualitative behavior of the system should be obtained from the simpler theory describing the dilute limit (Sections II-IV) that gives an intuitive understanding of demixing in diffusive systems. On the other hand, the existence of non-local terms in the chemical potentials at the third order in a power expansion in densities might indicate the emergence of new phenomena. It would be interesting to study these cases further including inhomogeneous terms, in order to bridge the gap between microscopic and coarse-grained models [53,54]. with flux components: Setting the flux J r = 0 for the steady-state solution results Eq.(4) in the main text. The remaining part only requires a uniform G(R) which satisfies the steadystate solution though J αβ R does not necessarily vanish for T α = T β .

Appendix B: Irving-Kirkwood method for calculation of internal stress
Following the stress equation for the interaction part, Eq.(18) in the main text up to second order in separation vector r while using Eq.(4) and keeping only nonvanishing terms upon integration, we rewrite the internal stress σ (v) ij given by Eq. (17) as: where p 0 = µ 0 A c A + µ 0 B c B − f 0 is the locally uniform pressure. Integrating by parts gives: where a summation on the k, l indices is performed using an Einstein summation convention, and the integral I ijkl is given by: Considering the symmetries and performing the sum, we finally reach an expression for the difference between stress obtained by two different methods: where σ ij is the result obtained by free energy deformation given in Eq. (15). Summing over α, β suggests that the addition of the gauge term discussed in Section II B, i.e., − 1 to the original free energy f , will recover Irving-Kirkwood formula.
Appendix C: Transformation of order parameters to normal coordinates around the critical point Close to the critical point {φ * A , φ * B }, the effective thermodynamic description becomes simpler if instead of using the volume fractions measured with respect to the volume fractions at the critical point δφ α = φ α − φ * α as variables, we make a linear transformation to the eigenvectors of the inverse compressibility matrix κ −1 p at the critical point. The inverse compressibility matrix is a symmetric matrix given by Eq. (23) and it can be written as The values of the coefficients are given by Eq.
αT +α ζ 1+α ζ β B . We denote the two eigenvalues of the matrix by ǫ and λ. The coefficients of the matrix are real and related to the eigenvalues by c = λ+ǫ 2 and (a 2 + b 2 ) 1/2 = λ−ǫ 2 . In the vicinity of the critical point, ǫ is small and vanishes as one approaches the critical point and the other eigenvalue λ remains finite at the critical point. The eigenvector associated with ǫ = 0 at the critical point, gives the direction in which the fluctuations diverge. We also define an angle θ such that a = − λ−ǫ 2 cos 2θ and b = λ−ǫ 2 sin 2θ. The diagonalization appears then as a rotation of angle θ and the eigenvalue matrix D ǫ,λ with diagonal entries ǫ and λ is obtained by rotation to the basis of eigenvectors D ǫ,λ = R(θ)κ −1 p R T (θ) using the rotation matrix R(θ).
In the volume fraction space, the eigenvectors are obtained using a rotation of the natural coordinates by an angle θ. The normal coordinates that we use are the coordinates along the eigenvectors of the inverse compressibility matrix at the critical point. At the critical point, the rotation angle of the eigenvectors is θ * that satisfies which can also be expressed in conjugate form in terms of φ * A using tan θ * = c * +a * b * . This angle gives the orientation of the tie lines close to the critical point. As a result, it provides an indication on the asymmetry of composition between the two phases. Thus, when tan θ * ≈ 1 , we would have two liquid-like phases, whereas for tan θ * ≫ 1, the critical point is towards B-rich side of the phase diagram with a solid-like and a gas-like phase coexisting and vice versa for tan θ * ≪ 1.
In the eigenbasis of the inverse compressibility matrix at the critical point, there is a coordinate ψ along the eigenvector associated to the vanishing eigenvalue and a coordinate η along the eigenvector associated to the finite eigenvalue λ. These two normal coordinates are related to the original volume fractions by the rotation matrix R(θ * ). Accordingly, we define: The coordinate ψ is the critical variable that we call the order parameter. Along these new coordinates, differentiation is performed as: (C4) Using these relations, we expand the free energy around the critical point and obtain: where the coefficients of the expansion are the derivatives of the free energyf 0 evaluated at the critical point: and we used the fact that η is a slowly changing variable. The Hessian matrix of the second derivatives of the free energy is equal to the inverse compressibility matrix. In the coordinates ψ, η the inverse compressibility matrix is D ǫ,λ evaluated at the critical point where ǫ = 0. The two derivatives The coordinates of the equilibrium phases a and b (the binodal line) are obtained by equating the chemical potentialsμ ψ andμ η and the pressure in the two phases. This leads to η a = η b and ψ 2 a = ψ 2 b = − k2 k4 η a . By construction of the normal coordinates, the tie lines, which are the straight lines joining the two phases at equilibrium correspond to lines of constant values of η.
The binodal line has therefore a parabolic shape with ψ a ≈ −ψ b . Each value of η a ≈ η b > 0 defines the coexisting phases ψ a and ψ b . It is convenient in the following to consider a symmetrized version of the order parameter ∆ψ ab = ψ b − ψ a , which satisfies ∆ψ 2 ab = − 4k2 k4 η a . The total free energy density is obtained by including the gradient terms which also can be transformed to the normal coordinates. In the vicinity of the critical point, as the non-critical variable η is much smaller than the order parameter ψ we only need to retain terms in (∇ψ) 2 . The total free energy density reads then where the coefficientL ψ is given byL ψ =L A cos 2 θ * − 2L AB cos θ * sin θ * +L B sin 2 θ * . In order to calculate the interfacial tension, we must define the two phases at equilibrium i.e. fix the values η a in the two phases in equilibrium. This also fixes the order parameter ∆ψ ab . We must then calculate what we called the tilted free energy around the critical point The constantsμ † ψ andμ † η are the chemical potentials calculated in the two phases at equilibrium. Note that the tilted free energy has a minimum and vanishes in the two phases in equilibrium. The profiles of ψ and η as a function of the coordinate z perpendicular to the interface are obtained by minimization of the tilted free energy. We first minimize the tilted free energy with respect to η. This leads to Inserting this result into the tilted free energy, we obtain the tilted free energy as a function of the order parameter ψ only, which we write as 2λ * . This free energy is minimal and vanishes at ψ = ψ a and ψ = ψ b . Finally, minimization with respect to ψ results: which shows that along the order parameter profile between the two phases, ∆f (ψ) = 2∆f 0 (ψ).
Appendix D: Scaling relations for equal-sized hard spheres when αT ≫ 1 As shown in the main text for equal-sized hard spheres where α v = 1, (and hence α ζ = 1), and ǫ A = ǫ B = β B = 8, the volume fractions at the critical point are φ * A = α −1 T , φ * B = 1/8 + (5/4)α −1 T when α T ≫ 1. Using these relations, the rotation angle at the critical point given by Eq. (C2) reads, We consider a mixture with average particle volume frac- We determine the normal distance η(φ 0 A , φ 0 B ) from the critical point. Then, as discussed in Appendix C, this fixes the order parameter ∆ψ 2 ab = − 4k2 k4 η(φ 0 A , φ 0 B ), and hence the compositions of the coexisting phases. In the limit where α T ≫ 1, the order parameter is: which is valid at (reasonably) finite volume fractions such that φ 0 A α T ≫ 1 where the first term dominates. In order to discuss the scaling of the surface tension, we also need the expressions of the coefficients k ψ and L αβ as functions of α T . For equal-size hard spheres, we have:L where σ D = 4 5 (6/π) 2/3 and v 0 is the molecular volume. We can then usz these results to obtainL ψ = L A cos 2 θ * − 2L AB cos θ * sin θ * +L B sin 2 θ * . The other coefficient k ψ = k 4 − k 2 2 2λ * can be calculated by taking into account the definitions below Eq.(C5) and differentiating Eq.(C4). For α T ≫ 1, k 4 ≈ 256 sin 4 θ * , k 2 ≈ −16α T sin 3 θ * and λ * ≈ α 2 T sin 2 θ * . Accordingly, k ψ ≈ 128 in this limit which is independent of the temperature ratio similarly toL ψ . These relations using (33) lead to Eq.(34) in the main text.
Appendix E: Calculation of D ψ and the droplet growth rate In order to determine the diffusion coefficient of the relevant order parameter ψ given in IV B, we need first to linearize time evolution equations (7) around one of the phases (say phase-a). This suggests: in which we define δφ α = φ α − φ * α and obtain Γ a by evaluating (23) at phase-a. Then, following Appendix B, we obtain: where Γ ′a = R(θ * )Γ a R T (θ * ). Using ∂η/∂t ≈ 0, we have D ψ = |Γ ′a |/Γ ′a 22 with |Γ ′a | = |Γ a | being the determinant of matrix Γ ′a and Γ ′a 22 is the second diagonal component. Around the critical point, using Eqs.(C4),(C7) and (C10), the diffusion coefficient is: where φ a A and φ a A are the values of volume fractions at phase-a. At high temperature ratios α T ≫ 1, we see that sin θ * ≈ 1 and hence: where we consider equal-size hard-spheres, α v = α ζ = 1 and ζ A = ζ B = ζ. Next, we estimate the growth rate of average droplet size r G such that R ∼ (r G t) 1/3 . Near saturation ∆ ∼ d 0 /R and hence r G ∼ D ψ d 0 . Using the above relation and the value of d 0 from the main text, we obtain r G ∼ T B ∆ψ ab v Appendix F: Third order expansion for pure hard-spheres

Virial expansion to third order
The strategy that we follow in this appendix is to define a potential of mean force at a finite concentration u tot αβ (r) for an αβ pair and to insert it into the Fokker-Planck equation (A1), (A2) by replacing u αβ → u tot αβ . This results in the steady-state value of the pair distribution function written in the Boltzmann form g ss αβ (r) = exp −u tot αβ (r)/T αβ where g αβ is the pair distribution function, u tot αβ (r) is the potential of mean force, and T αβ is the cross temperature between the two particles under consideration. In general, we can separate the potential of mean force into two parts, u tot αβ (r) = u αβ (r) + W αβ (r) where u αβ (r) is the bare pair interaction potential for two particles while W αβ (r) [55] which is due to the existence of other surrounding particles. The potential of mean force can be derived by geometrical considerations for hard-sphere mixtures.
When the excluded volumes between the two particles overlap such that a third particle γ cannot enter in the space between the two particles, there is a net attractive force between the two particles known as the depletion force [51]. The mean force between an αβ pair due to a third particle is given as F αβ = −p tot 0 S where p tot 0 is the total pressure of third body particles and S(r) is the cross-sectional area loss at the overlapping region, which depends on the distance r between the two particles of an αβ pair. For equal-size particles with diameter d, this surface is given by: The integration of this force F αβ yields the potential of mean force and the interaction W αβ (r) between the αβ pair mediated by a third particle within the region d ≤ r ≤ 2d. By imposing continuity of the interaction potential, we obtain: where w(r) is the overlap volume and Θ(x) is Heaviside step function. Note that since p tot 0 is a function of the third particle only, the mean force and the corresponding potential are independent of the αβ pair, i.e., W αβ (r) = W (r). In a more formal manner, this potential can be expressed as: Then by using the piecewise definition of the hard-sphere potential we write, up to first order in concentrations: where the pairwise hard-sphere potential u αβ (r) = ∞ for r < d and u αβ (r) = 0 elsewhere. We now use this pair distribution function to obtain thermodynamic properties. The pressure is obtained from the virial equation which becomes in the case of hard-spheres of equal-sizes: Inserting the expression of the pair distribution function Eq. (F4), we obtain the pressure as: The second and third virial coefficients are therefore B = 2π 3 d 3 and C = 2π 3 d 3 w(d) with w(d) = 5π 12 d 3 from Eq.(F2). Rewriting in terms of the virial coefficients gives Eq. (42) in main text. Note that the third order terms in the expansion are weighted by the temperature of the third particle. Thus, introducing different mobilities will not alter the form of this equation.
The next task is to implement the same strategy to the time evolution equations for the concentrations c A and c B , given by Eqs. (6), (7). In a similar way, we define the total mean forcef α on particles m, as the gradient of a potentialf α = −∇Φ α . Using Eq.(F4), the interaction part of the chemical potentials is expanded in powers of concentrations as Φ A = Φ (1) A + Φ (2) A . The lowest order term Φ (1) A has been obtained in terms of the second virial coefficients in Eq. (9). We obtain the next order contribution using Eqs.(F3) and (F4): T γ c γ (r 3 )Θ(d − r 13 )Θ(d − r 23 )dr 3 dr 2 .
By coordinate transformation, we get: where x denotes the unit vector along the direction of a vector x. The second integral is constrained over the overlap volume V ∩ . Taking into account the Dirac delta function δ(r − d), it can be expressed as: where y ′ = cos θ ′ . Finally, by Taylor expanding the concentrations around r 1 up to first order in displacement vectors, we obtain: V∩ c β (r 1 )c γ (r 1 ) + c γ (r 1 ) r · ∇ r1 c β (r 1 ) + c β (r 1 ) r ′ · ∇ r1 c γ (r 1 ) dr ′ dr .
For equilibrium systems with T A = T B , this result gives back the third-virial coefficients between hard-spheres. In the general case, for mixtures of particles with two different temperatures, we obtain Eq.(43) of the main text and ∇Φ (2) B = ∇Φ (2) A . Hence, it appears that this field is not integrable to obtain Φ (2) B . On the other hand, we observe that c A ∇µ A + c B ∇µ B = ∇p consistent with the result obtained from virial equation (F6).
This method could as well be implemented for mixtures of hard-spheres with different diameters d A = d B as well as different temperatures T A = T B . This would require to solve the general form of Eq.(F8) with varying contact distances, T γ c γ (r 1 + r ′ )Θ(d 13 − r ′ )Θ(d 23 − |r − r ′ |)dr ′ dr.

Phase separation and coexistence conditions
Following the procedure described in Section 3, we determine the phase diagrams at third-order in concentration. A first remark is that we find an instability when T A /T B 6.751 for equal-sized spherical particles. Hence, the third order terms delay the onset of phase separation compared to the second order expansion which leads to phase separation when T A /T B > 4. The phase coexistence conditions can be calculated by imposing that the particle fluxes to zero as well as the momentum flux. This last condition is obtained from c A ∇µ A + c B ∇µ B = ∇p = 0, which gives mechanical equilibrium, and therefore imposes that pressure is constant. The other conditions require that ∇µ A = 0 and ∇µ B = 0. However, the chemical potential gradients are not integrable. Nevertheless, since ∇Φ B , we can use the condition that ∇(µ A − µ B ) = 0 which would lead to two gradient free equations. Therefore, to summarize, we find two conditions for the the a and b phases to coexist where we defined with µ (1) α , α = A, B. The third condition required to construct the phase diagrams is not directly accessible in a general form due to the lack of well-defined chemical potentials. Yet, it is possible to derive an approximate condition in the vicinity of the critical point by linearizing the time evolution equations in concentrations.