Derivation and analysis of a phase field crystal model for a mixture of active and passive particles

We discuss an active phase field crystal (PFC) model that describes a mixture of active and passive particles. First, a microscopic derivation from dynamical density functional theory is presented that includes a systematic treatment of the relevant orientational degrees of freedom. Of particular interest is the construction of the nonlinear and coupling terms. This allows for interesting insights into the microscopic justification of phenomenological constructions used in PFC models for active particles and mixtures, the approximations required for obtaining them, and possible generalizations. Second, the derived model is investigated using linear stability analysis and nonlinear methods. It is found that the model allows for a rich nonlinear behavior with states ranging from steady periodic and localized states to various time-periodic states. The latter include standing, traveling, and modulated waves corresponding to spatially periodic and localized traveling, wiggling, and alternating peak patterns and their combinations.


I. INTRODUCTION
Active particles such as animals, bacteria, and artificial microswimmers, which transform energy into directed motion, are of significant importance to modern physics both due to their remarkable collective nonequilibrium behavior and their significant technological and biological relevance [1][2][3][4].Of particular interest in current research are mixtures of active and passive particles, which can exhibit remarkable dynamics.An example is motility-induced phase separation (MIPS) [5], i.e., liquidgas phase separation in a system of active particles that can occur even for purely repulsive interactions.The state diagram for MIPS has been found to be affected by the presence of passive particles [6].Further phenomena exhibited by active-passive mixtures include separation of active and passive particles [7], Hopf bifurcations leading to moving clusters [8], collapse in a circular shear cell [9], transitions from gas-like behavior to oscillations and collapses [10], and the existence of different demixing types [11].This has also led to an increased interest in theoretical models describing such mixtures [8,9,[12][13][14].Moreover, topological defects in active nematics can be modeled as mixtures of active and passive particles [15,16].The study of active-passive mixtures in polymer systems is of biological relevance as it allows one to describe phase separation in DNA strands [17,18].Research on such mixtures also has a variety of possible applications, such as controlling the properties of passive materials using active dopants [19][20][21][22] and effects of swimming microorganisms such as algae on their environment [23][24][25].Finally, whenever the active particles are immersed in a passive solvent -as is the case in almost all experimental or theoretical scenarios where active matter is considered -we are, strictly speaking, already dealing with a mixture of active and passive particles.Note that some of the mentioned properties of active-passive mixtures do also occur in "standard" active matter models without the passive component.For instance, although in active phase field crystal models [26] (see below) the onset of motion of clusters normally occurs via drift pitchfork bifurcations [27,28], Hopf bifurcations can also be responsible [29].
When considering PFC models for mixtures, it is interesting that microscopic derivations [75,76] give relatively complicated coupling terms, whereas models proposed on a phenomenological basis typically employ a very simple coupling [70,77,78,89].A similar observation can be made for active PFC models: The derivation of the active PFC model from DDFT in Ref. [84] is focused on the dynamical part, whereas the free energy was simply imported from phenomenological models [26].Microscopic treatments, on the other hand, show that the free energy in a PFC model with orientational dynamics can be very complex [80][81][82][83].Consequently, a systematic microscopic derivation of an active PFC model for a binary mixture is required to investigate both what the general form looks like and which assumptions are required to recover simple phenomenological models, thereby providing insights into their range of applicability and into ways to improve them.
The aim of the present work is twofold: First, we derive a general PFC model for a mixture of active and passive particles from a microscopic DDFT and obtain from it a simpler minimal model (which we refer to as model 1 ), thereby revealing both the necessary approximations and possible generalizations.Second, we investigate the active binary PFC model using linear and nonlinear methods.It is found to have interesting properties that differ from those of simple active PFC models and passive models for binary mixtures.
This article is structured as follows: The minimal binary active PFC model is introduced in Section II.A general PFC model is derived from DDFT in Section III.Approximations leading to the minimal model are discussed in Section IV.In Section V, we perform a linear stability analysis of the PFC model.Nonlinear results are shown in Section VI.We conclude in Section VII.

II. GOVERNING EQUATIONS
We start by introducing the central model of this work based on phenomenological arguments.In the active PFC model [26,27,29, 84], a system is described using a conserved scalar field ψ( r, t) (depending on position r and time t), which describes the dimensionless deviation from a mean particle number density, and a nonconserved vector field P ( r, t), which measures the polarization.They follow the dynamic equations Here, v 0 is the activity of the particles and F is a free energy functional.This functional can be written as where is a Swift-Hohenberg-type free energy [93,96] in d spatial dimensions and is the orientational contribution to F .For C 1 < 0 and C 2 > 0, the system exhibits spontaneous polarization.Most treatments only consider the case C 1 > 0 and C 2 = 0 where no spontaneous polarization occurs [27,84].In addition to the parameters C 1 and C 2 , the free energy depends on the parameters ψ (scaled shifted temperature), q 2 ψ (critical wavenumber), and ψ (dimensionless mean density).We mark these parameters with a subscript ψ since they can differ for the components of a mixture, for example if its components have different freezing temperatures [77].
For the case v 0 = 0, Eq. ( 1) reduces to the standard PFC model for passive systems [93,97].It has the form of a gradient dynamics, where the system evolves toward the minimum of a free energy.The passive PFC model shares this property with DDFT [50] and many other soft matter models [98][99][100][101][102].In the active case, however, the coupling terms −v 0 ∇ • P and −v 0 ∇ψ in Eqs. ( 1) and (2) break this structure since they are nonvariational, i.e., they cannot be obtained from a free energy.This is a typical feature of active matter models [32].
The active PFC model given by Eqs. ( 1) and ( 2) is a two-field model for a single species of active particles.The scalar field ψ and vector field P describe the density and orientational order (polarization) of the single particle species of a one-component system.A PFC model can, however, also describe mixtures of different particle species employing a larger number of fields.For a simple case of a binary mixture of passive particles, see Refs.[76,77].In contrast, here, we consider a binary mixture of passive particles with density φ and active particles with density ψ and polarization P .In this case, we have to add to Eqs. ( 1) and (2) the equation Note that no modification of Eq. ( 1) is required, since the coupling between the different particle species is purely variational (for a case of nonvariational coupling, see Section 4.2 of Ref. [89]).The variational coupling is introduced by adding a mixed term F coup [ψ, φ] to the free energy (3), which now takes the form It is common [70,77,78] to use the simplest nontrivial form where a is a coupling constant.Inserting Eqs. ( 4), ( 5), (7), and (8) into Eqs.( 1), ( 2), and ( 6), setting C 2 = 0, and carrying out the functional derivatives gives the equations of motion where by convention we have explicitly introduced the mean densities ψ and φ as parameters, i.e., d d r ψ = 0 and d d r φ = 0. From here on, we refer to the model given by Eqs. ( 9)- (11) as "model 1".It is closely related to active binary PFC models used in previous works [12,89], although there are small differences (Ref.[89] uses ψ = φ , C 2 = 0, and q ψ = q φ = 1, Ref. [12] uses a nonlinear instead of a linear coupling term).

III. DERIVATION OF GENERAL MODEL
After having introduced model 1 on a phenomenological basis, we now aim for a systematic derivation from the microscopic dynamics of the particles.In doing so, we will make use of a variety of previous results on the derivation of PFC models from DDFT.The "standard" case of a one-component system with no orientational degrees of freedom is discussed in Refs.[74,93,95].Binary systems were considered in Refs.[75,76], orientational degrees of freedom in Refs.[80][81][82][83]103], and active particles in Refs.[26,84].
The starting point is the DDFT where ρ ψ and ρ φ are ensemble-averaged one-body densities (this is why no noise terms are required [50,104,105]), β is the thermodynamic beta, with constants D and D ⊥ , dyadic product ⊗, and unit matrix 1 is the translational diffusion tensor of the (active) field ρ ψ , D R its rotational diffusion constant, D the translational diffusion constant of the (passive) field ρ φ , u(ϕ) = (cos(ϕ), sin(ϕ)) T with polar angle ϕ the orientation, F an equilibrium free energy functional, and v the active propulsion speed.The density ρ ψ gives the probability of finding a particle of species ψ with orientation u at position r at time t, whereas ρ φ gives the probability of finding a particle of species φ (assumed to have no orientational degrees of freedom) at position r at time t.Equation ( 12) is a standard active DDFT for uniaxial particles [52] in two spatial dimensions and identical to the one employed in Ref. [84].The more general case of active particles with arbitrary shape in three spatial dimensions was considered in Ref. [53].Equation ( 13) is a standard DDFT for passive particles.
As discussed in Refs.[46,48,53], DDFT can be derived systematically from the microscopic Langevin equations for active and passive particles using the "adiabatic approximation", where it is assumed that the pair correlations in the nonequilibrium system are equal to those of an equilibrium system with the same one-body density.
This assumption, which is required to close the equations of motion for the one-body density, is not exactly true, but it is a good approximation for many systems of interest [50].
To derive the corresponding PFC model, we need to approximate both the dynamical equations and the free energy.Although these approximations are, as shown by Archer et al. [95], intimately connected, it is helpful to consider them separately.We start with the free energy.In (D)DFT, it is (ignoring external potentials) given by Here, F id is the ideal gas free energy.It is known exactly and given by1 where Λ ψ and Λ φ are the (irrelevant) thermal de Broglie wavelengths, which may differ for the two species.
The second term in Eq. ( 14) is the excess free energy F exc .It is not known exactly and has to be approximated.
A systematic way of doing this is a functional Taylor expansion [43,50,82,93] around a homogeneous state up to fourth order, which gives We have used the direct correlation functions [93] with i ∈ {ψ, φ}, as well as the deviations from the homogeneous reference densities ρ i .Irrelevant zeroth-order and first-order contributions have been ignored.(Zeroth-order contributions have no effect at all, first-order contributions shift the chemical potential by a constant [106], but vanish after applying the gradient operator.)Next, we have to define the fields ψ, P , and φ in terms of the microscopic densities ρ ψ and ρ φ .In principle, there are two possibilities for a physical interpretation of the order parameters in a binary PFC model: The first option is to introduce a shifted rescaled total density field n 1 and a rescaled density difference field n 2 , defined as n 1 = (ρ φ + ρ ψ − ρ)/ρ and n 2 = (ρ ψ − ρ φ )/ρ, where ρ = ρψ + ρφ is the total homogeneous reference density [73].In this case, n 1 describes the structure and n 2 the composition of the particle distribution [93].Here, however, we choose the second option [76] of letting ψ and φ describe the shifted rescaled densities of the two different particle species, i.e., the different fields in the theory directly correspond to the two different particle types.This has the advantage that we can clearly distinguish between an active field ψ and a passive field φ.The existence of these two options shows the importance of providing a microscopic derivation for the binary PFC model, since otherwise the physical interpretation of the fields and therefore the implications of theoretical results for experiments are not clear.
For the field ρ ψ , we use the Cartesian orientational expansion [107] where is the rescaled shifted density of the active particles and is their local polarization.The expansion (19) is equivalent to a Fourier series [107,108], which is widely used in derivations of field theories for active particles [8,37,38].By truncating the expansion (19) after the polarization term P , we are neglecting the contribution of the nematic order to the free energy (as is common in active PFC models).More general expressions including the nematic order can be found in Refs.[81][82][83].For ρ φ , where an orientational expansion is not required since this field has no orientational dependence, we simply define The next step is the expansion of the ideal gas free energy.If we insert Eqs.(19) and (22) into Eq.( 15) and expand the result up to fourth order (to allow for stable crystals) in both fields, we find where we have dropped irrelevant terms of zeroth and first order.The excess free energy is more difficult to treat.We follow Ref. [83] and start with a Fourier expansion of the direct correlation functions.Our treatment of symmetries parallels that of Refs.[8,[36][37][38][39][81][82][83], with differences arising from the fact that not all particles have an intrinsic orientation.Let us consider the direct correlation functions c 1,..., n ({ r i }, {ϕ i }), where 1 , . . ., j = ψ and j+1 , . . ., n = φ.Translational and rotational invariance implies that they can be parametrized as Note that there is no dependence on φi for j = 0 (if i = φ for all i in Eq. ( 17), then no functional derivatives with respect to functions depending on ϕ are taken on the right-hand side, such that the direct correlation function cannot depend on ϕ) or j = 1 (since then there is no angle ϕ 1 and thus Eq. ( 26) is undefined).Moreover, translational and rotational invariance implies that c φφ also does not depend on φRi (there is no angle ϕ R2 in this case, such that, since j = 0, Eq. ( 25) is undefined).
For our purposes, it is sufficient to consider the terms with m i = 0 and m i = 0, ±1 in Eq. (27).The reason for why we only need m i = 0, ±1 is straightforward: We only consider orientational order parameters of zeroth and first order, which is equivalent to considering only the zeroth-and first-order contributions in a Fourier expansion of the one-body density [107].Due to the orthogonality of the basis functions in a Fourier expansion, higher-order contributions of the direct correlation functions will then vanish after performing the angular integrals in Eq. (16).
We now replace ∆ρ ψ → ρψ (ψ + P • u) and ∆ρ φ → ρφ φ in Eq. ( 16).If we also insert Eq. ( 27) into Eq.( 16) and perform the integrals over ϕ, ϕ 1 , . . ., every term will have a prefactor c m,m 1 ,..., n ({R i })e −im•ϕ R .The next step (which is helpful to illustrate why we can restrict ourselves to m i = 0) is a gradient expansion, which is used to obtain a local free energy.The general procedure is explained in Appendix A. A gradient expansion corresponds to a Taylor expansion of the spatial Fourier transformations of these prefactors.This expansion is truncated at zeroth order, except for terms resulting from the second-order terms in Eq. ( 16) that do not involve P (we assume gradients of the polarization to be small compared to density gradients), which are truncated at fourth order.It is common in the literature to gradient expand the secondorder contribution in the functional Taylor expansion up to fourth [93] and the third-and fourth-order contributions up to zeroth order [95], and it is common in active PFC models to truncate the gradient expansion at zeroth order for the polarization terms [26].If we have m i = 0 for at least one i, the zero-wavelength component of the expansion coefficients is where R = (R 1 , . . ., R n−1 ), ϕ R = (ϕ R1 , . . ., ϕ Rn−1 ), k = ( k 1 , . . ., k n−1 ), and r = (R 1 u(ϕ R1 ), . . ., R n−1 u(ϕ Rn−1 )) are vectors (with the wavenumbers k i ).The only terms for which nonzero wavenumbers are considered are contributions of the second order in the functional Taylor expansion (16) (i.e., from c ψψ , c ψφ , and c φφ ) that do not involve P .The contributions resulting from c ψψ and c ψφ vanish when the integral over ϕ is performed if any m i = 0, and c φφ has no angular dependence anyway.Consequently, we can restrict ourselves to m = 0.All superscripts of the expansion coefficients of the direct correlation function from now on are entries of m .
Combining the Fourier and the gradient expansion and performing the angular integrals, we find The expansion coefficients A 1 -A 25 are defined in Appendix B. They are defined with a prefactor 1/n for terms of n-th order in the fields.Defining also B 1 = 2πβ −1 ρψ and B 2 = β −1 ρφ , we find the complete free energy Next, we wish to rewrite the free energy (31) in a more familiar form.In particular, we want to eliminate terms proportional to ψ 3 and φ 3 .For this purpose, we introduce the rescaled fields ψ = ψ − ∆ψ, (32) φ = φ − ∆φ (33) with the shifts Note that Eqs.(34) and (35) have their complicated form only due to the presence of third-and fourth-order contributions in the functional Taylor expansion (16).Otherwise, we would simply have The general free energy (31) can then be written in the form with the scaled shifted temperatures ˜ i , the wavenumbers qi (with i = ψ, φ, coup), and the rescaled expansion coefficients Ã3 -Ã25 .All coefficients are listed in Appendix C. Note that the coefficients A 20 -A 24 are not affected by the shifts (such that, e.g., Ã20 = A 20 ).Terms of zeroth and first order in fields that arise from the redefinition of the fields have been ignored as they do not contribute to the dynamics.We have also slightly changed the ordering of the terms compared to Eq. ( 31).The advantage of using the form ( 37) is that we can see some structure: The first terms up to Ã25 φ4 /4 are the familiar contributions from model 1 (with an additional contribution Ã19 P 4 /4) with a free energy that is an even polynomial in ψ, φ, and P .Then, there are terms coupling ψ and φ at linear order in the equation of motion, thereby including higher-order spatial gradients.We have written them in such a way that they have the same (Swift-Hohenberg-like) structure as the other linear terms (but without a prefactor 1/2, as this prefactor would not vanish after the functional derivatives).Note that ˜ coup , unlike ˜ ψ and ˜ φ , has no ideal gas contribution (see Eqs. (C18)-(C20)) and can therefore, strictly speaking, not be interpreted as a scaled shifted temperature.Finally, we have terms that lead to nonlinear couplings between the fields in the equations of motion.These arise from the ideal gas term (coupling between ψ and P ) and from higher-order contributions in Eq. ( 16) (coupling of all fields).Next, we turn to the dynamical part.For Eq. ( 13), this is rather simple: We insert the parametrization (22) into Eq.( 13) and make a constant mobility approximation, i.e., we replace the prefactor Dβ(1 + φ) by Dβ to arrive at which is identical to Eq. ( 6) up to rescaling.The situation is more complex for Eq. ( 12), where we can follow Refs.[83,84].We insert the parametrization (19), and use the fact that [83] Moreover, we make (as done in Ref. [84]) the approximation D ≈ D ⊥ = D 0 , such that D T = D 0 1.Finally, we make the constant mobility approximations for the passive terms in Eq. ( 12).Then, Eq. ( 12) yields Integrating Eq. ( 40) over ϕ and using ∂ 2 ϕ (δF/δψ) = 0, ∂ 2 ϕ u = − u, and 2π 0 dϕ u = 0 gives Similarly, multiplying Eq. ( 40) by u and integrating over ϕ using 2π 0 dϕ u ⊗ u = π1 gives Equations ( 41) and ( 38) are equivalent to Eqs. ( 6) and ( 7) in Ref. [84] and, up to rescaling, identical to Eqs. ( 1) and ( 2) in the present work.Note that the coupling term v 0 ∇ • (ψ P ) employed in Refs.[12,90] does not appear in the derivation presented here.
What is left to do now is nondimensionalization.If we write t = t 0 t, r = r 0 ˜ r, ψ = ψ 0 ψ, P = P 0 ˜ P , φ = φ 0 φ, and F = β −1 F with dimensionless time t, dimensionless position ˜ r, dimensionless free energy F , rescaled fields ψ, ˜ P , and φ, and constants and drop all tildes and underlines, we obtain the free energy and the dynamic equations Definitions of the nondimensionalized coefficients i , q i , G i , C 1 -C 7 , a 2 -a 6 , v 0 , and M φ are listed in Appendix D. Equations ( 48)-( 51) constitute the most general PFC model considered in this work, which will be referred to as model 3a.Essentially, the rescaling gives us five free parameters (r 0 , t 0 , ψ 0 , P 0 , and φ 0 ) that we can use to set five coefficients to one.We have chosen the prefactors of ψ 4 /4, φ 4 /4, and ∇ 4 ψ 2 /2 in the free energy and the mobilities of ψ and P .Thereby, our equations resemble as closely as possible those used in previous work [89].We have absorbed the prefactors 1/n into the rescaled coefficients (except in cases where the standard passive PFC model [93] also contains them).
When comparing model 3a to model 1 (Eqs.( 9)-( 11)), we can note that model 3a only contains ψ and φ, whereas Eqs. ( 9)-( 11) additionally contain parameters ψ and φ measuring the total particle number.This is a matter of notation and not a physical difference.The replacements ψ → ψ + ψ and φ → φ + φ, which are required to obtain Eqs. ( 9)-( 11), can, mathematically, be made without any problem since they simply correspond to a redefinition of the fields ψ and φ that (since we drop terms up to first order) only affects the nonlinear terms in Eq. ( 48) (of which, in model 1, only ψ 4 /4 and φ 4 /4 are left).Physically, however, there is a lot to unpack here.Often, it is argued that the parameter ψ is a measure for the total particle number of the species ψ [27,29].This, however, is in need of further justification if the parameter ψ simply arises from a redefinition of ψ that, by itself, does not necessarily have any physical meaning.
The key to understanding this aspect is the physical interpretation of the constants ρψ and ρφ in Eqs.(20) and (22).There are two options, which, for simplicity, we discuss using Eq. ( 22).First, we can choose for ρφ the homogeneous density of the liquid state.In this case, the total number of the particles of type φ is measured by the parameter ρφ , whereas the field φ only measures their spatial distribution.The advantage of this approach is that the density of the homogeneous liquid state is certainly the most natural choice for ρφ .However, the disadvantage is that, in this case, a change of the total particle number affects almost all parameters of the PFC model.In many cases, it is desirable to have only one parameter that corresponds to the total particle number, since this facilitates to study the effects of varying it by means of, e.g., numerical continuation [77].What we can make use of here is that ρφ does not have to be the actual liquid density [93].In principle -and this is the second optionwe can make any choice that is convenient (although we should ensure that it does not deviate too much from the actual density to ensure that the functional Taylor expansion ( 16) can be truncated).It has also been argued [106] that one can use ρφ to eliminate the φ 3 term from the PFC free energy functional.Namely, if one writes the free energy as with a local free energy density f , we can Taylor expand f up to fourth order around ρ φ = ρφ and then choose ρφ in such a way that the third derivative f (3) (ρ φ ) vanishes.This procedure, however, is less general than the rescaling of φ employed here since it is not guaranteed that there is a density ρφ for which f (3) (ρ φ ) = 0.For example, in the case of an ideal gas, we have f (3) (ρ φ ) = −β −1 /ρ 2 φ , which does not vanish regardless of the choice of ρφ .
If we use ρψ and ρφ as adjustable parameters, we could in principle use Eq.(D18) to set M φ = 1, namely by choosing them in such a way that We will, however, consider the more general case M φ = 1 first (as one generally has to if one uses the physical densities), and only later set M φ = 1 to obtain model 1 introduced in Section II.
Once we have fixed ρφ , the total particle number is controlled by d 2 r φ.In particular, we can substitute φ → φ + φ, where φ is chosen in such a way that d 2 r φ = 0.In this case, the parameter φ is a measure for the total particle number.The same considerations hold for the field ψ.Further rescalings of the fields ψ and φ that have been made during the derivation to simplify the free energy functional do not affect these aspects in principle, although they are, of course, relevant for the quantitative link between the parameters ψ and φ and the total particle numbers.

IV. APPROXIMATIONS AND MODEL HIERARCHY
In Section II, we have introduced the minimal model 1 on phenomenological grounds, whereas the microscopic derivation in Section III has led to the very general model 3a.In this section, we will introduce, step by step, the approximations that are required to obtain model 1.Thereby, we can get insights into their physical significance.Moreover, we obtain a hierarchy also containing intermediate models which are more general than model 1, but less general than model 3a.These arise if only some, but not all approximations are made.For ease of notation, we will use ψ and φ instead of ψ + ψ and φ + φ.
When comparing the general free energy of model 3a given by Eq. ( 48) with the minimal model 1 introduced in Section II, we can note that model 3a has six properties not present in model 1: • A: Nonlinear terms coupling ψ and P .
• B: A nonlinear term proportional to P 4 .
• D: Gradient terms in the linear coupling of ψ and φ.
• F: A mobility M φ of the field φ.
Apart from the last one, these all affect the free energy.We will now discuss the approximations that remove these features and thereby lead from model 3a to model 1, thereby also classifying intermediate models.In general, models with properties A and B will be denoted by "a", models with property B (but not A) by "b", models without properties A and B by "c", models with properties C and D by "3", models with property D but without property C by "2" and models without properties C and D by "1".For example, a model with nonlinear coupling but without nonlinearities in P would be "model 3c".This naming scheme is illustrated in Table I.Combinations that do not fit into this classification (such as a model that has a nonlinear coupling of ψ and φ but no gradients in the linear coupling) will not be considered, since their derivation would require keeping some nonstandard third-and fourth-order terms while dropping some standard second-order ones.Notably, properties A-F all concern the passive part of the transport equations and not the nonvariational active term.Therefore, our discussion is also relevant for passive PFC models and will be made with a particular emphasis on the correct passive limit.This gives our results additional significance for applications in materials science (for example in models for passive liquid crystals doped with spherical particles [109,110]).Properties E and F do not give rise to additional terms, such that we will not use them to specify an additional class of models (apart from the fact that we reserve the name "model 1" for the case First, we consider the step from models of type 3 to models of type 2. For this purpose, we make the Ramakrishnan-Yussouff approximation [111], which is very common in the derivation of PFC models from (D)DFT and typically made right from the beginning [93].In this approximation, the functional Taylor expansion ( 16) is already truncated at second order.In other words, the expansion coefficients A 11 -A 25 are all set to zero.The general free energy (48) then simplifies to Moreover, as can be seen from the microscopic definitions of the coefficients listed in Appendices B-D, the microscopic expressions for the expansion coefficients simplify drastically.For example, the coefficients C 2 -C 4 now only arise from the ideal gas free energy.The free energy (54) constitutes model 2a.It is a model of type 2 by the classification introduced above as it has property D but not property C, and it is a model of type a as it has properties A and B. The main differences to model 3a are • There is no nonlinear coupling between the fields ψ and φ (previously encoded in the coefficients a 2 -a 6 of model 3a).
• There is no coupling between P and φ.
This means that both corresponding aspects of model 3a only arise from higher-order contributions in the functional Taylor expansion, and can thus not be obtained within the Ramakrishnan-Yussouff approximation that is usually considered.This shows three advantages of our approach: First, by considering these higher-order terms, we were able to show that there can be a direct coupling of P and φ in a binary mixture of active and passive particles if φ is the passive field, namely through higher-order effects.Second, we show that models with coupling of P and φ but without nonlinear coupling of ψ and φ are not reasonable since these couplings have the same physical origin (higher-order terms in Eq. ( 16)).Third and most importantly, the assumption made in model 1 that fields describing two different particle species couple only linearly (through quadratic terms in the free energy) is quite robust, since nonlinear terms only arise if we go beyond the approximations PFC models are usually based on.This also justifies the phenomenological coupling terms used in simpler models such as the one from Ref. [77].
We now move further from models of type 2 to models of type 1.For this purpose, we have to get rid of the higher-order terms in the gradient expansion.To recover model 1, we have to set A 9 and A 10 to zero, while keeping all other terms that are of second order in the functional Taylor expansion (16).A "naive" argument for dropping these two coefficients while keeping all the others can be found from an expansion in a smallness parameter.Let us assume that the fields are slowly varying in space, such that terms of order n in gradients are of the order ε n gr , where ε gr is a small dimensionless parameter.Moreover, we assume that the ratio of the various correlation functions is given by where ε cor is another small dimensionless parameter.These parameters now have to be tuned in such a way Level of approximation With ψ-P coupling (a) With P 4 term (b) Without P 4 term (c) With nonlinear interaction terms (3) Model 3a (Eq.( 48)) Model 3b (Eq.( 67)) Model 3c (Eq.( 70)) Ramakrishnan-Yussouff approximation (2) Model 2a (Eq.( 54)) Model 2b (Eq.( 68)) Model 2c (Eq.( 71)) Coupling without gradients (1) Model 1a (Eq.( 59)) Model 1b (Eq.( 69)) Model 1c (Eq.( 72)) that terms proportional to ψ ∇ 4 ψ (of order ε 4 gr ) and φψ (of order ε cor ) are kept in the free energy, whereas terms proportional to ψ ∇ 6 ψ (of order ε 6 gr ) and φ ∇2 ψ (of order ε 2 gr ε cor ) are dropped.This can be achieved considering the distinctive limit characterized by the condition To elucidate the physical meaning of the smallness parameter ε cor , we consider the random phase approximation [93] where U ij is the interaction potential between particles of species i and j.Thus, in the simplest case, the direct correlation function depends on the strength of the interaction, i.e., the assumption (55) implies that the interaction of the particles of different species is weaker than between particles of the same species.This would imply that model 1 describes a mixture of active and passive particles where the interaction between different particle types is weaker than the interactions of the particles of one species among themselves.The reason for why this argument is somewhat naive is that the quantitative predictions obtained for the PFC model parameters by a derivation from DFT can be quite inaccurate, since the assumptions made in this derivation (ψ and φ are slowly varying in space) are not really justified in the case of crystallization [112,113].Nevertheless, one can still learn something from such derivations, since the qualitative predictions and the mathematical structure of the PFC models so obtained can still be quite accurate [95], and it is the qualitative structure that we are interested in in the present work.However, it should be noted that in practice, one typically obtains the coefficients of a PFC model by fitting a fourth-order polynomial to the Fourier-transformed direct correlation function ĉ( k) [112] depending on the wavenumber k.Consequently, dropping the gradients in the linear coupling is a good approximation as long as it is a good approximation to replace ĉψφ ( k) by a constant while ĉψψ ( k) and ĉφφ ( k) are fitted with fourth-order polynomials.This, however, is consistent with the basic conclusion of our "naive" argument, namely that we can drop the gradient terms if c ψφ is small compared to the other correlation functions.The reason for this is that, if c ψφ is small, any errors we make in fitting this function (such as replacing it by a constant) are less relevant for the overall accuracy of the free energy functional we derive.
The direct correlation functions (more precisely: the static structure factors S which can be calculated from them) are required as an input not only in PFC models, but also in the mode coupling theory of the glass transition [114].A mode-coupling theory for a mixture of active and passive particles has recently been developed by Feng and Hou [13], who obtained the structure factors via Brownian dynamics simulations.They found that, while the overall shape of the partial structure factor S ψφ (k) (k = k ) is similar to that of the partial structure factors S ψψ (k) and S φφ (k), S ψφ (k) is smaller (around zero and even negative for some k).
The Fourier-transformed direct correlation functions are related to the structure factors by [13] Equation (58) shows that, if S ψφ is small, then (other things being equal) ĉψφ will also be small.Consequently, the simulation results from Ref. [13] support the assumption made in model 1 that gradient terms in the linear coupling are less important than the other gradient terms.
Dropping the gradient terms in the coupling in Eq. ( 54) thus gives the free energy of model 1a, i.e., with the linear coupling parameter Model 1a still has nonlinear terms proportional to ψ P 2 , ψ 2 P 2 , and P 4 that are not present in model 1.Compared to the interaction-dependent terms, it is more difficult to see how one can eliminate them given that they arise also from the ideal gas free energy (23).Thus, it looks as if they are present regardless of the assumptions we make about interactions.However, the situation is more intricate.If we make the assumption that the interaction does not depend on the particle orientation as it is the case for Brownian spheres, in the passive limit, the free energy should not depend on P .In the active case, we can add a term from the ideal gas term).This is unproblematic for the passive limit since, for v 0 = 0, P does not couple to ψ.Thus, to ensure the correct passive limit, we should (depending on other choices we make in modeling the interactions) use model 3c, 2c, or 1c.
What is confusing here is that the limit v 0 → 0 based on models of type a does not give us a model of type c even for orientation-independent interactions, since the terms nonlinear in P and, in particular, the terms coupling ψ and P are still present.In the ideal gas limit, passive models of type a give Note that taking the ideal gas limit is problematic in the rescaled model also because we have forced the prefactor of the term ψ ∇ 4 ψ/2 in the free energy, which should be zero in the ideal gas limit, to be one.We have simply dropped this term in Eqs. ( 61) and ( 62), but we now will turn to the model given by Eqs. ( 38), (41), and (42) (i.e., the "original" model without rescaling and nondimensionalization), to discuss the ideal gas limit.Equations ( 61) and (62) show that there is still a coupling between P and ψ (which should not be the case for spheres in the passive limit), and also some nonlinear terms which are obviously unphysical.Physically, however, we should simply have a linear diffusion equation for ψ.This is also the result that we would get from DDFT.To see the origin of this problem, it is helpful to take one step back and consider the single-component PFC model for passive spheres.This was discussed in much detail by Archer et al. [95], who identified the Taylor expansion of the logarithm and the constant mobility approximation (CMA) as crucial steps that lead to potential unphysical behavior.Here, we will briefly discuss the implications of this issue for the dynamics.For the passive field φ, the model (38) with CMA and free energy (23) gives The result ( 63) is not correct since, for an ideal gas, the dynamics of φ is simply given by the standard diffusion equation.If we use the nonconstant mobility 1 + ψ obtained from DDFT instead, we find and dropping terms of order φ 4 .This is the correct description of the dynamics.If, however, we make a CMA, the only way to recover the diffusion equation for noninteracting particles is to set Similarly, if we describe the field ψ with CMA using Eq. ( 41), insert the ideal gas free energy (23), and set v = 0, we find Without a CMA one finds 2 where in the last step we have ignored terms of fourth order in the fields.Here, terms coupling ψ and P disappear as they should.
One might now ask why we do not just stop the Taylor expansion of the logarithm after the quadratic term, in which case the PFC model gives the correct equation of motion for the ideal gas limit (i.e., the diffusion equation).The reason is that we also have to ensure the correct static results in the passive limit: For the passive PFC model, as for the passive DDFT, the system evolves toward the state that minimizes the free energy functional.Consequently, whether the phase transitions predicted by the model are correct depends on the accuracy of the free energy functional, which increases if we take into account higher orders in the expansion of the logarithm.
A particularly important reason why the higher-order terms can be important is stabilization.As an example, let us assume that we can treat interactions in the Ramakrishnan-Yussouff approximation.In this case, the coefficient C 1 has contributions from the ideal gas and excess free energy, while the coefficient C 2 arises solely from the ideal gas free energy.Consequently, although we certainly have C 2 > 0, we might have C 1 < 0 if the interactions are sufficiently strong.In this case, if we drop higher-order terms in the free energy, the polarization would grow without bounds, which is obviously unphysical.This can be prevented by introducing a term 1  4 C 2 P 4 , resulting from the ideal gas free energy.In this case, instead of an (unphysical) blow-up, we get the well known phenomenon of self-polarization [27], which, as is evident from the microscopic derivation, is only possible for 2 To get the first equality of Eq. ( 66), replace D 0 ∇ 2 by D 0 ∇ • (1 + ψ + P • u) ∇ in Eq. ( 40) (with v = 0) and integrate over ϕ.The term ( ∇δF/δ P ) • P is, in index notation with summation convention, given by (∂ i δF/δP j )P j .
orientation-dependent interactions.Consequently, if the terms in Eq. ( 48) that are quartic in the free energy become negative due to interactions, it might be necessary to include even more terms in the ideal gas expansion.
For this particular purpose, however, it is sufficient to keep the term 1  4 C 2 P 4 in addition to 1 2 C 1 P 2 .Terms coupling ψ and P , which (in the Ramakrishnan-Yussouff approximation) only arise from the ideal gas free energy, are almost always ignored in PFC models.While this is typically done without any justification, we can now give a physical motivation for such models: They are a compromise between models of type a, where all nonlinear terms in P are kept (giving a stable and accurate free energy functional and an incorrect diffusion equation), and models of type c, which are appropriate for active Brownian spheres.Since they are in between these two extremes, we call them "models of type b".Type b models cannot be obtained via a systematic expansion in a smallness parameter, they rather arise as a minimal model allowing for orientation-dependent interactions.They also come in three forms, with the most complicated one being model 3b: Terms coupling φ and P have also been dropped (since it is somewhat inconsistent to keep them if we drop terms coupling ψ and P ).A simpler variant, obtained using a Ramakrishnan-Yussouff approximation, is model 2b: Finally, by dropping gradients, we can also get model 1b: + a 1 ψφ + 1 2 Let us now assume that there are no orientationdependent interactions at all.In this case, C 1 is guaranteed to be positive, and there is no need for terms of higher order in P to ensure stabilization.On the other hand, as discussed above, these terms give unphysical results in the ideal gas limit.If we, based on these considerations, drop all terms involving P apart from C 2 P 2 in the free energy (48) of model 3a as it is appropriate for active Brownian spheres, we are left with model 3c: Also making a Ramakrishnan-Yussouff approximation then gives model 2c: Finally, if we also drop gradient terms in the interaction, we get the simple model 1c: Model 1 introduced in Section II is identical to model 1c apart from the fact that G φ = 1 and M φ = 1 in model 1b.In principle, model 1 is "only" a special case of model 1c recovered by setting the coefficients G φ and M φ to one.Recall that we can get M φ = 1 also by adjusting the reference densities (see Eq. ( 53)).Setting M φ = 1 and G φ = 1 and dropping the subscript "1" from a 1 , Eqs. ( 49)-( 51) and (72) give + aψφ + 1 2 To summarize: Model 1 is a PFC model for a mixture of active and passive particles, which assumes that • nonlinearities in P resulting from the ideal gas free energy can be ignored (generalization: models 3a, 3b, 2a, 2b, 1a, and 1b); • the excess free energy functional can be obtained from the Ramakrishnan-Yussouff approximation (generalization: models 3a-3c); • gradient terms in the coupling can be ignored (generalization: models 3a-3c and 2a-2c); • the mobilities of both fields are equal; • the factor G φ can be set to one (i.e., the prefactor of the gradient terms is equal for both fields).
We have thus achieved a systematic derivation of the minimal model 1 and identified the approximations required for obtaining it.Moreover, we have derived and classified more general models (3a-3c, 2a-2c, and 1a-1c), which provide a very general description of mixtures of active and passive particles.An overview over the model hierarchy is given in Table I.

V. LINEAR STABILITY ANALYSIS
Next, we consider the linear stability of uniform states for a mixture of active and passive particles using the above derived model 1 as given by Eqs. ( 9)- (11).The limiting case of a mixture of passive particles can easily be obtained by setting v 0 = 0. From now on we restrict ourselves to one spatial dimension.
The real parts Re(λ) of dispersion relations λ j (k) with j = 1, 2, 3 are displayed in Fig. 1 for four exemplary parameter choices.In Fig. 1(a), a case of low activity v 0 = 0.05 is given where a monotonic small-scale instability occurs, i.e., above instability onset, where a single mode of wavenumber k c = 0 (the critical value) becomes unstable, a finite band of unstable wavenumbers exists centered about k c .In a direct time simulation this produces a resting crystal (not shown).A decrease in the effective temperature ψ will widen the band of unstable wavenumbers.A case of higher activity v 0 = 0.4 is presented in Fig. 1(b).Again, a small-scale instability arises (finite band of unstable k centered about k c = 0), however, the instability is now oscillatory.A time simulation produces a traveling periodic state, i.e., a traveling crystal (not shown).Both cases presented in Figs.1(a) and 1(b) are also encountered in one-component active PFC models [27].This does not apply to the remaining cases of Fig. 1.In particular, Fig. 1(c) shows a case of small activity v 0 = 0.05 where two monotonic small-scale modes are unstable.This is also found in the limiting case of a PFC model for a mixture of passive particles [77].Increasing the activity to v 0 = 0.4 renders one of these unstable modes oscillatory, as presented in Fig. 1(d).
Note that an increase in the linear coupling between the densities ψ and φ generally results in a further destabilization.This can already be anticipated inspecting the dispersion relation in the passive limit that reads [77] There, increasing the coupling strength a, always results in an increase of the largest eigenvalue λ + .

VI. NONLINEAR STATES
In this section, we present selected fully nonlinear states as obtained by numerical path continuations [115,116], as well as results of direct time simulations.Thereby, the focus is on effects that are not seen in the standard one-component active PFC model.For the path continuations we employ the package pde2path [117] while all time simulations are performed using a semi-implicit Euler scheme with adaptive step size.All results from time simulations are presented after initial transients have decayed.In the following, we only consider the particular case of equal temperature parameters ψ = φ , and equal critical wavenumbers q ψ = q φ = 1 for the active and passive particles, and drop the subscripts accordingly.
First, we determine resulting states in the passive limit v 0 = 0.The corresponding phase diagram in Fig. 2 is quite similar to the one presented in Ref. [77].Pairs of thermodynamically stable coexisting phases lie on the black binodal lines.Note that, for clarity, in Fig. 2 we only present thermodynamically stable coexistence and omit metastable and linearly unstable coexistence.Particular coexisting states with equal chemical potentials and pressure on two binodals are connected by gray tie lines.Within the hatched region, uniform states are unstable w.r.t.phase separation resulting in coexistence of the two bordering phases.The two gray triangles mark FIG. 2. Phase diagram of the PFC model for a mixture of passive particles in 1D displayed in the ( ψ, φ)-plane at = −1.5 and v0 = 0. Four phases (liquid, alloy, ψ-crystal, φcrystal) can be distinguished.The binodal lines are marked in black.Coexisting states on the binodals are connected by gray tie lines.The thus hatched areas correspond to two-phase coexistence of adjacent phases.Three-phase coexistence is indicated by triangular gray shaded areas.The binodals and tie lines are only presented for thermodynamic coexistence, not for metastable or linearly unstable coexistence.Along the horizontal blue line, the bifurcation diagrams in Figs. 3,  5, and 6 are presented, while the diagonal blue line denotes the path corresponding to Fig. 7.The remaining parameters are as in Fig. 1.
regions where phase separation into three phases occurs.Overall, Fig. 2 shows four thermodynamically stable liquid (uniform) and crystal (periodic) phases.At low densities ψ and φ, the liquid state is the globally stable state.In the top left and bottom right corners of the shown range, one of ψ and φ is high while the other one is low, resulting in a crystal state of the particles with the high density, while the other particles are in a weakly modulated liquid state.When ψ and φ are both large, the resulting state is a crystalline alloy: both fields show peaks of similar amplitudes that are in phase.For a deeper discussion of such phase diagrams, see Ref. [77].Using the phase diagram in Fig. 2 for the passive mixture as reference, we analyze the bifurcation behavior in the active case: Fixing the activity at v 0 = 0.23, 0.3, or 0.4, we consider two paths through parameter space.First, we vary ψ at fixed φ = 0 (Figs.3-6) and second we consider the path defined by φ = ψ (Fig. 7).
All resting states w = (ψ, P, φ) T are characterized by  II.The insets magnify regions where the branch of resting periodic states changes stability.Hopf bifurcations occurring on this branch are marked by diamonds.Figure 4 shows selected space-time plots at the locations marked by letters "a"-"c".The corresponding symbols are enlarged and red.The remaining parameters are = −1.5, q = 1, a = −0.2,Dr = 0.5, and C1 = 0.1.The domain size is fixed at L = 16π.
their L 2 norm with the domain length L. Traveling states are characterized by the temporal average of w .At small activity v 0 , the overall structure of the bifurcation diagrams is very similar to the one in the passive limit (see, e.g., Figs. 8 and 9 of Ref. [89]).With increasing activity, the parameter range where localized states (LSs) exist slightly decreases [89].
Here, we first consider states at fixed mean density φ = 0 where the field φ is in the periodic crystalline state.That is, it provides a periodic background on which, for increasing mean density ψ, LSs, i.e., crystallites in ψ, grow.Although the periodic background influences where the active particles prefer to be, it is not entirely static as it is also influenced by the active particles.The particular results for v 0 = 0.23 are summarized in Fig. 3.At low densities ψ, resting periodic states exist.These  are destabilized at ψ ≈ −0.753 in a Hopf bifurcation, marked by the leftmost diamond symbol in the upper left inset of Fig. 3.After further destabilization by another Hopf bifurcation at ψ ≈ −0.750, a subcritical pitchfork bifurcation occurs at ψ ≈ −0.735.There, two branches of symmetric LSs emerge with respective odd and even number of peaks.Both branches of symmetric LSs are initially unstable.On each branch a series of Hopf bifurcations occurs before the first saddle-node bifurcation (not shown) rendering the LS less unstable.The branches of symmetric LSs are interconnected by ladder branches of asymmetric LSs that emerge in pitchfork bifurcations.As one moves along the branches of symmetric LSs consecutively, pairs of additional peaks are added.Note that, in contrast to the behavior in the passive limit at the same temperature (not shown), this does not involve saddle-node bifurcations.Similar slanted snaking behavior is observed at higher effective temperature in the passive limit of the one-component PFC model (see, e.g., Refs.[89,97]).Between ψ = −0.743and ψ = −0.537,both branches of symmetric LSs undergo Hopf and pitchfork bifurcations.As a result, they change stability in such a way that always at least one of them is linearly stable.
When the whole domain has filled with density peaks, the branches of symmetric LSs reconnect at ψ = −0.564 to the branch of periodic states.This branch regains stability at ψ = −0.512after undergoing three further Hopf bifurcations, visible in the lower right inset of Fig. 3. Further following the branch, at ψ = −0.234 it becomes unstable again in a drift pitchfork bifurcation where a branch of steadily traveling periodic states emerges.The panels in Fig. 4 show space-time plots of the density of the active species ψ and the passive species φ after transients have decayed.In particular, we show three states that have not yet been described and may not exist for the one-component active PFC model [27,29].At φ = 0, the passive background φ(x, t) is periodic, in the considered domain this corresponds to a periodic solution with n = 8 peaks.Figures 4(a Figure 4(a) shows an oscillating LS in ψ on a periodic background in φ.Thereby, each individual peak oscillates like a standing wave, with neighboring peaks oscillating in anti-phase.With other words, active particles alternately co-occupy neighboring sites where passive particles are located.The behavior is strongest at the center of the LS while oscillation amplitudes become smaller as one moves toward the outside tails of the LS.Overall, the structure resembles a spatially localized space-time checkerboard pattern, somewhat similar to the modulated standing waves in Ref. [118].We call it an alternating LS.
Figure 4(b) shows an alternating LS at higher ψ that accordingly has a larger spatial extension than the one in Fig. 4(a) and nearly fills the domain.Another consequence is the much smaller amplitude of the oscillations at the center of the LS as compared to the larger amplitude at the fringes.With other words, the three central sites of the passive background are always co-occupied by active particles, forming a nearly steady core, while the sites further from the center show oscillating occupancy by active particles.Note that in Fig. 3 one also finds fully spatially periodic (not localized) variants of such an alternating state.They are marked by filled squares, and a typical space-time plot is given in Fig. 5(e) (see below).
Figure 4(c) presents a qualitatively different state that features a slowly traveling passive background ψ(x, t).At the same time, the LS in φ(x, t) alternately oscillates and travels at a faster speed into the opposite direction than the passive background.Similar to Fig. 4(b), the (now traveling) LS almost fills the domain and the amplitude of the oscillations is much smaller at the center than at the edges of the LS.The LS travels by gaining density on one side and losing it on the other.The sites further away from the center of the LS again show alternately oscillating occupancy.
The various described patterns dominate in different regions of the snaking structure in Fig. 3. Thereby, time simulations result in steady LSs in the central region in the range −0.698 ψ −0.581, while closer toward the ends of the snaking structure oscillatory states dominate.This can be seen in the magnifications shown as insets in Fig. 3.
Results for a slightly higher activity v 0 = 0.3 are presented in Fig. 5.For this value, all branches of steady LSs have vanished, i.e., we are in a region analogous to the one above the final cusp in Fig. 12 of Ref. [89].The periodic steady state is linearly stable at low densities as before.It becomes unstable via a Hopf bifurcation at ψ ≈ −0.717.Furthermore, at ψ ≈ −0.512 a branch of traveling periodic states emerges in a supercritical drift pitchfork bifurcation.It is initially unstable but eventually stabilized through a series of Hopf bifurcations with the final one occurring at ψ ≈ −0.434.Employing numerical path continuation of steady states, here, we only access resting and steadily traveling states.In the ψrange where all of these states are unstable, we determine the system behavior by direct time simulations.As a result, we find five types of standing, traveling, and modulated traveling wave states, exemplified in the space-time plots in the insets of Fig. 5.
At high densities, two kinds of traveling states exist and coexist: Inset 5(a) shows a steadily traveling periodic state, also obtained by continuation.The second kind, presented in inset 5(b), consists of oscillating directionreversing periodic states.The wiggling back and forth motion of density peaks is overlaid by a slow drift of the entire pattern.This is akin to the oscillating directionreversing LSs found in Ref. [29] (see their Fig.11(a)), albeit there they do not show an additional drift.Most likely, the corresponding branch emerges from the branch of steadily drifting states via a Hopf bifurcation.Note that there is a range of multistability.Another type of wiggling state is presented in inset 5(c).There, two patches of wiggling peaks (having three peaks each and wiggling in anti-phase) are separated at the center and at the (periodic) boundary by a respective localized oscillating one-peak state.Thereby, the two one-peak states communicate via the wiggling patches such that together they form an alternating space-time pattern, similar to neighboring peaks in Fig. 3.At first sight, the space-time plot in inset 5(d) seems to show a similar state.However, there, the fast interrelated oscillating and wiggling pattern is furthermore modulated by a slower oscillation, i.e., there is a second, lower frequency, i.e., a Hopf bifurcation is the most likely transition scenario.Finally, inset 5(e) represents a fully space-time periodic checkerboard pattern where each peak oscillates and neighbors alternate in their oscillation.
Note that for all time-periodic states shown in the insets 5(a)-(d), the density of the passive particles φ normally closely follows the one of the active species ψ.An exception are nearly static density peaks in φ whose sites are co-occupied by the oscillating peaks in ψ in insets 5(c)-(e).
Next, in Fig. 6 we present the results at the highest here considered activity, v 0 = 0.4.As in Fig. 5, branches of resting and traveling periodic states are obtained by numerical continuation.The resting states are linearly stable at low densities ψ, and become unstable via a Hopf bifurcation at ψ ≈ −0.656.Twelve further Hopf bifurcations occur along the branch (not shown).Then, at ψ ≈ −0.486 a branch of unstable steadily traveling periodic states emerges in a subcritical drift pitchfork bifurcation.In contrast to the case of Fig. 5, this branch remains unstable.Again, this leaves a large ψrange where all steady states are unstable and we resort to time simulations.At low densities, the time simulations produce the expected steady states.Directly beyond the first Hopf bifurcation, alternating periodic states are found similar to the ones presented in Fig. 5(e).Above ψ = −0.6,two kinds of traveling states are found not known from the cases of lower activity.Both have in common that the passive background performs a wiggling back-forth motion, overlaid by a slow drift of the whole pattern, similar to Fig. 5(b).The active particles, however, move at a much larger speed in the same direction.The wiggling motion of the background is then due to a time-periodic pull the active particles (ψ) exercise onto the passive particles (φ).At densities between ψ ≈ −0.6 and ψ ≈ −0.44 the active particles (ψ) organize into localized density peaks.Inset 6(a) shows such a state with three localized density peaks of ψ (in red) that repeatedly move through the domain, and the background of passive particles (φ, in blue) moving much slower.The number of peaks in ψ increases with ψ, until the whole domain is filled with density peaks.An example of such a traveling periodic state with different speeds of the active and passive particles is presented in inset 6(b).For both, localized and periodic traveling states with two speeds, the speeds increase with increasing mean density ψ.
Finally, in Fig. 7 we return to the lower activity of v 0 = 0.3 and analyze states occurring along the line through the phase diagram (Fig. 2) defined by ψ = φ.Contrary to the previous case where we had fixed φ = 0, now the passive species no longer provides a periodic background.Instead, it shows similar density profiles as the active species.At low densities ψ = φ, the uniform liquid state is linearly stable.It is destabilized in a supercritical pitchfork bifurcation at ψ = φ ≈ −0.721,where the branch of resting periodic states with n = 8 peaks emerges.Close to onset, these states are linearly stable, but very soon loose stability in a secondary pitchfork bifurcation at ψ = φ ≈ −0.72.There, two branches of symmetric resting LSs emerge subcritically, with an odd and even number of peaks, respectively.Both branches are initially unstable.The branch of odd LSs gains stability in a saddle-node bifurcation where it folds back toward higher densities.The branch of even LSs gains stability after a saddle node bifurcation and a pitchfork bifurcation, where the first of the ladder branches of asymmetric LSs emerges.Similar to the passive limit, the two branches of symmetric LSs undergo slanted snaking involving a series of saddle-node and pitchfork bifurcations.The latter produce in total five branches of asymmetric LSs.The branches of resting symmetric LSs reconnect to the branch of resting periodic states at ψ = φ ≈ −0.486 when the whole domain is filled with peaks.In Figs. 3,  5, and 6, all LSs consist of crystalline patches of the active species (ψ) on a periodic background of the passive species (φ).Here, all LSs are crystalline patches of an alloy of the passive and active species in a liquid background.
In contrast to the passive limit, at the chosen activity only the sub-branches of symmetric LSs with one, two, and three peaks are linearly stable, as further up each branch undergoes more than 30 Hopf bifurcations.The highest density for which resting LSs exist is ψ = φ ≈ −0.741.Also the steady periodic states undergo further bifurcations: at ψ = φ ≈ −0.598 a branch of steadily traveling periodic states emerges in a supercritical drift pitchfork bifurcation.Initially it is unstable, but after undergoing a total of ten Hopf bifurcations it gains stability in a final one at ψ = φ ≈ −0.485.There is also a period-doubling pitchfork bifurcation on the branch of steadily traveling periodic states, where a branch of drifting states emerges, that are never stable and not shown here.
As before, this leaves a density range, where no resting or steadily traveling states exist.Resorting to time simulations, at low densities up to ψ = φ = −0.76we recover first steady periodic states and then steady odd LSs.At higher densities, there are alternating LSs (not shown here), similar to those in Figs. 3 and 4, albeit on a liquid background.At a slightly higher mean density, inset 7(a) shows a state combining two localized patches of different states encountered before: At the center of the domain is a single wiggling peak performing a back-forth motion without net drift, while on the (periodic) boundary a single peak of the active species oscillates, i.e., behaves like a standing wave.Similar to the behavior of their domain filling variants, seen in Fig. 5, the passive density peak at the center closely follows the wiggling motion, while the passive density peak on the boundary oscillates with a very low amplitude only.
Another combination of two localized patches of alternating and wiggling LSs is presented in inset 7(c).There, a patch of two wiggling peaks coexists with a patch of an alternating LS with three peaks.This is furthermore overlaid by a very slow drift of the entire pattern.Note that the two outer peaks of the alternating LS are not at an equal distance to the wiggling LS.Because of this, both patches of LS are asymmetric, in contrast to those in inset 7(a) that show a distinct spatio-temporal symmetry.Inset 7(b) shows an example of the irregular, potentially chaotic states shown as star symbols in the main panel in Fig. 7. Here, we do not analyze these states further.Lastly, inset 7(d) gives a steadily traveling state consisting of a combination of two different LSs, with two and four peaks, respectively.The whole pattern moves at a constant speed with the passive density closely following the active one.Similar to the passive limit, at high densities the entire domain is filled with density peaks.There, only steadily traveling periodic states are found.

VII. CONCLUSIONS
We have discussed a hierarchy of active PFC models that represent continuum descriptions for a mixture of active and passive colloidal particles at different levels of approximation.First, we have presented a systematic microscopic derivation of a very general PFC model (model 3a) from a DDFT, i.e., a microscopic continuum description that itself may be derived from a microscopic particle-based description [53].The derivation we have presented here includes a systematic treatment of the relevant orientational degrees of freedom.Thereby, our particular interest has been on the establishment of the nonlinear and coupling terms.Then, we have employed a series of approximations that have a clear physical meaning and motivation to simplify the general PFC model.Passing through these approximation steps, a hierarchy of models (models 3a-1c) has been established that allows for interesting insights into the microscopic justification of constructions used in various PFC models for active particles and mixtures, the approximations required for obtaining them, and possible generalizations.The minimal model (model 1) indeed corresponds to the one constructed and analyzed in Ref. [89] by purely phenomenological means.
Note that, in passing, the presented derivation has also recovered a number of related models and contributed to the understanding of the employed approximations and their limitations.Namely, when eliminating the passive species from model 1 one obtains the active PFC model for a single species of active colloidal particles derived and analyzed in Refs.[26-29, 84, 91].Furthermore, a PFC model for a mixture of passive particles with a linear coupling between the two fields [77] (also see Section 4.1 in Ref. [89]) is recovered when taking model 1 in the limit of vanishing activity.Note that this has allowed us to use the phase diagram obtained in Ref. [77] as a reference for the analysis of the parameter space.This implies that the various models 2 and 3 may in the passive limit be used as more exact models for mixtures of passive particles, the most precise being model 3a with v 0 = 0. (Note, however, that a coupling to P should not be present in the passive limit for spherical particles, see Section IV.)In this way one may, e.g., obtain a model closely related to the one derived in Ref. [76].There, however, four-point correlations were not included.In contrast, the present model 3a contains entropic and enthalpic fourth-order terms in the free energy functional.
In the second part of this work, we have restricted our attention to the derived minimal model 1.Linear stability analysis of the trivial uniform state at different mean concentrations has shown that it may become unstable to either a monotonic or an oscillatory small-scale instability similar to the one-component active PFC model [27].The difference is that now there can be two instability modes at similar or different wavenumbers active at the same time.These may be two monotonic modes or a monotonic and an oscillatory mode.
Beside the linear considerations, we have employed numerical continuation and time simulation methods to investigate the fully nonlinear regime.Thereby, the use of the former has allowed us to determine branches of periodic and localized steady and steadily traveling states.These can form the intricate intertwined slanted snaking structures expected for Swift-Hohenberg-type systems [119,120] with conservation laws [89,97,121].However, in contrast to classical (conserved) Swift-Hohenberg-type systems, that are normally variational and therefore do only allow for steady states, here activity supports several types of standing, traveling, and modulated periodic and localized wave patterns.For instance, we have described direction reversing traveling periodic and localized states that we have called "wiggling states" as in Ref. [29] (where also other names and references for these states are given).Another state are spatio-temporal patterns where individual peaks behave like a standing wave and neighboring peaks show alternating (or anti-phase) oscillations.This resembles the modulated standing waves described in Ref. [118] and has to our knowledge not yet been described for PFC systems.Here we have encountered such states in spatially periodic and localized versions.A recent study of nonreciprocally coupled Cahn-Hilliard equations reports on localized states whose outer peaks oscillate asymmetrically or symmetrically either in phase or in anti-phase [122].Furthermore, the different regular spatio-temporal patterns can coexist in different regions of space, thereby giving rise to more intricate behavior.Transitions to seemingly irregular behavior have also been observed.
A systematic investigation of the transitions between the various time-dependent states has been outside the scope of the present work, but forms a formidable future challenge.Further possible extensions of this work include detailed comparative investigations of the more general models.An important remaining question is whether the minimal model 1 is already able to describe all states occurring in model 3a including their stability and sequence of appearance.Moreover, one could extend the derivation by incorporating also the nematic order parameter (as done in previous derivations for passive models [80][81][82][83]), particle inertia (as done for a singlespecies model in Refs.[87,88]), or by considering a mixture of two different active species.

DATA AVAILABILITY
The raw data corresponding to the figures shown in this article are available as Supplementary Material [123].

FIG. 3 .
FIG. 3. Bifurcation diagram showing branches of onedimensional states of the PFC model 1 (9)-(11) for coupled active and passive particles.Shown is the L 2 norm w as a function of the mean density ψ at fixed φ = 0 and activity v0 = 0.23.Solid [dashed] lines indicate stable [unstable] states.The blue line corresponds to crystalline (periodic) states with n = 8 peaks.The intertwined light and dark purple lines represent the slanted snaking of branches of LSs with odd and even number of peaks, respectively.They are interconnected by branches of asymmetric states (black lines).The branch of steadily traveling periodic states is given as orange line.The various states obtained by time simulations are indicated by symbols according to TableII.The insets magnify regions where the branch of resting periodic states changes stability.Hopf bifurcations occurring on this branch are marked by diamonds.Figure4shows selected space-time plots at the locations marked by letters "a"-"c".The corresponding symbols are enlarged and red.The remaining parameters are = −1.5, q = 1, a = −0.2,Dr = 0.5, and C1 = 0.1.The domain size is fixed at L = 16π.

FIG. 4 .
FIG. 4. Panels (a) to (c) show space-time plots of the densities of active (ψ(x, t)) and passive (φ(x, t)) particles at the locations indicated by corresponding letters in Fig. 3 (after transients have decayed).Panels (a) and (b) show alternating LSs in the interval t = [0, 500].Due to the slower dynamics, the interval in panel (c) is t = [0, 1000], showing alternating traveling LSs.
) and 4(b) show time-periodic states where the passive background is nearly steady: it only shows small oscillations as reaction to the large amplitude oscillations in the LS in ψ.

FIG. 5 .
FIG.5.Bifurcation diagram as a function of the mean density ψ at fixed activity v0 = 0.3.The model, solution measure, remaining parameters, symbols, and line styles are as in Fig.3, while the insets show space-time plots as in Fig.4(a), focusing on the density of the active particles ψ.

FIG. 6 .
FIG.6.Bifurcation diagram as a function of the mean density ψ at fixed activity v0 = 0.4.The two insets show space-time plots of the active and passive densities (upper panel: ψ(x, t); lower panel: φ(x, t)) at the marked locations in an interval t = [0, 500] after initial transients have decayed.The model, solution measure, remaining parameters, symbols, and line styles are as in Fig.3, while the insets are as in Fig.4(a).

FIG. 7 .
FIG. 7. Typical bifurcation diagram as a function of the mean densities ψ = φ at fixed activity v0 = 0.3.The light blue line corresponds to liquid (uniform) states.The model, solution measure, remaining parameters, symbols, and line styles are as in Fig. 3, while the insets are as in Fig. 4(a), focusing on the density of the active particles ψ.

TABLE I .
Overview over the various models.

TABLE II .
The various observed states and the symbols by which they are marked in the bifurcation diagrams in Figs. 3 and 5-7.