Thermodynamic stability and critical points in multicomponent mixtures with structured interactions

Theoretical work has shed light on the phase behavior of idealized mixtures of many components with random interactions. However, typical mixtures interact through particular physical features, leading to a structured, nonrandom interaction matrix of lower rank. Here, we develop a theoretical framework for such mixtures and derive mean-field conditions for thermodynamic stability and critical behavior. Irrespective of the number of components and features, this framework allows for a generally lower-dimensional representation in the space of features and proposes a principled way to coarse-grain multicomponent mixtures as binary mixtures. Moreover, it suggests a way to systematically characterize different series of critical points and their codimensions in mean-field. Since every pairwise interaction matrix can be expressed in terms of features, our work is applicable to a broad class of mean-field models.


Main text
Determining the phase behavior of mixtures is an important goal of Statistical Physics. But while the thermodynamics of mixtures with few components is well understood theoretically [1], most functional mixtures are made up from a large number of distinct components, and the principles underlying the phase and critical behavior of such multicomponent mixtures are less clear. There have been substantial steps towards understanding these systems, but only in limiting cases. Sear and Cuesta [2] and subsequent followups [3,4] determined conditions for phase separation in idealized mixtures with random, independent pairwise interactions. Taking a very different limit, Sollich and others [5][6][7] have made progress for polydisperse mixtures interacting through a continuous distribution of attributes. Different from these theoretical studies, many physical examples are made up of defined components whose interaction structure is governed by the physical details that underpin them. In phase-separation prone lipid membranes, though there are thousands of chemical species, interactions are thought to be primarily driven by just a few features -interactions between headgroups, the degree of acyl-chain saturation, and the mismatch between hydrophobic height [8][9][10][11][12]. In protein condensates, interactions are likely mediated by a combination of specific motifs like repetitive binding domains, and less specific electrostatic and hydrophobic interactions [13][14][15]. This observation suggests that in both cases the resulting effective pairwise interaction matrix is non-random in a particular way: its rank, given by the number of independent features, can be considerably smaller than its dimension, given by the number of components. In other examples it might be less clear what features mediate interactions, but an effectively low-rank interaction matrix is likely common to most mixtures made up of many components. One class of examples are fluids like petroleum for which an approximation in terms of lower-dimensional interaction parameters has been successfully applied [16]. To systematically investigate the role of such a low-dimensional interaction structure for phase behavior, in this work we develop a theoretical framework to study the phase behavior of mixtures with many components but structured interactions. We show that the stability of phases and critical behavior can be understood in a "feature space", which is typically much lower-dimensional than the space of component densities. Mean-field model. We specifically consider a family of multicomponent models with a pairwise interaction matrix of variable rank (see Fig. 1). The mixture is made up of N different component types. Component type i is characterized by a "feature vector" s i composed of R real features s J (2) J (3) J (4) s ( where xy is the sum over all neighboring sites x, y on the lattice and the spins take the values σ (α) (x)=s (α) i ∈R if site x is occupied by component type i (see [17] for a related model with a single feature). In a mean-field approximation, our system is described by a Flory-Huggins-like free energy density per k B T [18,19] where the densities ρ i are subject to the incompressibility constraint i ρ i =1. The interaction matrix is given by for a lattice coordination number z. For R features, the interaction matrix is of rank r≤R. While this decomposition into features may be motivated by the physics of interactions, any real and symmetric interaction matrix χ ij can be decomposed in this way, with R≤N , eigenvectors s (α) i and eigenvalues C (α) |s (α) | 2 ; see [16] for a related (eigen)decomposition in the context of petroleum. A (precise) way to think about the features is thus as eigenvectors of the interaction matrix. Furthermore, as long as interactions are pairwise and meaningfully described by mean-field, our choice of representing components in terms of additive Ising-like features is entirely general. For now, we assume all eigenvalues to be positive, but discuss the general case in the SM and briefly below. Thermodynamic stability and critical points. The main challenge in working with mixtures with N 1 components is that they are embedded in a very high-dimensional space of densities. We now develop an analytic framework, wherein the mixtures are instead represented in the corresponding, potentially much lower-dimensional feature space. To this end, we use matrix inversion techniques and will then successively derive conditions for local thermodynamic stability and the occurrence of different series of critical points. In general, the thermodynamic behavior of the multicomponent mixture is determined by the free-energy landscape in the N -dimensional space of densities ρ (N ) . Due to the incompressibility constraint, the densities are not independent, ρ N =1− N −1 i=1 ρ i , leaving the free energy density f a function of N −1 densities and temperature. The mixture is (locally) particular feature α.
thermodynamically stable if the Hessian matrix (i, j=1,. . .,N −1) is positive definite. Here r is the rescaled and shifted feature vector. While in the presence of a solvent, component N is most straightforwardly associated with this solvent, all physical results are ultimately independent of the choice of reference point. At high temperatures (large T → small C (α) → small r (α) ), the system is dominated by entropy (H≈K) and all eigenvalues of the Hessian matrix are positive; the system is thermodynamically stable to perturbations which are local in composition. At lower temperatures one or more of the eigenvalues can become negative, implying that the system can spontaneously lower its free energy by phase separating along the corresponding eigenvector. The boundary of local thermodynamic stability is called the spinodal. It corresponds to the submanifold of compositions and temperature where the smallest eigenvalue of the Hessian matrix is zero. Since a matrix is invertible if and only if all of its eigenvalues are non-zero, the spinodal is also the submanifold where the matrix becomes singular for the first time, starting from all positive eigenvalues. Here we take advantage of the fact that H is the sum of a positive definite matrix K with inverse K −1 ij = δ ij ρ i − ρ i ρ j and a lower rank contribution U U T arising from interactions (Eq. 2). Analogously to Ref. [20], we use the Woodbury matrix identity [21] on K and U U T to invert H and find that the Hessian matrix is invertible if and only if 1 − U T K −1 U =: 1 − Cov is invertible; see SM. Here, Cov is the covariance matrix between the (rescaled) features: where the averages are taken with respect to the probability measure given by the mixture composition ρ (N ) : X The rank of the covariance matrix R Cov corresponds to the maximal number of linearly independent feature vectors, R Cov ≤ R; see SM. Thermodynamic stability These results imply that the mixture becomes unstable when the largest eigenvalue of the covariance matrix λ (1) =1. Importantly, this condition is independent of the number of component types, including as limits the two-component mean-field Ising model and the infinite-component limit as discussed in the context of polydisperse systems [5,22]. Direction of instability In order to find the initial direction of phase separation at the spinodal we next determine the eigenvector corresponding to eigenvalue 0 (1) of the Hessian (covariance) matrix: If 1−Cov is invertible, the inverse of the Hessian matrix is given by H −1 = K −1 +K −1 U (1−Cov) −1 U T K −1 . Using the eigendecomposition of the covariance matrix in terms of its (descending) eigenvalues λ (γ) , γ = 1, ..., R and corresponding orthonormal eigenvectors V (γ) (whose dependency on ρ (N ) we drop for conciseness) Cov αβ = R γ=1 λ (γ) V (γ) α V (γ) β , the inverse of the Hessian is Close to the spinodal (λ (1) ≈1), the dominant term is 1 1−λ (1) e (1) e (1) T . Correspondingly, on the spinodal e (1) is the eigenvector of the Hessian with eigenvalue 0 2 and coincides with the direction of instability He (1) = 0 (see also Ref. [5] for polydisperse systems). Equation 4 implies that the relative enrichment δρ i /ρ i ∼ e (1) i /ρ i = E (1) i of the feature vector (relative to the mean) onto the 1st PC of Cov, see Eq. 4. Coarse-graining the multicomponent mixture as an effective binary mixture that preserves the location of the system with respect to the spinodal and critical manifold (by conserving the second and third cumulant along the 1st PC of Cov) leads to a composition as shown by the large, translucent disks; see SM for details. For a critical mixture (A), the composition of the binary mixture is symmetric, for a non-critical one (B), the densities of the two components are different. Note that all results are independent of global rotations or translations of the feature vectors or reflections s i ∀i; see SM. Right: For the multicomponent mixture to be critical, the skewness of the distribution of relative enrichments has to be zero (mixture 1 in A1 vs. mixture 2 in B1), see Eq. 5. make contact are critical points (cp). At a usual critical point ρ (cp) , two phases become indistinguishable, corresponding to two minima and one maximum of the tilted Landau free energy f →f − i ρ i ∂ i f | cp merging into one minimum. This merging occurs when the change in free energy along the direction of instability δf =f (ρ (cp) + e (1) )−f (ρ (cp) ) is zero up to order ≤ 3 in . The first order term of the tilted free energy is zero by definition and (∂ i ∂ j f )e (1) i cp =0 as ρ (cp) lies on the spinodal, yielding the following additional condition for the critical point: where the average is with respect to the density at the critical point ρ (cp) ; see SM (compare also [5]). Thus, the third cumulant (skewness) of the partition coefficient needs to be zero at a critical point; see Fig. 2. This condition on the third cumulant extends and substantiates the notion that binary systems and systems composed of ideal random copolymers are critical if their mixture composition is symmetric [23,24]. Notably, the conditions themselves (λ (1) =1 on spinodal, skewness =0 at a critical point) are valid irrespective of the mixture composition or feature distribution. To illustrate this generality, for Fig. 2 we randomly generated the components' features (R=2) either (A) via a multivariate Gaussian with zero mean or (B) as two independent features following a Poisson distribution (plus Gaussian noise) with non-zero mean and a Gaussian, respectively. In both cases, the mixture composition is drawn from a uniform distribution over the N −1-simplex; see SM. This procedure results in interaction matrices of rank 2, while "usual" random matrices have full rank [25].
Higher-order critical points At an n-th order critical point ρ (cp) , n phases become indistinguishable. For a single order parameter (density) ρ, this condition corresponds to the merging of n minima and n−1 maxima into a single minimum of the tilted Landau free energy. The free energy expansion around the critical point is then of the order 2n: δf ∼O δρ 2n . In a high-dimensional density space, the phases that become indistinguishable when crossing the n-th order critical point ρ (cp) do not necessarily lie on a straight line. Instead, the phases merge along a more general smooth curve ρ i ( )=ρ (cp) i +δρ i ( ) in density space, parameterized by ; see also [5]: k =0 to conserve the incompressibility constraint. The (tilted) free energy change δf ( )=f ( ρ c +δ ρ( ))−f ( ρ c ), whose first order term vanishes, should be of order 2n in : the same time, the system has to be stable against fluctuations in orthogonal directions. Thus, we determine the curve δ ρ around ρ (cp) in a way that it minimizes the free energy change δf up to the respective order; see SM. Systematically minimizing and setting the coefficients in front of m to zero, we find the following conditions for an n-th order critical point in terms of the partial exponential Bell polynomials B m,l ; see SM: which only depend on the vectors Υ (m) =:ρ (cp) Ω (m) , m=1, ..., n−1 determined recursively: Ω (1) = E (1) and whereB m,l := (−1) l (l−1)!B m,l (Ω (1) , ..., Ω (m−l+1) ). In the case of a single feature, R=1, this recursion is solved by , and the conditions for an n-th order critical point reduce to κ (s) m is the m-th cumulant of the spin s with respect to ρ (cp) . Thus, the more cumulants (order m ≥ 3) of the spin distribution are zero, the more phases become indistinguishable and the higher the order of the critical point can be. Merging of phases necessarily happens along the single direction of instability and there is only one series of higher-order critical points -the one just discussed. This is not true for R>1, which we discuss next. Multiple directions of instability A series of critical points distinct from the previously discussed higher-order critical points occurs when the largest eigenvalue 1 of Cov is D-fold degenerate; see SM. To ensure stability along any direction in the corresponding D-dimensional subspace of eigenvectors, the third cumulant of all vectors in the subspace then needs to equal 0. For example, a system has a critical point with two unstable directions if it has a two-fold degenerate maximal eigenvalue, λ (1) =λ (2) =1, and if the four distinct third cumulants κ (αβγ) : The emergent symmetry of these critical points is reminiscent of the order-parameter symmetry in q-state Potts models, with q=D+1. In two dimensions there are known to be critical transitions in the q-state Potts model for q≤4 [27], but in mean field these transitions are first order [28,29] except for the case of the Ising model, q=2 [30,31]. Our results show that it may be possible to have critical transitions in mean field models that have the symmetry of q-state Potts models, but only in models with sufficient flexibility. For q=3, an example of an N =6 component model with R=2 features is given explicitly in the SM. Model extensions So far we have focused on a mean-field free energy whose interaction matrix has positive eigenvalues and whose components are of the same size. We now briefly discuss how the previously discussed conditions generalize to the case of negative interaction strengths, and comment on components with different sizes as in the original Flory-Huggins theory [18,19] in the SM. If the feature interaction strengths J (γ) are all positive, the pairwise interactions satisfy 2χ ij −(χ ii +χ jj )<0 ∀i, j, and interactions between alike components are always energetically preferred compared to dislike components. To resolve this limitation, we consider the general case with R + "positive", attractive features and R−R + "negative", repulsive ones: J (α) >0 ∀α=1, . . ., R + and J (α) <0 ∀α=R + +1, . . ., R. Performing a conceptually similar but more intricate analysis as before, the spinodal criterion is λ C is the largest eigenvalue of the real, symmetric matrixC=C (++) −C (+−) (1+C (−−) ) −1 C (−+) , which is determined by the covariances among the subsets of positive (+) and negative (−) features; see SM.C has dimensions R + ×R + and can be interpreted as representing a multicomponent system with R + positive features and effective, reduced interactions. The extent to which the negative features influence the phase behavior depends on the relative correlations between all features. If for each dominant positive feature, there is a highly correlated negative feature of similar strength, their effects will roughly cancel and the mixture will not phase separate. Conversely, if the dominant positive features driving phase separation correlate weakly with the negative features, thermodynamic stability is barely modified by the presence of the latter. At the spinodal, the direction of instability isē in terms of the first eigenvector φ (1) ofC and the deviations of positive (π) and negative (ν) features from the mean; see SM. We observe that this direction of instability again corresponds to (a combination of) feature deviations from the mean, projected onto the first principal component φ (1) (now of the "effective covariance matrix"C). Roughly speaking, the relative sign of the contributions of the negative and positive features depends on whether they are correlated or anti-correlated (negative or positive sign). Finally, performing the same analysis as for the original model, we find an analogous condition for the ordinary critical point: i . Discussion In this work, we consider a general mean-field model for multicomponent mixtures with an arbitrary pairwise interaction matrix χ ij of variable rank which we decompose in terms of different "features" mediating additive interactions between the components. The analytic conditions we derive for the spinodal and (higher-order) critical points only depend on the distribution of components in feature space. Specifically, the spinodal and submanifold of ordinary critical points are determined exclusively by the variance and third cumulant of the component distribution projected along the first principal component of the feature covariance matrix; Fig. 2. This representation in feature space is reminiscent of the dimensional reduction obtained for polydisperse systems whose excess free energy only depends on a few generalized moments of the attributes [5][6][7]. While the derivation of the "moment free energies" relies on either a division of density space into a subspace of moments and its "transverse" space or on combinatorial arguments [5,6], here we instead exploit that the condition for the Hessian matrix to become singular only depends on an R-dimensional matrix originating from the interaction structure. A related simplification of the spinodal condition in terms of a lower-dimensional matrix has been achieved for Flory-Huggins models with an excess free energy depending only on a finite number of moments of the molecular weight distribution [32][33][34]. The representation in feature space also suggests a principled method for finding coarse-grained binary mixtures with similar properties. By choosing the composition and interaction strength of the binary mixture so as to preserve the second and third cumulant along the first principal component, the coarse-grained binary mixture maintains the location of the multicomponent system with respect to the spinodal and critical manifold; see Fig. 2 and SM. In addition, our analysis allows for a systematic identification of the codimension of different series of critical points in multicomponent systems; see also [35,36]. For instance, we find that, in the absence of symmetries, a tricritical point has codimension four in mean-field. Furthermore, higher-order critical points with symmetry reminiscent of the q-state Potts model require tuning of q 2 + q+1 3 parameters. For the q=3-states Potts model, this counting suggests a codimension of seven for the critical point, which is larger than the one accessible with just N =3 components but feasible for a mean-field model with N =6 components and R=2; we explicitly construct such a 3-states-Potts-like model containing a critical point; see SM.
Our results offer an appealing avenue towards understanding intracellular liquid-liquid phase separation [15] and the critical phase behavior observed in cell-derived plasma membranes [37]. These mixtures are composed of thousands of proteins (and lipids), and depending on the conditions, small domains form spontaneously. The number of coexisting domains appears to be orders of magnitude smaller than the number of components and is thus well below the limit set by Gibb's phase rule [38]. In cell-derived plasma membranes, while true phase separation occurs when cooling them below the critical temperature [39], nanoscopic domains observed at physiological temperatures [40] have been suggested to be critical fluctuations close to a thermodynamic critical point in the 2D Ising universality class [37]. Strikingly, specific lipids and proteins robustly partition into specific phases -seemingly under fairly broad conditions [41]. Our work offers an interpretation of this experimental observation: phase behavior is determined by just a few important features. Looking for such a low-dimensional feature space representation might help to make sense of the growing amount of experimental data generated by proximity-labeling techniques [42], and should provide important insights into the physical characteristics underlying intracellular phase separation. In these biological systems, effects of finite dimension (two or three) and sequence-dependent interaction patterns [43, 44] will likely quantitatively, but not qualitatively change the mean field picture we present here. Finally, our analytic theory only makes predictions about local thermodynamic properties but cannot now make statements about the global phase behavior, which would require knowledge of the full free energy landscape [1,3,4,45,46]. Whether global phase behavior can be understood in feature space is an interesting question for future research.
[37] S. L. Veatch, P. Cicuta, P. In this Supplemental Material we present in detail the derivations of the conditions for the spinodal and critical points as outlined in the main text and explain how we generated the distribution of features in the multicomponent mixtures depicted in Fig. 2 of the main text. We first discuss the mean-field approximation of the lattice Hamiltonian in terms of a free energy of Flory-Huggins type. In this approximation, we then derive the condition for the spinodal and the intial direction of phase separation, the direction of instability. To find the conditions for (higher-order) critical points, we further systematically expand the mean-field free energy around the critical density. After deriving these conditions, we look at two extensions of our original model, systems with components of different sizes and mixtures with partially negative and partially positive interaction strengths, and derive the corresponding conditions for the spinodal and ordinary critical points. To show that our result do not depend on choosing the features linearly independently, we then comment on the case of linearly dependent features. In a similar spirit, we further show that the system is invariant under global rotations and translations of the rescaled feature vectors. While our focus is on the case where the largest eigenvalue of the covariance matrix is non-degenerate, we then dedicate one section to the case with several directions of instability. Based on this discussion, we further propose an effective mean-field model for the two-dimensional q=3-states Potts model which exhibits a critical point, in contrast to the regular mean-field theory of the Potts model. We then briefly discuss our idea how to coarse-grain multicomponent mixtures in terms of an effective binary mixture. Finally, we provide details about how we generated the multicomponent mixtures illustrated in Fig. 2 of the main text.
Note that a glossary containing the most important symbols and definitions can be found at the end of the SM.
where γ = 1, . . . , R denotes the different features, x and y label the sites on the lattice and xy denotes the set of all nearest neighbors on the lattice.
i . Since we assume that each lattice site is occupied by exactly one component (particle), the overall number of lattice sites equals the number of particles. In a mean-field approximation, we neglect correlations between the different sites on the lattice and the different sites are not distinguished anymore. The probability that component type i occupies any site x is ρ i . Denoting by z the coordination number of the lattice, there are z/2 pairs of nearest neighbors per particle. For each pair of nearest neighbors xy with component type i on x and j on y, the energetic contribution is j . The probability for this configuration is ρ i ρ j . Thus, the MF contribution per particle is (z/2) j . Overall, the energetic contribution per particle is where M denotes the total number of particles/lattice sites. Combining the energetic contribution with the entropic contribution, k B T N i=1 ρ i log ρ i , per particle, we find the following mean-field free energy per particle and k B T with The subscript N indicates thatf is a function of all N densities ρ (N ) .

S1.1 Flory-Huggins interaction matrix
Comparing the energetic contribution from the pairwise interactions with a Flory-Huggins term, − N i,j=1 ρ i χ ij ρ j , we identify the interaction matrix as (S5) S1.2 Effective reduction to an N −1-component system Due to the incompressibility constraint, not all the densities can be changed independently. Here, we directly impose this constraint by replacing in the free energy. Instead off N , we consider the "restricted" free energy which only depends on the first N −1 densities ρ (N −1) . Note that as expected and as we will see explicitly below, the "physical results" do not depend on this choice of integrating out component N (instead of any other or a combination of them).
S2 Derivatives of the free energy f N −1 and the Hessian matrix Taking the derivative of the restricted free energy f N −1 , Eq. S7, with respect to the density ρ i , i = 1, . . . , N − 1 yields For the Hessian matrix H ij , i, j = 1, . . . , N − 1, we find where we define 1 ij := 1 ∀i, j in order to keep track of sums over repeated indices.
corresponds to the Hessian of an incompressible N -component system without interactions. Furthermore, denotes the rescaled and shifted features. Similarly, the higher-order derivatives n ≥ 3 are

S3 Spinodal condition and direction of instability
For a non-interacting multicomponent mixture, which is entirely dominated by entropic effects and accordingly does not phase separate, the Hessian matrix equals K. All eigenvalues of the Hessian are thus positive: where in the last step, we used that the variance of v with respect to the probability measureρ i = ρ i /(1−P), i = 1, . . ., N −1 for the N −1-component system is always larger or equal to 0. Note that Introducing (positive) interactions among the component types, the system becomes more and more prone to phase separate and the eigenvalues of the Hessian decrease. The system becomes marginally thermodynamically stable if the smallest eigenvalue of the Hessian crosses 0. Since a square matrix is invertible if and only if all eigenvalues are non-zero, marginal stability is equivalent to the Hessian matrix becoming singular for the first time.

S3.1 Inverse of Hessian matrix
From Eq S9, we observe that the Hessian corresponds to a rank ≤ R correction of the matrix K. The existence of the inverse of such a rank correction to K and its explicit form has been determined by Woodbury [21]: In this case, the inverse is given by and W = U T , we find that the Hessian matrix, Eq S9, is invertible if and only if K and 1−U T K −1 U are invertible. As we will see next, K is always invertible and 1 − U T K −1 U is rewritten as 1 − Cov, with the covariance matrix as defined in the main text: Inverse of K K itself is written as a rank-1 correction to a diagonal matrix M : Since P ≤ 1, this condition is always satisfied and the inverse of K is Here and in the following, averages . . .
are with respect to the densities of the N -component system (for simplicity, we will often omit the N and the vector notation): (S18) Taken together, we find that the Hessian matrix is invertible if and only if 1 − Cov is invertible. Note that the Hessian matrix is an For positive interaction strengths, Cov is a true covariance matrix, which is real, symmetric and positive semi-definite. The eigenvalues of 1 − Cov are thus all real and ≤ 1. For high temperature or weak/no interactions, the rescaled features r (α) are small and thus 1 − Cov has only positive eigenvalues. Decreasing the temperature or increasing the attractive interactions, the mixture becomes marginally stable when 1 − Cov becomes singular for the first time, i.e. when the largest eigenvalue of Cov is 1. Writing the covariance matrix in terms of its eigenvectors V (γ) and eigenvalues λ (γ) , γ = 1, . . . , R: with λ (α) ≥ λ (β) ∀α ≤ β, this condition for marginal stability (spinodal) is

inverse of the Hessian is given by
or written in terms of the eigenvectors and eigenvalues of the covariance matrix, where e (γ) i and we used that Approaching the spinodal, the inverse of the Hessian is dominated by the term ∼ 1/(1 − λ (1) ) 3 , suggesting that on the spinodal e (1) is the eigenvector of the Hessian to eigenvalue 0. Using that r (γ) In particular, He (α) = 0 if λ (α) = 1. On the spinodal, where λ (1) = 1, e (1) thus corresponds to the eigenvector of the Hessian to eigenvalue 0:

S3.2 Degenerate maximal eigenvalue of Cov: Dimensionality of the plane of instability
If Cov has a D-fold degenerate maximal eigenvalue of 1, the corresponding eigenvectors can be chosen orthonormal since Cov is real and symmetric. Therefore, they span a D-dimensional submanifold in feature space. One may wonder if also the corresponding plane in component space, spanned by the different directions of instability e (1) , . . . , e (D) , is D-dimensional, or if these directions of instability might be linearly dependent. Here we briefly show that the (generally non-orthogonal) vectors e (1) , . . . , e (D) are indeed linearly independent: . . , D and, correspondingly, the e (α) are linearly independent, thus spanning a D-dimensional plane of instability in component space. This argument is indeed true for all sets of vectors {e (γ) }, for which all corresponding eigenvalues λ (γ) = 0. As we will see later, λ (γ) = 0 only for linearly dependent features, and we conclude that the plane in component space spanned by {e (γ) } γ=1,...,R has dimension equal to the number of linearly independent features.

S4 Systematic expansion of the free energy along a path in density space
The spinodal corresponds to the submanifold in phase space where a homogeneous phase becomes unstable with respect to local fluctuations and the system spontaneously phase separates into two (or more) phases. Critical points on the spinodal occur if these phases become indistinguishable. Here, we define an n-th order critical point as a point where n different phases become indistinguishable and merge into one phase. In a one-dimensional system, an n-th order critical point at densities ρ (cp) manifests as a minimum of order 2n−1 in the tilted Landau free energy There n minima and the n − 1 maxima in between merge into one minimum. We generalize this idea to the multicomponent system by requiring that for an n-th order critical point, the free energy change ∆f N −1 (with respect to ρ (cp) ) along the path of minimal change exhibits a minimum of order ≥ 2n − 1 [5]: for all paths in density space, with equality along an "optimal" path in density space for which ω (m) = Ω (m) ∀m. Here, the density vectors (for which we omit the vector notation) are to be understood in the N − 1 dimensional space. The incompressibility constraint, Similarly, the changes in density need to satisfy if extended to N dimensions. Note that in Eq. S29 we have defined the "derivatives" ω (m) in a way that factors out the densities ρ (cp) i . Furthermore, ω (m) generally depends on the density ρ (cp) . For simplicity, we omit this dependency in the following.
Since the first derivative of the tilted free energy is zero, expanding ∆f N −1 ( ) in terms of the (higher-order) derivatives of where we implicitly sum over repeated indices i 1 , . . . , i l and implicitly assume that A (l) is evaluated at the critical density ρ (cp) here and in the following. Grouping these terms in orders of gives are invariant with respect to permutations of the indices {i j } j=1,...,l , we do not need to explicitly keep track of the indices in ω (mj ) and simplify the previous expression as where we defined the full contraction and furthermore υ (m) := ρ (cp) ω (m) , or in index notation The last sum in Eq. S35 contains products with l factors of υ (mj ) , whose superscripts m j add up to n, i.e. terms of the form υ (1) p1 . . . υ (n) pn with p 1 , . . . , p n ∈ N 0 , subject to the constraints n j=1 p j = l and n j=1 jp j = n. The number of combinatorial possibilities to get a term of the form υ (1) p1 . . . υ (n) pn is l!/(p 1 ! . . . p n !), corresponding to the number of distinct permutations with repeated elements. As a result, Eq. S35 is rewritten as The last sum exactly corresponds to the partial exponential Bell polynomial B n,l (υ (1) , υ (2) , . . .) and we find with the implicit contraction between A (l) and the Bell polynomial B n,l . To find conditions for the occurrence of (higher-order) critical points, we successively determine the optimal vectors Ω (m) by minimizing the coefficients in front of n (only necessary for even n; see below) and setting them to zero one after the other (up to order 2n − 1 for an n-th order critical point as discussed before).

S4.1 Critical point
For the usual critical point n = 2, we have to second order It follows that the Hessian H = A (2) needs to be positive semi-definite and the system is metastable thermodynamically (each critical point lies on the spinodal). Correspondingly, the tangent to the optimal path ρ (cp) Ω (1) coincides with the direction of the instability: Note that, in principle, there could be several directions of instability if the smallest, zero eigenvalue of the Hessian is degenerate. Here, we focus on the non-degenerate case with a unique direction of instability, but we will briefly comment on the degenerate case below in a separate section. Using Eq. S42 with the optimal υ (1) = Υ (1) , we find for the third order where in the second step we used Eq. S42 and the definition of the direction of instability, A ij Υ (1) i = 0, and in the third step Eq. S334. Taken together, apart from critical points lying on the spinodal, the third cumulant (or moment) of the direction of instability needs to be zero at a critical point: Before discussing higher-order critical points, we observe that the third order coefficient, Eq. S43, is independent of any vector ω (i≥2) not yet determined previously, and therefore does not provide a condition on the next higher vector ω (2) .
Instead setting it to zero yields a further condition on ω (1) without the need for minimization. As we will see in the next subsection, this pattern that minimizing the coefficient in front of 2n defines the vector ω (n) and that the coefficient of the next higher order 2n+1 only depends on the previously determined vectors ω (m) , m ≤ n, is true for all values of n. Correspondingly, increasing the order of the critical point by 1, n → n + 1, imposes two more conditions on the critical point and simultaneously defines one additional vector.

S4.2 Higher-order critical points
More precisely, for an n-th order critical point we find the following recursive equation defining vectors Υ (m) , m ≤ n − 1: Here, (A pseudo ) −1 is the pseudo-inverse of the Hessian H = A (2) (which is singular on the spinodal where all critical points lie). Specifically, we define/choose it as where -as noted previously -we assume a non-degenerate maximal eigenvalue of 1 for the covariance matrix, i.e. that λ (γ) < 1 ∀γ ≥ 2. The vectors b (m) , m = 2, . . . , n − 1 satisfy the "orthogonality" relation The conditions of the n-th order critical point, Eqs. S30 and S40, can then be rewritten in terms of the vectors Υ (m) as where we used Eq. S334 to write the full contraction for l ≥ 3 in terms of averages. We will prove the recursion relation, Eq. S45, together with the orthogonality relation, Eq. S47, by induction.
Induction step: n → n + 1 Suppose now that the recursion, Eq. S45, holds for a critical point of order n. We will show that it then also holds for a critical point of order n + 1. The argument will be done in three steps:

1) Orthogonality
We show that if the recursion relation and the orthogonality condition hold for all m = 2, . . . , n − 1, then the orthogonality condition holds for n as well: 2) 2n -term: Independence of higher-order vectors and minimization We demonstrate that the orthogonality condition implies that the coefficients in front of the 2n term, which is the first additional one considered for a critical point of order n + 1, does not depend on any vector υ (i) with i ≥ n + 1. By minimizing the coefficient for 2n with respect to υ (n) , we then recover the next step in the recursion, relating Υ (n) to the previously determined vectors Υ (i) , i ≤ n − 1.
3) 2n+1 -term: Independence of higher-order vectors From the recursion relation for Υ (n) together with the orthogonality condition, it then follows that also the coefficient in front of the 2n+1 term is independent of any vector υ (i) with i ≥ n + 1.
We will go through them one-by-one.
Suppose that Eqs. S45 and S47 hold for a critical point of order n.

1) Orthogonality
To begin with, we use the definition of b (n) to write where we used Eq. S333 and Ω (1) cp = 0 to rewrite the partial contraction of A (l) and B n,l−1 . The partial Bell polynomials satisfy the recursion relation Using this relation, we rewrite Ω (1) b (n) cp as where, for simplicity, we skip the index cp here and in the following. We have n ≥ 2, or n + 1 ≤ 2n − 1, and so n+1 l=3 according to Eq. S48. Thus, where we used that A (2) Υ (1) = 0 (spinodal condition) and the definition of holds by induction and we have Employing Eq. S338 for the product of the Hessian A (2) with its pseudo-inverse and using that b where we used that Ω (1) b (k) = 0 ∀k = 2, . . . , n − 1 by induction. Thus, The recursion relation, Eq. S45 together with the definition of the pseudo-inverse, Eq. S46, implies that ∀i, j = 2, . . . , n − 1 where the last step is due to symmetry i ↔ j. As a result, we have In the second step, we used the variable transformation k → n + 1 − k on the second term. Overall, since 2) 2n -term: Independence of higher-order vectors and minimization: next recursion step Consider now the coefficient in front of 2n , Collecting all terms containing υ (n+1+s) for 0 ≤ s ≤ n − 2 (B 2n,l only includes vectors υ (i) with i ≤ 2n − l + 1 since all its terms contain l factors whose sum of superscripts needs to equal 2n), we find where we used that any vector υ (n+1+s) with s ≥ 0 cannot be in a product together with any other vector υ (n−1−s+t) with t > 0 (otherwise the sum of the superscripts exceeds 2n). In particular υ (n+1+s) with s ≥ 0 cannot occur with a power > 1. Furthermore, since we factored out υ (n+1+s) , in the last sum the number of factors l gets reduced by 1 and the sum of the superscripts by n + 1 + s (effectively setting p n+1+s = 1). We observe that the last sum is exactly B n−1−s,l−1 /(n − 1 − s)!, giving where terms in brackets [] are to be interpreted as vectors with single index. Plugging in the optimal vectors as determined from the conditions for the n-th order critical point, υ (i) = Υ (i) ∀i = 1, . . . , n − 1, we have ∀ 0 ≤ s ≤ n − 2: Furthermore, the recursion relation for the n-th order critical point: can be rewritten as using Eq. S338 and that Ω (1) b (n−1−s) ρ (cp) = 0 for all 0 ≤ s ≤ n − 2 due to the orthogonality relation for the n-th order critical point. Thus, τ n+1+s = 0 (S70) and the coefficient in front of 2n is independent of υ (n+1+s) for s ≥ 0. Correspondingly, the 2n -coefficient is only a function of υ (n) and the previously determined Υ (i) , i = 1, . . . , n − 1, and by minimizing it with respect to υ (n) , the optimal vector Υ (n) can be determined. Specifically, the coefficient in front of 2n , c 2n contains a quadratic term in Υ (n) , a linear term in Υ (n) and a constant term: Using the orthogonality condition for b (n) analogously as in Eq. S56 for b (k) , k = 2, . . . , n−1, where we used the normal matrix and vector [] notation ([] T is the transpose). Since the Hessian A (2) is a positive semi-definite matrix,c 2n (or, equivalently, c 2n ) has a minimum for for any α ∈ R. Here, we choose α = 0 for simplicity. This choice, however, does not influence any of our conditions for critical points as it just corresponds to a rescaling of the path in density space (see below). We thus recover the next step of the recursion relation: The minimized coefficient c min 2n must then equal zero and we find the first additional condition when going from an n-th to and n + 1-th order critical point: As discussed before, this expression is a function of Υ (1) , . . . , Υ (n) only. All terms with higher-order vectors cancel each other. The same is true for the coefficient in front of the 2n+1 term, as we will show next. Therefore, for an (n+1)-th order critical point, only the vectors Υ (1) , . . . , Υ (n) need to be (or can uniquely be) specified.
Taken together, we find the following to additional conditions for an n + 1-th order critical point as compared to an n-th order critical point: which both only depend on the vectors Υ (i) , i = 1, . . . , n given by the recursion relation It follows that and for i = 2, . . . , n where we used Eqs. S333, S335 and abbreviated B i,l (Ω) := B i,l (Ω (1) , Ω (2) , . . .). Thus, we find the following vectors for an n-th order critical point: where we used Eqs. S329, S334 and that B m,2 (Ω) = 1 2 m−1 k=1 m k Ω (k) Ω (m−k) . Note that choosing a different pseudo-inverse by adding a term like βe (1) i e (1) j to our specific choice of the pseudo-inverse, Eq. S46, does not influence our results: Ω (1) is orthogonal to b (n) for any choice of β. Furthermore, note that in order to conclusively show that a certain mixture composition corresponds to an n-th order critical point, it is also necessary to show that the next-higher-order term in the expansion of the free energy is positive for all possible continuations of the optimal path; see Section S9.2 for a discussion in the context of the suggested mean-field q=3-states Potts model.

S4.3 Recursion relation in terms of derivatives of a "Hamiltonian"-like function
Instead of expressing the recursion relation in terms of the partial exponential Bell polynomials, they can also be rewritten in terms of derivatives of a "Hamiltonian"-like function. To see this equivalence, we define the analytic function G( ) with G(0) ≡ 1 (for convenience, see below) via its derivatives at = 0: or alternatively, with G (0) := G and Ω (0) := 1. The function G needs to fulfil or equivalently, N . The recursion relation for an n-th order critical point, Eq. S45, is then rewritten as using Eq. S333. Importantly, the partial Bell polynomials satisfy Faà di Bruno's formula: We use this formula for the function F (x) = 1 − x + log x, whose derivatives fulfil F (k) (G( )) =0 = F (k) (1) = (−1) k+1 (k − 1)! ∀k ≥ 2 and F (k) (G( )) =0 = 0 otherwise. This allows us to express Υ (m) as follows: where in the last step we used Eq. S335. Thus, ∀m = 2, ..., n − 1. (S103) We now express in terms of a Hamiltonian-like function h with h(0) = 0. This specific form satisfies G(0) = 1 and G cp = 1 as desired.
Eq. S103 then requires that the Hamiltonian fulfils Assuming that h is analytic, we thus find (S106) Evaluating this expression and its derivative at = 0 determines the constants of integration: where we used that (S109) This form of the recursion relation has to be combined with the conditions for the critical points, Eqs. S48, which are rewritten as While for general number of features, we did not manage to find an explicit solution to these equations, we did solve them for R = 1. We discuss this solution in the next subsection.

S4.4 Single feature R = 1: Solution of the recursion relation and explicit conditions for higher-order critical points
For a single feature R = 1, the recursion relation for the Hamiltonian h, Eq. S109, reduces to where we abbreviated r (1) = r and have used that

This recursionindeed for all orders n -is solved by
where c is a constant in terms of the vector, i.e. independent of the vector index. Without loss of generality, we choose c = 0 (any constant cancels in G = e h+c / e h+c cp = e h / e h cp ) and we find G = e r e r cp . (S113) In terms of G (or the Hamiltonian h = r), the optimal path is given by (S114) Here, it becomes apparent why we chose G(0) = 1: This choice allows to directly express the optimal path as Taylor series of G in . This explicit representation of the optimal path (vectors) in terms of G also enables us to derive explicit conditions for critical points of arbitrary order n. To that end, we rewrite the conditions for higher-order critical points, Eq. S48, in terms of G: x + x log x, whose derivatives areF (k) (G(0)) =F (k) (1) = (−1) k (k − 2)! ∀k ≥ 2 andF (k) (G(0)) = 0 otherwise, we use Faà di Bruno's formula again and find Using that which can be shown by induction, the conditions for an n-th order critical point, Eq. S118, are given in terms of the cumulants of r: Here, κ i is the i-th cumulant of r. Solving these equations recursively yields: which we briefly demonstrate by induction: Base case n = 2 Eq. S124 gives: Since κ 2 corresponds to the variance of r = √ 2Cs and since κ 2 = 0 can be excluded as this case would correspond to all component types being equal, we conclude that κ 2 = 1 and κ 3 = 0.
Overall, we thus find the following conditions for an n-th order critical point: depending solely on the cumulants κ i of the rescaled spin value r = √ 2Cs. Rewriting them in terms of the cumulants κ = k B T zJ (S137) κ (s) m = 0 ∀m = 3, . . . , 2n − 1. (S138) In the limit n → ∞, the spins s i need to be distributed according to a Gaussian with variance 1/(2C), implicitly requiring an infinite number of component types N → ∞.

S4.5 Choice of α in the recursion relation
In Eq. S75, we have seen that in principle there are several solutions for the recursion relation, namely all recursions of the form where α (n) is a constant and α (1) = 0 (otherwise, all Υ (n) = 0). So far, we have chosen α (n) = δ n1 , and one may wonder whether other choices lead to the same conditions for the spinodal and critical points.
In this section, we show that the choice of α (n) (with α (1) = 0) does not change our results. Instead it just corresponds to an effective rescaling of . To see this, we perform the following steps: The partial exponential Bell polynomials can be defined via the following series expansion: (S147) Thus, B n,k (x (1) , x (2) , . . .) is given by (S148) Using this expression and abbreviating B n,k (x (1) , x (2) , . . .) =: B n,k (x), we rewrite . (S150) As a next step, we observe that the boundaries for j can be extended to 0 and ∞, respectively. For j ≥ k + 1, the term ∞ r=1 α (r) t r r! j is of order O(t k+1 ) since r ≥ 1, and yields a zero contribution when the k-th derivative ∂ k /∂t k is m is of order O(s j+1 ) and is zero when the derivative ∂ j ∂s j is evaluated at s = 0. We thus find Taken together, where we used the definition of Z (i) for i = 1, . . . , n − 1 and that, since m ≥ 2, if s ≥ k in one factor, the terms are of order O(t k+1 ) and yield zero after taking the derivative and setting t = 0, thus concluding part 1. Suppose Y and Z satisfy the recursion relations, Eqs. S143, S145, respectively. We will show by induction that then Base case n = 2 (S158) Induction step: n → n + 1 Suppose Z (i) = i j=1 Y (j) B i,j (α) ∀i = 1, . . . , n − 1, then since Z satisfies recursion S145, or, equivalently, since B n,1 (α) = α (n) , concluding our proof by induction and thereby part 2.
As a last step, we show that Y and Z yield the same conditions for critical points and the spinodal. Suppose that which is proven by a simple induction.
In terms of the optimal path in component space, the choice of α (n) just corresponds to a certain parameterization of the curve: Thus, going from Y to Z leads to a reparameterization of as (S173)

S5 Components of different sizes
In Eq. S7, the densities ρ i can be interpreted as either the volume or the number fraction of component i in the mixture. These are equivalent if all components have the same size and occupy the same number of lattice sites. In realistic mixtures of lipids, proteins and DNA, this assumption might not be satisfied. So how does the size of components modify the phase and critical behavior of the mixture? Addressing this question requires modifying Eq. S7 to include the sizes/lengths l i of components i. Here we do this in a purely mean-field way that neglects correlations between lattice sites arising due to the finite length of the components, as done in Flory's and Huggins' pioneering work for two-component systems [18,19]: where ρ i is now the volume fraction of component i. The number fraction is given by (ρ i /l i )/ j (ρ j /l j ).
Expressing the free energy in terms of the volume fractions ρ i , i = 1, . . . , N − 1 (using that ρ N = 1 − N −1 i=1 ρ i due to the incompressibility) and differentiating the free energy twice yields the Hessiañ Using Woodbury's matrix identity [21] on L, we find for its inverse: Using similar arguments as before, the Hessian is invertible if and only if 1 − U T L −1 U is invertible. This matrix can be rewritten as Due to the factors of l appearing here, the second part is not as straightforwardly written as a covariance matrix, at least not with respect to ρ. However, we can define a new probability measure φ by which satisfies φ m ≥ 0 and N m=1 φ m = ( N m=1 ρ m l m )/l = 1, as desired. In terms of φ, the matrix 1 − U T L −1 U is and so where Cov (φ) is the covariance matrix of the rescaled and shifted features r (γ) with respect to the new probability measure φ: Here, all averages are taken with respect to the weighted volume fractions of the different components: Overall, we find that the Hessian matrix is invertible if and only if 1 −lCov (φ) is invertible. By analogous arguments as before we thus conclude that the condition for the spinodal is whereλ (1) is the variance of the first principal component of Cov (φ) . The direction of instability is determined by looking at the inverse of the Hessian matrixH −1 = L −1 + L −1 U (1 − lCov (φ) ) −1 U T L −1 . We rewrite its components as The covariance matrix Cov (φ) is a real and symmetric matrix and therefore can be decomposed into an orthonormal set of eigenvectorsṼ (γ) corresponding to the descending eigenvaluesλ (γ) : (S186) Using this decomposition the inverse of the Hessian is On the spinodal, whereλ (γ) = 1/l, the inverse of the Hessian is dominated by the term j ), suggesting that the direction of instability is given bỹ This can also be verified by explicitly calculating which is zero on the spinodal, whereλ (1) = 1/l. Looking at the (tilted) free energy change along a path ρ( ) = ρ (cp) + δρ( ) = ρ (cp) + ∞ n=1 n n!Υ (n) centered on a point ρ (cp) on the spinodal, we find up to second order in and all terms are evaluated at ρ (cp) . Analogously to before, if we require that two minima and one maximum merge along the path, we identifyΥ (1) =ẽ (1) .
To determine the condition for the critical point, we then set the third order term in to zero, giving 0 =Ã where we used that N i=1Υ as condition for the critical point.

S6 Negative interaction strengths
So far, we have implicitly assumed that the rescaled features r (γ) are all real and that Cov is a true covariance matrix, which is symmetric and positive semi-definite and has an orthonormal set of eigenvectors. However, these assumptions are only true if the interaction strengths J (γ) between the features are all positive. Our previous analysis is thus restricted to the case where all eigenvalues of the interaction matrix ∼ J (γ) s (γ) 2 are positive and where the pairwise interactions interactions between alike components are always energetically preferred as compared to interactions between dislike components. To resolve this limitation, we now discuss the general case with a combination of positive and negative interaction strengths. Without loss of generality, we assume that the first R + "positive" features are attractive and the remaining R − = R − R + "negative" ones are repulsive: (S196) The rescaled and shifted features can thus be written in terms of real numbersr (γ) ∈ R ∀ γ = 1, . . ., R as r (γ) =:r (γ) ∀γ = 1, . . . , R + r (γ) =: ir (γ) ∀γ = R + + 1, . . . , R, (S197) and the "covariance" matrix Cov (which, in the general case, is in fact a pseudo-covariance matrix) as with the real matrices These matrices quantify the covariances between the subsets of positive (+) and negative (−) features, respectively: Cov In particular, C (++) and C (−−) are true covariance matrices restricted tor (γ) for γ = 1, . . ., R + and γ = R + + 1, . . ., R, respectively, and are positive semi-definite: C (++) , C (−−) 0. The matrix 1 − Cov occurring in the condition for the Hessian to be invertible is thus a 2 × 2 block matrix whose lower diagonal element is 1 + C (−−) 1 and hence invertible. Therefore, 1 − Cov is invertible if its Schur complement with respect to 1 + C (−−) , denoted by (1 − Cov)/(1 + C (−−) ) =: 1 −C, is invertible and has no zero eigenvalue. HereC is given bȳ which is a real and symmetric matrix since C (++) and C (−−) are symmetric and (C (+−) ) T = C (−+) .
Spinodal For small absolute values of the interaction strengths, for which the system should not phase separate, the matrices C (±±/±∓) are small and, consequently, 1 −C is close to 1 and has only positive eigenvalues. The condition for the spinodal is thus that the smallest eigenvalue of 1 −C is 0, or, equivalently, that the largest eigenvalue ofC is 1: where {λ (α) C } α denote the eigenvalues ofC in descending order. Note that the matrix dimension of the Schur complement, 1 −C, Eq. S200, is R + × R + . This illustrates that phase separation relies on features with positive interaction strengths (positive eigenvalues of the interaction matrix). Then again, features with negative interaction strengths modify the propensity to phase separate. More specifically, since (C (+−) ) T = C (−+) and 1+C (−−) (or, equivalently, (1+C (−−) ) −1 ) are positive semi-definite, ≥ 0 for any real vector v. For two real, symmetric, positive semi-definite matrices A and B, the largest eigenvalue of A − B is smaller than the largest eigenvalue of A: Let v be the normalized eigenvector corresponding to the largest eigenvalue λ max . Thus, the negative features always tend to lower the eigenvalue ofC and thereby prevent phase separation. The extent to which they influence the phase behavior depends on the relative correlations between all features: Suppose v (++) is the normalized eigenvector of C (++) corresponding to the largest eigenvalue. Or, in other words, assume that the dominant combination of positive features that drives phase separation is R (++) := As a result, if the dominant combinations of positive features driving phase separation correlate only weakly with the negative features, r (β+R + ) R (++) ρ ≈ 0 ∀β, the thermodynamic stability is barely modified by the presence of the latter: . Conversely, if for each dominant combination of positive features, there is a highly correlated negative feature, the two effects counteract each other and, depending on the relative strengths and correlations, the effective interactions between the components may be small -the mixture will not phase separate. In a similar spirit, depending on the correlation structure, the addition of repulsive features may shift the relevance of the different attractive features with respect to their role for phase separation. Taken together,C can be interpreted as representing a multicomponent system with R + positive features with effective, reduced interactions integrating the interplay of the positive and negative features in the original system.

Direction of instability
If the Schur complement 1 −C of 1 − Cov with respect to 1 + C (−−) is invertible, the inverse of 1 − Cov exists and is given by For our purposes, this expression is conveniently rewritten as The inverse of the Hessian matrix then exists and is given by where for all i, j = 1, . . . , N − 1 we have defined the deviations of the rescaled features from the mean and their weighted sum Since C (++) and C (−−) are symmetric, C (−+) = (C (+−) ) T and all matrices are real,C is real and symmetric and thus has an orthonomal set of R + eigenvectors φ (γ) corresponding to the descending eigenvaluesλ (γ) : β and we find Close to the spinodal, the inverse of the Hessian is dominated by the term 1 j . Similarly to before, we thus identifyē (1) with the eigenvector of the Hessian to eigenvalue 0 at the spinodal: Critical point Expanding the free energy along a path ρ( ) up to third order in , we find an analogous condition for the critical point as for the original model: whereē i .

S7 Linearly dependent features
In this section, we show that adding linearly dependent features does not change the results of our analysis. We restrict our discussion to the case of only positive interaction strengths and find that there is an additional zero eigenvalue λ (Γ) = 0 (all others remain the same) but the corresponding additional eigenvector V (Γ) leads to a zero vector E (Γ) = 0, thus not changing our predictions.
To simplify notation, we assume without loss of generality that 2C (γ) = 1 (absorbing the interaction strengths into the spin values) and that s (γ) N = 0 (shifting the spin values; see also Section S8). Thus, r (γ) = s (γ) . Suppose now that the R-th feature can be written as a linear combination of the other features: with ν ∈ R R−1 and s (γ) = (s where we defined the matricesS αi = s (α) i for i = 1, . . . , N , α = 1, . . . , R − 1 and M αβ = δ αβ + ν α ν β for α, β = 1, . . . , R − 1. The latter matrix is real and symmetric and can thus be written in terms of a diagonal matrix D and an orthogonal matrix Q: Further, since M is a rank-1 correction of the identity matrix, the eigenvalues of M are 1 + ||ν|| 2 and 1 (R−2 times degenerate): The corresponding normalized eigenvectors are thus ν/||ν|| =: w (1) and a set of R−2 normal vectors w (γ) , γ = 2, ..., R − 1 with w (γ) ⊥ w (δ) for all γ = δ. In terms of these, the matrices are Moreover, we can express the interaction matrix as where and the square of the diagonal matrix is We thus have two representations of χ, one in terms of R linearly dependent features and the other one where we have integrated out the last feature: where S γi = s (γ) i has dimensions R × N and P δi =: p (δ) i has (R−1) × N . Do these two representations yield the same conditions for the spinodal and the critical points? To answer these questions, we consider the respective covariance matrices Cov (s) and Cov (p) and their eigenvalues and eigenvectors next.
To conclude, we will check that these different sets of eigenvectors give the same vectors E (γ) in the space of components.
To this end, we observe that for γ = 1, . . . , R − 1 Using the definitions ofQ and P and that Q is an orthogonal matrix, we find For β = 1, . . . , R − 1, this expression simplifies to and, similarly, for β = R, we have i . Overall, we conclude that Combining this with Eq. S239, we findẼ showing that indeed both representations lead to the same vectors in the space of components. Taken together, we find that the addition of linearly dependent features does not change the eigenvalues of the covariance matrix, nor the resulting vectors in component space. Our results for the spinodal and the critical points thus do not depend on whether or not the features are chosen as linearly independent.
Finally, we observe that if the covariance matrix has a zero eigenvalue λ (R) = 0, then the corresponding vector E (R) = 0. Furthermore, either one features is deterministic, r 1 ∀i, and therefore does not lead to any effective interaction between the components, or the features are linearly dependent. To see this, we note that since the covariance matrix is real and symmetric, we can write it as where T is orthogonal and ∆ is a diagonal matrix with ∆ RR = λ (R) = 0. As a result, where Since the variance of z can only be zero if z is deterministic/constant, we have For the vector E (R) we then find Since V (R) is an eigenvector of Cov, it is not equal to the zero vector and, unless r (δ) − r (δ) ρ = 0 for some δ, we conclude that the (shifted) features {r (γ) − r (γ) ρ } γ=1,...,R are linearly dependent. Note that, in general, the space spanned by the feature vectors of a mixture with N components is at most N − 1dimensional: N points always lie in an N − 1-dimensional Euclidean subspace. Thus, the covariance matrix has at most N − 1 non-zero eigenvalues, in accordance with the Hessian matrix having dimensions (N − 1) × (N − 1).

S8 Invariance under global rotations and translations of the rescaled feature vectors
In this section, we show that the system is invariant under global rotations or translations of the rescaled feature vectors. We proceed in two steps: First, we demonstrate that the tilted free energy only changes by a constant and thus describes the same physical system. We then show that, accordingly, the conditions for the spinodal and (higher-order) critical points do not change under global rotations or translations of the rescaled feature vectors.
To simplify notation, without loss of generality we consider the case when the interaction strength is absorbed into the spin values and we have 2C (γ) = 1∀γ. In this case, the rescaled and shifted feature vector just corresponds to and the system is invariant to global rotations and translations of the rescaled feature vector: where O is an orthogonal rotation matrix, OO T = O T O = 1 and det(O) = 1, and τ is a constant shift 4 . To see this, note that the interaction matrix is where we used that , the second term is linear in ρ and the third term corresponds to a constant. Neither term thus contributes to the tilted free energy and the thermodynamic behavior of a system with feature vectors s i is the same as for a system with feature vectors s i . Correspondingly, the conditions we derive are the same in both cases: The covariance matrices are related via the rotation matrix (the translation cancels since a constant shift in the mean does not change the covariance): Since the rotation matrix is orthogonal, this implies that the covariance matrices have the same eigenvalues and their eigenvectors are related via As a result, the relative enrichment is the same: where we used that r (α) = O αβr (β) and so O αγ r (α) = O αγ O αβr (β) = δ γβr (β) =r (γ) (implicit sum convention). Since, furthermore, Eq. S91 implies by induction that Ω (k) =Ω (k) ∀k and, thus, according to Eq. S93 all conditions are the same. Note that the same is true for reflections s

S9 Several directions of instability
In this section, we briefly comment on a second series of higher-order critical points, arising from several directions of instability, i.e., in the case when on the spinodal the largest eigenvalue of the covariance matrix is degenerate. 4 If C (α) = C (β) for some α, β, then the system is invariant to

S9.1 Codimension of the critical submanifold
Suppose that the covariance matrix has a degenerate maximal eigenvalue of 1: Then, according to Eq. S26, we have with the directions of instability ρ as defined previously in terms of the orthonormal set of eigenvectors V (γ) of the covariance matrix. Note that since the covariance matrix is real and symmetric, eigenvectors corresponding to degenerate eigenvalue can still be chosen orthogonal, which we assume in the following. The Hessian is a linear operator and so we have for any linear combination demonstrating that the directions of instability are linearly independent. They thus span a D-dimensional subspace of the component density space along which the second derivative of the free energy is zero.
What are the conditions that m ≤ D + 1 local minima of the free energy merge along this D-dimensional subspace? One possibility is that simultaneously along each of m − 1 linearly independent directions of instability, two minima merge, i.e. that a usual critical point arises simultaneously along m − 1 linearly independent linear combinations D γ=1 µ (i) γ e (γ) , i = 1, . . . , m − 1 of the directions of instability. In analogy to our previous analysis, where we considered the optimal path in the density space, this requires that the third cumulant of these linear combinations has to be zero: However, since it should not be possible to minimize the free energy along any other direction (in the D-dimensional subspace where the Hessian is zero), the third cumulant needs to be zero for all linear combinations of the D directions of instability: and therefore there will be simultaneous ordinary (or higher-order) critical points along all D directions of instability, thus making the critical point order D + 1 (or higher). The above conditions are only fulfilled if for all possible combinations of γ 1 , γ 2 , γ 3 = 1, . . . , D. Since the distinct combinations correspond to drawing a number in 1, . . . , D three times without considering the ordering, there are 3+D−1 3 conditions that need to be satisfied. What is the codimension resulting from D eigenvalues being 1? As discussed in [26], for real symmetric matrices (such as the covariance matrix or the Hessian), the set of matrices with a D-fold degenerate eigenvalue has codimension (D+2)(D−1)/2. Since in our case, these eigenvalues have to be equal to 1, this gives a codimension of (D + 2)(D − 1)/2 + 1 = D(D + 1)/2 = from the third cumulants being 0 5 . Note that for D directions of instability (at least) D + 1 phases will become indistinguishable at a critical point: either all third cumulants in the D-dimensional subspace of directions of instability are zero (and at least D + 1 phases merge) or there is a direction along which the free energy can be lowered and the point is not thermodynamically stable. Of course, one could also imagine combining this series of higher-order critical points arising from merging of minima along several directions of instability with the series as discussed before, where several minima merge along a single direction of instability. This combination will, however, not be discussed here. To conclude this section, we will briefly touch upon one idea for an effective mean-field theory for the two-dimensional q=3-states Potts model with two direction of instability.

S9.2 A mean-field version of the q=3-states Potts model
The regular mean-field version of the Potts model can be expressed in terms of our framework by a N = 3-component system with the following R = 3 features: where i, γ = 1, 2, 3. Here and in the following, we directly incorporate the interaction strengths J (γ) (and other constant factors) into the spin values and only keep the dependency on the temperature T , effectively setting zJ (γ) /k B = 1. Then, the covariance matrix is given by At the symmetry point ρ 1 = ρ 2 = ρ 3 = 1/3, we have with eigenvalues λ (1) = λ (2) = a 2 /3/T and λ (3) = 0. The spinodal is thus located at a 2 = 3T and the system exhibits two directions of instability with corresponding (orthonormal) eigenvectors of Cov Importantly, not all combinations of third cumulants of these directions of instability are zero. While two are indeed zero, the other two are equal to ±a 3 /3/ √ 6, thus leading to thermodynamic instability along certain directions in the 2-dimensional subspace of directions of instability. Hence, as noted previously [28,29], the symmetric state does not exhibit critical behavior, in contrast to the two-dimensional q=3-states Potts model [27]. Ultimately, this lack of a continuous transition is due to the fact that for D = 2, seven conditions need tuning while any N = 3, R = 3 model only exhibits five degrees of freedom, two ratios of densities and three interaction parameters: The three feature vectors necessarily lie in a plane, thus making the effective number of features =2 (a rotation of the feature vectors onto e.g. the x-y-plane does not While this regular mean-field version of the q=3-states Potts model does not exhibit a critical point, here, we propose a slightly more complex version with N = 6 and R = 2, which preserves the S 3 symmetry group of the q=3-states Potts model but represents each "lattice component type" by two "effective component types".
More specifically, we consider a mixture with the following feature vectors: corresponding to two equilateral triangles centered around the origin (see Fig. S1). For simplicity, we absorb all the interaction strengths and constants into the spin values, and assume that 2C (1) = 2C (1) = 1/T , retaining only the dependency on temperature. Furthermore, we focus on the symmetric state where all component types belonging to the same equilateral triangle have the same densities: For this symmetric case, the average of all spins vanishes, s (1) = s (2) = 0, and the covariance matrix is given by The eigenvalues can be read off as with corresponding eigenvectors The largest eigenvalue is thus two-fold degenerate and on the spinodal, there are two directions of instability. According to our previous analysis, we thus need the third cumulant of any linear combination of to be zero in order for the system to be at its critical point (where then the 3 symmetric phases merge into a single phase): In addition to the condition for the spinodal, we thus find that at the critical point. Since x and 1 − x denote the combined densities of component types 1, 2, 3 and 4, 5, 6, respectively, (and since we have seen previously that a three component system is not enough,) we need x ∈ (0, 1), yielding the additional constraint Choosing the density as in Eq. S281, the spinodal condition is where r = a/b is the ratio between the linear dimensions of the two "feature triangles". In principle, equality S283 is satisfied by a one-dimensional submanifold of (b, r)-space. We do, however, need to ensure that the symmetric state is thermodynamically stable (at least to local fluctuations) and that up to the next order in the free energy expansion (the fourth order), the free energy change is positive for all possible paths in density space. The fourth order term in the free energy expansion is (see Eq. S40) Here, we have used that the paths of minimal change in free energy have a tangent within the plane of instabilities: Thus A (2) (ρE) = 0 and the fourth order term is rewritten as This quadratic form exhibits a minimal value although A (2) is only positive semi-definite: Analogously to the case with a single direction of instability we can define the pseudoinverse of A (2) by It satisfies In particular, where we used Eq. S334 and that the third cumulants (E) 2 E (δ) = 0 ∀δ. As a result, F has a minimum, which is given by Here κ is the fourth cumulant of E and we used Eq. S334 together with S333. Overall, we thus find that the fourth order term is positive for all optimal paths only if κ (E) 4 < 0. Rewriting the fourth cumulant gives Since ab < 0, this expression is negative only if We thus find the following condition for the ratio r between a and b: We note that this range is invariant under r→1/r, as expected from relabelling a↔b. Taken together, this N =6/R=2 model with threefold symmetry exhibits a critical point where three phases merge if b and r = a/b are chosen as We suggest that any such model might therefore be a candidate for a mean-field version of the two-dimensional q=3-states Potts model which retains a continuous phase transition. Interestingly, we first tried a version with R = 3 features and N = 6 component types whose spin vectors don't lie in one plane but instead in two parallel planes: In this case, it is not possible to simultaneously satisfy the spinodal and third cumulants condition and to have a positive fourth order term for all possible paths. Ultimately, this is due to the fact that in this case where the features space is three-instead of two-dimensional, the free energy can be lowered by choosing a vector of second derivatives ω (2) that leaves the plane of instability. Looking ahead, it will be interesting to see whether a related model with R = 3 and three directions of instability could serve as a mean-field version of the two-dimensional q=4-Potts model, which also exhibits a continuous transition in two dimensions but not in regular mean-field theory. One possibility might be to consider feature vectors spanning a regular tetrahedron and rotations thereof. Furthermore, it would be enlightening to understand whether such types of models are not possible for D > 4 (for the two-dimensional q>4-states Potts models, the transitions are indeed discontinuous).

S10 Coarse-graining a multicomponent mixtures as a binary mixture
In this section, we discuss one idea how to coarse-grain a multicomponent mixture as a binary mixture of two components with effective features, while conserving two local thermodynamic properties, the spinodal line and the line of (ordinary) critical points. As we have seen before, the conditions for the spinodal and ordinary critical points only depend on the first principal component of the covariance matrix (note that we restrict our discussion to the case of a non-degenerate maximal eigenvalue of the covariance matrix). Hence, if only those two characteristics are to be preserved, only the distribution of features along the first principal component is of relevance and it might be enough to consider a binary system with the same first principal component. We therefore suggest to coarse-grain the multicomponent mixture (denoted by regular symbols) as a binary mixture (denoted by tildes) whose components i = 1, 2 exhibit feature vectors along the 1st principal component of the original system: where x i is the distance of component i along the 1st principal component and γ = 1, . . . , R denotes the different features (note that it would equally be possible to do this coarse-graining procedure with a single feature, corresponding to a rotation of this feature vector onto the first axis). This binary mixture has a covariance matrix where ϕ i is the density of component type i = 1, 2 of the binary mixture: ϕ 2 = 1 − ϕ 1 . We observe that Cov is a rank-1 matrix and, as desired, the first principal component of the original multicomponent mixture, V (1) , is also an eigenvector of Cov (corresponding to the only non-zero eigenvalue, . Similarly, for the third cumulant of the direction in component space corresponding to this eigenvector we find where we used thatr Using that x ϕ = ϕ 1 x 1 + ϕ 2 x 2 = x 2 + ϕ 1 (x 1 − x 2 ), the third cumulant of x (and thus of the potential direction of instability) is rewritten as One possibility to maintain the location of the spinodal and critical line while coarse-graining is to preserve the second and third moment along the first principal component, i.e. to require that where . . . ρ denotes the average in the original space of densities of the multicomponent mixture. As a result, which is solved by For x 1 − x 2 we then find Taken together, we find the following mapping For the density we chose the positive sign without loss of generality.
To proceed we note that x 1 − x 2 is related to the Flory-Huggins interaction strengthχ: which leads to the following conditions for the composition ϕ := ϕ 1 and the interaction strength: At the critical point, where µ 3 = 0 and so τ → ∞, we have We hence find for the distance of the interaction strength from its critical value and for density difference ∆ϕ := ϕ 1 − ϕ 2 = 2ϕ − 1 where the approximation is valid close to the critical point. Taken together, coarse-graining procedures based on such (or similar) arguments can conserve local properties of the free energy like the spinodal and critical points. Whether (and in what cases) the same can be said for global characteristics like the binodal is, however, not clear a priori and is an interesting question for future research. S11 Generation of the multicomponent mixtures illustrated in Fig. 2

of the main text
To illustrate the role of the skewness of the distribution of the relative enrichments δρ i /ρ i for the occurrence of critical behavior, we generated the multicomponent mixtures as displayed in Fig. 2 of the main text from random distributions with zero and non-zero skewness, respectively. More specifically, we followed the procedure as outline next (for simplicity, we assume 2C (γ) = 1 ∀γ).

S11.1 Mixture 1: zero skewness
We first drew N =1000 two-dimensional feature vectors (component types) from a multivariate Gaussian distribution with zero mean and covariance matrix σ 2 1 0 0 σ 2 2 with σ 1 = 1, σ 2 = 0.3. We then randomly generated the corresponding densities of the component types by drawing N =1000 exponentially distributed random variables (with rate λ = 100) and by normalizing the vector so that its sum (1-norm) equals 1 6 . In order to enforce a variance of 1 along the first principal component and to illustrate the case when the first principal component is not aligned with a feature axis, we proceeded as follows: 1. We rotated all feature vectors in a way that the first principal component aligns with the x-axis (by multiplying them by the matrix of eigenvectors of the covariance matrix).
3. We rotated all feature vectors (and thereby the first and second principal component) counterclockwise by an angle ϑ = π/6, thus aligning the first principal component V (1) with cos(ϑ) sin(ϑ) .
The relative enrichment δρ i /ρ i of component type i was then determined by taking the scalar product of the difference of the respective feature vector from the mean of the feature vectors (weighted by the densities) with cos(ϑ) sin(ϑ) . In mixture 1 displayed in Fig. 2 of the main text, the skewness is µ 3 ≈ 0.022, thus very close to zero (it is non-zero due to the finite size of the sample). Note that for the discussion there, we treat it as being equal to zero.

S11.2 Mixture 2: non-zero skewness
For mixture 2, the densities of the component types were generated as for mixture 1 from an exponential distribution with rate λ = 100 and proper normalization, corresponding to a uniform distribution of the mixture composition on the N −1-simplex. In order to achieve a non-zero skewness for mixture 2, we randomly generated the feature vectors as follows: 1. Feature 1 (the first component of the feature vectors) was generated as a sum of two independent random variables, a Poisson distribution with variance σ Poisson = 1.5 and a Gaussian distribution with zero mean and variance σ Gaussian1 = 0.1.
2. Feature 2 was generated independently of feature 1 from a Gaussian distribution with zero mean and variance σ Gaussian2 = 0.03.
3. To ensure a variance of 1 along the first principal component V (1) and to align it with the x-axis, we then rotated and rescaled them as described for mixture 1: We multiplied them by the matrix of eigenvectors of the covariance matrix (obtained after steps 1 and 2) and then rescaled the first and second components of the feature vectors by σ 1 /σ PC1 and σ 2 /σ PC2 , respectively (again σ 1 = 1 and σ 2 = 0.3).
The relative enrichment δρ i /ρ i of component type i was then determined analogously to mixture 1 by taking the scalar product of the difference of the respective feature vector from the mean of the feature vectors (weighted by the densities) with cos(ϑ) sin(ϑ) , where ϑ = 0 for mixture 2. Mixture 2 displayed in Fig. 2 of the main text exhibits a skewness of µ 3 ≈ 0.99.

S11.3 Invariance under constant shifts and rotations
Note that all the physical results are invariant under constant shifts or rotations of the set of rescaled feature vectors, see Section S8. We tried to illustrate this invariance by choosing a principal component not aligned with any feature axis for mixture 1 and by generating a mixture with non-zero (weighted) mean for the feature vectors for mixture 2.

S11.4 Binary system
The binary systems for both mixtures were determined as described in the previous section with µ 2 = 1 and µ 3 as determined from the empirical distributions of the relative enrichments (weighted by the densities of the component types). In particular, we chose the feature vectors as where τ = 4(µ 2 ) 3 /(µ 3 ) 2 . Furthermore, µ 1 is the (weighted) mean of the projections of the feature vectors of the original N =1000 mixture along the first principal component; µ 1 is non-zero for mixture 2. Note that, as noted previously, the constant shift by µ 1 V (1) leaves the physical characteristics of the binary system invariant. Finally, the densities were chosen as

S11.5 Histograms
For both histograms, we chose a bin width of 0.1. S12 Useful formula S12.1 Properties of the Hessian Hessian acting on a vector For the Hessian A (2) , Eq. S9, and a vector Z with Z (N ) ρ = 0 we find Proof: since r implicitly summing over all repeated indices on the left-hand side.
Proof: Using Eq. S326 we have (S332) S12.2 (Partial) Contraction of higher-order derivatives of the free energy where we used that ρ N = 1 − P and Z (k) (N ) ρ = 0. where we used the partial contraction, Eq. S333, in the first equality, and that Z (n) (N ) ρ = 0.