Chromo-dynamic multi-component lattice Boltzmann equation scheme for axial symmetry

We validate the chromo-dynamic multi-component lattice Boltzmann equation (MCLBE) simulation for immiscible fluids with a density contrast against analytical results for complex flow geometries, with particular emphasis on the fundamentals of the method, i.e. compliance with inter-facial boundary conditions of continuum hydrodynamics. To achieve the necessary regimes for the chosen validations, we develop, from a three-dimensional, axially-symmetric flow formulation, a novel, two-dimensional, pseudo Cartesian, MCLBE scheme. This requires the inclusion in lattice Boltzmann methodology of a continuously distributed source and a velocity-dependent force density (here, the metric force terms of the cylindrical Navier–Stokes equations). Specifically, we apply our model to the problem of flow past a spherical liquid drop in Re  =  0, Ca regime and, also, flow past a lightly deformed drop. The resulting simulation data, once corrected for the simulation’s inter-facial micro-current (using a method we also advance herein, based on freezing the phase field) show good agreement with theory over a small range of density contrasts. In particular, our data extend verified compliance with the kinematic condition from flat (Burgin et al 2019 Phys. Rev. E 100 043310) to the case of curved fluid–fluid interfaces. More generally, our results indicate a route to eliminate the influence of the inter-facial micro-current.

Keywords: lattice Boltzmann, multiphase flows, hydrodynamics (Some figures may appear in colour only in the online journal)

Introduction
Since 1991, when Gunstensen and Rothman [1,2] invented the technique, a range of increasingly sophisticated multi-component lattice Boltzmann (lB) equation variants have evolved. Still, Gunstensen and Rothman's method remains a milestone of statistical physics [3]. Current multi-component lattice Boltzmann equation (MCLBE) variants depart substantially from Gunstensen's. They may be classified by the physical content of their fluid-fluid interface algorithm. Broadly, where the kinetics of phase separation feature, free-energy methods [4,5] and their thermodynamically consistent extensions, due to Wagner and co-workers [6][7][8], is a natural choice. In continuum, isothermal hydrodynamics, inter-facial fluid physics is defined by dynamic and kinematic conditions [10]. Here, MCLBE simulation may be performed using the phase field or, hereafter, chromo-dynamic method (which is a combination of the algorithms of Lishchuk et al [11] and d'Ortona et al [12]).
To induce fluid-fluid boundary effects, chromo-dynamic MCLBE uses an immersed boundary [13] (henceforth Lishchuk) force which requires appropriate corrections be applied to the velocity [14] and a computationally-efficient, analytic component segregation [12] (which is isotropic, mass-conserving and Galilean invariant). This synthesis is robust and transparent in its encapsulation of Laplacian interfacial tension and its no-traction condition [11]. It also has a very low micro-current which allows direct parametrization of inter-facial tension and width [15]. We note that Reiss and Phillips [16] developed inter-facial perturbation operators in place of the Lishchuk force which, arguably, offer the most consistent encapsulation of inter-facial tension in a kinetic-scale, distribution function-based technique. In its original form, chromo-dynamic MCLBE does not allow for a density contrast between immiscible fluids. Density-difference MCLBE methods of the chromo-dynamic class, based on multi-relaxation time (MRT) collision schemes, due to Ba et al [17], Liu et al [19], Wen et al [20] and Spendlove et al [21] have been developed and benchmarked, in applications such as the evolution of three-dimensional (3D) convective instabilities [17] and, recently, by careful comparisons with theoretical predictions in two-dimension (2D) with restricted symmetry [21,24]. While this previous work certainly confirms the utility of the method, questions relating to its fundamental physical accuracy remain. Here, we aim to assess the extent to which 3D MCLBE simulation complies with the hydrodynamic boundary conditions of mutual impenetrability and the no-traction conditions using, note, comparison with appropriate analytic results for complex flow in the presence of non-uniform curvature. To achieve this, the noise introduced by the inter-facial micro-current and computability present significant difficulties. To overcome the influence of the inter-facial micro-current, we employ a phase-field freezing method, to overcome computational limitations we develop a pseudo-Cartesian methodology, which, itself, requires two extensions to lB to incorporate (i) spatially variable sources/sinks and (ii) velocity-dependent forces. The data we generate to validate boundary condition compliance implicitly supports our practical method for removing the influence of the inter-facial micro-current.
We organise this article as follows. In section 2 we outline chromo-dynamic MCLBE methodology, in section 3 we derive a pseudo Cartesian formulation suitable for our intended applications and develop a suitable MRT scheme, in section 4 we outline its application and in section 5 we present and interpret our data. We conclude in section 6. Details of our analysis of the MRT scheme and high-order accurate lattice stencils are provided in appendices A and B respectively.

Background method
We express the chromo-dynamic variant in colour-blind form, designating immiscible fluid components red and blue, to be described by distribution functions R i (r, t) and B i (r, t) where: An MCLBE for immiscible fluids with a large density difference, based upon an MRT collision scheme for a fluid subject to an immersed boundary, or body force, F α (r) may be formulated as follows: where the density-difference supporting equilibrium, f i (ρ, v), redistributes mass away from the rest (i = 0) link via term φ i : In equation (1), A ij is an element of the collision matrix and source terms F 1i ,..,F 3i correct the dynamics for large density contrasts, external body force densities, D α , and flow sources/sinks respectively. In section 3, we find appropriate forms for F 1i ,..,F 3i . We return to parameter φ i shortly.
An external force density D α is, in general, a sum of distinct contributions. We shall be concerned, here, with velocity-dependent forces which we expose by writing: in which G might represent e.g. gravity, or the Lishchuk force (see equation (8) below). This separation is motivated by the definition of force-adjusted macroscopic observables: Apparently, when a contribution to D depends upon velocity v, equation (3) will define an implicit problem to be solved for v. Above, ρ, ρ R , ρ B , i, δ t , c iα , t i and c s denote overall nodal density, red fluid nodal density, blue fluid nodal density, link-index, time step, the α component of the ith lattice basis vector, the weight for link i and the colour-blind speed of sound (second order tensor lattice isotropy constant). A significant part of our methodology will be to develop explicit solutions for the third of equation (3) for v, when F α (v) depends upon spatial velocity gradients. In implementation, many quantities in equations (20), (21) and (8) will clearly require the numerical computation of gradients, as described below.
Returning to the mass activation parameter now: where α R and α B are chosen to control the fluids' contrast in density: and, for D2Q9, k = 9 5 [24]. In equation (5), ρ 0C is the density deep within the coloured component C = R, B and the third part of the equality reflects mechanical stability of a flat interface. The speed of sound in the red (blue) component is therefore We find κ by noting that α R → α B implies Λ → 1; however, considering equations (2) and (4), we see the traditional equilibrium [18] is recovered for α R = α B = 4 9 : in the latter, the speed of sound (in lattice units, in D2Q9) is 1 , which gives κ = 3 5 and hence c sC = 3(1−αC)

5
. Accordingly, the pressure step across the interface of a red drop suspended in a blue fluid is: Immiscible species are identified in simulation by a generalised colour index, or phase field [17]: in terms of which the Lishchuk force, G (which carries inter-facial tension effects [11]) is: Here, σ is the surface tension parameter and K the interface's mean curvature. ρ N is considered to be continuous: ρ N ∈ [−1, 1] with ρ N = 1(−1) indicative of pure red (blue) fluid and changes occur rapidly in the inter-facial region. Since it identifies the fluid component, the value of index ρ N may be used to control the value of e.g. kinematic viscosity, by introducing into collision matrix eigenvalue λ 3 an appropriate functional dependence upon ρ N : The functional form of λ 3 (ρ N ) has been shown to have significant consequences [22,23]. Kinetic-scale, post-collision colour species segregation, or re-allocation, is an adaptation of the method of d'Ortona et al [12] , which may be written as: where superscript + (++) denotes a post-collision (post re-colour) quantity, β is a chosen segregation parameter [12] and the + (−) sign is used for the red (blue) component. We note that this segregation rule is mass-conserving, simple to implement, local (given a constant director, n) and perhaps most significantly, 'bottom-up', i.e. a kinetic scale postulate which is justified a posteriori. Indeed, it has recently been shown that equation (9) is consistent with the following continuum-scale kinematics (expressed relative to the red component) [24], for uniform fluid motion: In the above, the last term on the right hand side originates in the dynamics correction term, F 1i (see equation (20)) and determines the effect on the model kinematics of the dynamics [24]. When the term 2δ t c 4 s ∂ α ∂ β ρRF αβ ρ in equation (10) is small, Burgin et al [24] find the following, physically correct chromo-dynamics, obtained by solving equation (10): with equivalent behaviour for the blue component. This result represents interface advection within fluids moving at uniform velocity. Note that, with this solution, the corresponding colour field profile is amenable to numerical differentiation i.e. ρ N (r, t) → tanh [βn · (r − vt)], so that the chromo-dynamic field is approximately a material invariant at order δ 2 t . On the other hand, the error term in equation (10) constitutes an issue even with pure advection, i.e. it is present even in uniform flow, which is shown to restrict applicability of the method according to approximate relation Λu = constant [24]. Taking the order δ t approximation in equation (10), DρC Dt ≈ 0 equivalent, Burgin et al [24] arrive at the following, at lowest order of Chapman-Enskog theory: from continuum fluid dynamics. The resulting formulation motivates two, essentially distinct, developments to lB methodology. First, a MCLBE MRT scheme (for a system with density contrast) able to handle spatially distributed flow sources/sinks and second, an advance to the methodology for incorporating velocity-dependent forces in MCLBE. We therefore establish a pseudo 2D representation of 3D, axially symmetric flows, then develop a suitable MCLBE scheme after Dellar et al and finally extend to velocity-dependent forces.
3.1. Two-dimensional representation of axially symmetric flows Figure 1 illustrates the reduction from three to two dimensions, of a multi-component system with axial symmetry. A pseudo 2D system may be derived by discounting ignorable cylindrical polar co-ordinate φ, with an equivalent Cartesian equivalent system (x, z) parametrizing the φ = 0 plane. The Navier-Stokes and continuity equations for incompressible flow, expressed in cylindrical polar co-ordinates, are [10]: 1 r In axial-symmetry, with the following replacements: this system reduces to a 2D, quasi-Cartesian form wherein metric terms of the acceleration equations transform to velocity-dependent accelerations and the continuity equation acquires a source: Above, ∇ = ∂ ∂x , ∂ ∂z and a (A) denotes an acceleration (source/sink): After straightforward algebra, an equivalent velocity-dependent body force density (with units of N m −3 ) and continuity source may be identified for the weakly compressible formulation characteristic of lattice Boltzmann: Since we occupy the φ = 0 plane of a cylindrical polar co-ordinate system with axial symmetry, the Mean curvature in equation (8), for the Lishchuk force, should, however, be computed using the cylindrical polar divergence: Asterisks identify discrete, lattice quantities. Diagrams are not to scale. In (a) the sphere represents a liquid drop of radius R. The 2D system in (b) is described using cylindrical polar co-ordinates (r, z), i.e. φ is ignorable. Cartesian equivalent system (x, z) parameterizes the φ = 0 plane. Our domain is part of the region x ∈ [0.5, ∞), z ∈ (−∞, ∞), i.e. the symmetry boundary r = 0 lies outside simulations. Boundaries are implemented straightforwardly, to achieve the following: (i) an origin of lattice coordinate x * = 0 corresponding to x = 1 2 containing the effect of an axis of symmetry at x = 0, (ii) at large |z * |, periodic replicas of the system and (iii) at large x * mid-link bounce-back, to represent a no-slip boundary located at large x.
Clearly equations (18) and (19) contain spatial gradients. In simulation, such gradients are computed using stencils based upon lattice link vectors, c i . For example, by exploiting lattice tensor isotropy it is possible to devise a compact, third-order stencil : ∂f ∂xα ≈ 3 i t i f (r + δ t c i )c iα , where f denotes the scalar function to be differentiated and δ t is the simulation time step. In appendix B we derive the high order stencils necessary in part of this work.
The continuity equation source/sink, A, introduced above, may be understood in terms of the system's geometry. Consider an annular volume element with volume drdz(rdφ). With v r = constant and ρ = constant, there still exists a mass flux differential between the volume element's curved exterior and interior surfaces, the relevant fluxes being ρv r (r + dr)dφdz and ρv r rdφdz, respectively. The difference in mass flux is equivalent to a source or sink per unit volume of the 2D, pseudo Cartesian domain: Henceforth, we denote by (x * , z * ) the discrete lattice co-ordinate corresponding to (x, z) respectively.

MRT scheme
From section 3.1, we see a quasi-Cartesian lB model for axi-symmetric flow must encompass velocity dependent forces and sources. Here, we derive a MRT scheme using the method of Dellar, able to incorporate geometrical forces, F(v) and sources, A, alongside other simpler body forces, G and also to treat immiscible fluids with a density contrast, Λ.
Dellar developed the most logically consistent approach to an MRT scheme for single comp onent flow [25]. It avoids explicit assignment of collision matrix elements, A ij . The collision matrix is, instead, implicitly defined by its eigenvalues and eigenvectors, a majority being assigned naturally in the Chapman-Enskog process. By a Chapman-Enskog analysis of the kinetic equation (1), using the framework of Guo et al [14], we derive, in appendix A, an MRT scheme-based collision model. Our analysis is performed in D2Q9 with the lattice link vectors c i and indexing defined in figure 2. The scheme is derived colour-blind, so as to clarify the coupling between it and the model kinematics [24]. Put another way, we wish to incorporate the effects of the chromo-dynamics MCLBE segregation or re-colour step in the macroscopic description. Based on an analysis of the f i , we demonstrate correct macroscopic dynamics for the following kinetic sources: Above, λ 3 is an eigenvalue of matrix A, which determines the kinematic viscosity of the lattice fluid, the tensor: Clearly, to use equations (1), (20), (21) and (23) to evolve the f i requires the elements of A. To avoid this expense, we follow Dellar's approach and use a complete set of linearlyindependent macro-scopic modes, m (p ) , p = 0, 1, .., (Q − 1), defined in table 1. A majority of these are physical i.e. m ( p) = i h ( p) i f i represent observables like momentum. However, for closure, a subset of 'ghost modes' require a second lattice link weight function, g i : where, e.g. the designation 'odd' indicates a long lattice link, indexed by an odd value of subscript i (see figure 2) Conveniently, the m (p ) turn out to obey scalar relaxation equations, the parameters for which, we will show, effectively define the elements of A.
We define a projection matrix, comprised of linearly independent left row collision matrix eigenvectors, h ( p) , each the projector of a mode, m (p ) : such that (see table 1). Above, the column vector f ≡ ( f 0 , f 1 , ..., f 8 ) T and M is a column vector, each element of which is a left eigenvector. We define all the h ( p) as polynomial expressions in the lattice basis. Let us project evolution equation (1) using left multiplication by M: where F is the column vector whose elements are the source terms, Using equation (24), the above projected evolution equation decomposes into a set of forced scalar relaxation equations, one for each mode: In equation (25) it is immediate from the properties of the whence we obtain the stated set of scalar, modal relaxation equations. Note that zero eigenvalues are associated with physical modes subject to conservation principles. The problem of developing an MRT scheme now reduces to one of computing appropriate equilibria, m (0)(p ) , and source terms S (p ) . These computations are straightforward when performed within the framework of Guo et al, [14]. In appendix A, we use these authors' approach to write our kinetic equation source (F 1i + F 2i + F 3i ) in terms of the following tensor: Then, in terms of C αβ , above, we have: We note that the source term F i has no projection onto the non-hydrodynamic modes N, J x , J y . For the modal projections of the particle distribution function equilibrium, f i (ρ, v), we have: With sources and equilibria defined, it is possible to write down the full set of modal evolution equations. These are stated for all the Q = 9 modes ρ, (ρu x ), (ρu y ), P xx , P yy , P xy , N, J x and J y in appendix A. An inversion, from mode space, directly to obtain the post-collision distribution function may be performed: Here, a distinct advantage of Dellar's approach is that projection matrix M may be inverted based upon lattice isotropies. Define the components of column vectors k ( p) , p = 0, 1, .., 8 as follows: It is straightforward, using the usual lB simulation lattice properties of isotropy of even tensors (see e.g. equations (A.3)-(A.5)), to show that h ( p) · k ( p ) = δ pp and hence it follows: whence, using equation (28):  (26)), (ii) no explicit collision matrix, A, is constructed and (iii) the post-collision distribution function is constructed directly from post-collision modes, m (p )+ , using equation (28).
Of course, species or colour is finally re-allocated according to f + i , using equation (9).

Velocity-dependent forces
For velocity dependent body force densities, the relationship between the first moment of lB's distribution function and the lattice fluid's macroscopic velocity requires attention. A body force and fluid velocity are generally related according to ρv = i f i c i + 1 2 F. If a component of F depends on v, this relationship may be become implicit. Previous approaches to this general problem have involved approximate 'predictor-corrector' type methods [35] and algebraic solution [36]. Here, we begin by writing: where G denotes any contribution to the total body force density which is velocity-independent, as in equation (34). Here, for F(v) given in equation (18), we find equation (30) may be solved. For the metric force density in equation (18), we obtain from equation (30) two decoupled partial differential equations (PDEs) for v as follows: From equation (31), after a little algebra, we can write a pair of PDEs to be solved for v x and v z : where we have defined: so that S includes any contribution, G, to the body force density which is velocity-independent. Relevant examples of the latter are the Lishchuk force for inter-facial tension effects and a buoyancy force: We seek v x and v z from PDEs (32) using integrating factors. We neglect spatial variation in S in order to obtain a tractable scheme. We will justify this approximation a posteriori. For the first of equation (32), an integrating factor is 1 x exp − x 2 ν and for the second, exp − x 2 ν . Accordingly, we can re-cast the first of equation (32) ν dx and substituting limits, using transformation u = x √ ν and simplifying, we obtain for the pseudo-Cartesian lattice velocity x-component: Numerical integration to determine I(ν, x) will be discussed shortly. Similarly, the second of and selecting φ(y) = 0, we have for the pseudo-Cartesian lattice velocity y -component: The factor I(ν) in equation (35) is evaluated, using Simpson's rule, for a range of x (transverse co-ordinate). Specimen data for ν = 1 6 are shown in figure 3. Note the integrand in I(ν), in equation (35), decays very rapidly indeed. This fact justifies our neglect of the spatial variation in the quantity Sx S0 , above. For a given lattice position (x * , z * ), the measured value of v x is multiplied by I(ν, x * ). Naturally, scale-factors I(ν, x * ) may be pre-compiled for efficient computation. Note that, from equation (36), no adjustment to v z is required in our pseudo Cartesian system.
The discussions above deal with the assignment of an appropriate velocity in a quasi 3D scheme. To simulate, one must of course be able to apply relevant boundary conditions. Appropriate kinetic scale simulation lattice closure rules and representation of the vertical boundaries x = 0 and x → ∞ will be considered in section 4.

Simulations
We note that our pseudo-Cartesian methodology, reported above, reduces the computational expense of the complex flow validation simulations we describe below by at least two orders of magnitude.
To validate our pseudo Cartesian, MCLBE method (and with it the foundations of the chromo-dynamic variants' inter-facial kinematics and dynamics), we shall compare simulation data with analytical solutions. The latter exist, for curved interface configurations, only for steady flow. We first consider steady flow at Re = 0 past a tethered, spherical drop, considered by Rybczynski [31] and Hadamard [32]. (See also [33] for a self-contained and contextualized treatment.) Second, we shall consider the perturbed solution of this flow, solved by Taylor and Acrivos [34], in which the drop deformation in flow is computed as a perturbation expansion in Weber, Wb, and Reynolds, Re, numbers. Figure 1(b) shows the simulation geometry of our test-bench solutions, which both require an unbounded external (blue) fluid.
Throughout this section, discrete lattice positions are indicated by use of an asterisk. The flows represented in schematic figure 1(B) were simulated as follows. Reported data correspond to a simulation lattice with x * = 0,1,..,199 and z * = 0,1,..,299, corresponding to x = 0.5, 1.5, .., 199.5 and z = 0, 1, .., 299. Horizontal periodic boundary conditions were applied between sites z * = 0 and z * = 299, ∀x * . The first lattice nodes in the domain, i.e. x * = 0, were considered to lie at x = 0.5, to avoid a potential singularity in the source, A, at x = 0. Off the lattice boundary x = 0 is a symmetry condition, co-located with an inviscid, solid boundary condition, v n = 0. Both these physical conditions were encapsulated in the applied lattice closure rule of mid-link specular reflection. A constant, positive body force density −G 0 ρ Rêz was applied. This body force produced acceleration, even for a neutrally buoyant drop, note.
For all data presented and discussed, the non-hydrodynamic, ghost, modes of our Dellartype scheme were relaxed directly to equilibrium, by setting λ i = 1 (i = 3) and the fluids' kinematic viscosity was ν = 1 3 = constant, (corresponding to the MRT relaxation parameter λ 3 = 2 3 ). To facilitate comparison, our data will not depict the external flow. Of course, all velocity field data reported actually correspond to 3D, axially-symmetric flow, confined within (say) meridional planar section φ = 0, taken from the equivalent 3D system.

Flow past a tethered spherical liquid drop
For the simulations performed on non-deforming drops, we shall utilize the fact that our testbench analytical solution gives the flow field at all points. Accordingly, we present data for the entire flow field, which will facilitate our examination of inter-facial boundary conditions.
In the Re,Ca =0 regime, the anticipated micro-current activity will be relatively large, since micro-current amplitude is proportional to surface tension σ and the physical flow (value of U 0 ) must be small, to avoid deformation. Nevertheless, this artefact can be subtracted, to reveal the physical flow [35], even when it seemingly dominates. For clarity, we shall present data on internal flow within spherical drop and assume the kinematic viscosity of the system is uniform. Hence, the ratio of the shear viscosities is determined by the density contrast: The flow is most conveniently computed using a spherical polar co-ordinate system [33]. The analytical calculations [31,32] represent a stringent test-bench and it makes explicit reference to key dynamic and kinematic boundary conditions, the representation of which in MCLBE simulation is our present concern. Note that the form of the theory quoted below does not contain explicit reference to inter-facial tension i.e. it refers only to kinematic and viscous, no-traction conditions at the interface and assumes surface tension forces are so strong as to enforce spherical shape. These constraints are raised when we come to consider a lightly deformed drop. Expressed in spherical co-ordinate system with the origin at drop centre, the internal velocity field is: v s r = cos(θ) Here, U 0êz is the constant speed of the blue fluid at a large distance from the red drop and r s = √ r 2 + z 2 is spherical polar distance from the origin. Note the stagnation point -the centre of a single internal vortex-is always located at θ = π 2 , r s = R √ 2 or, equivalently, x = R √ 2 , z = 0. This flow feature does not change position when Λ changes. The above solution is plotted as a vector field for U = 1.0, R = 20, in figure 4 (panel (A)).
In simulation, the steady-state of flow was identified by a constant velocity residual and the surface tension parameter, σ, used was always as small as possible for a front-back symmetric drop. A Galilean transformation to the rest frame of the drop was applied to the z velocity component (only), once steady-state was reached. This was achieved by subtracting the average velocity of the red fluid, the latter being computed as a sum of annular mass increments: Hence, the simulated flow field was translated to the approximate rest frame of drop, translated z * co-ordinates being rounded. Front-back symmetry and spherical drop shape were used to confirm that data correspond to Re ≈ 0. (Formally, the computed Re = <vz>R ν < 3.0 × 10 −4 .) The phase field of the steady state was then frozen, the applied body force removed and the corresponding micro-current flow allowed to evolve to its steady state. Finally, this flow was measured and subtracted. For all data reported in the next section, the value of phase field segregation parameter was β = 0.67 and the buoyancy parameter was G 0 = 8.0 × 10 −9 lattice units. Data reported correspond to R * = 20 lattice units and were checked for convergence to the narrow interface limit and for general finite size effects as follows. Our check for finite size effects was performed by taking a system with size parameters R * = 20 lattice units, x * ∈ [0, 199], z * ∈ [0, 299] and repeatedly doubling x * , z * whilst keeping R * fixed; our check for the narrow interface (continuum) limit was performed by repeatedly doubling all of x * , z * , R * .

Flow past a tethered deforming drop
For these simulations, we shall present data for the drop shape deformation alone, for which analytic expressions exist. This validation is therefore more challenging to interpret that that reported in section 4.1 above.
Taylor and Acrivos use a singular perturbation solution of the antisymmetric equations of motion to predict the shape (and drag) of a slightly deforming drop [34]. We shall be concerned with Taylor and Acrivos' equation (30) and will use their notation. With our restriction of equal kinematic viscosities, the deformed drop radius at spherical polar zenithal location, θ, may be written: with expansion co-efficients [34]: where P 2 (x) etc is the second order Legendre polynomial [29] and: No confusion should arise from our use of symbol λ above, as the MRT collision eigenvalues, λ p , p = 1, ..., 3. In simulation, we compute dimensionless Reynolds and Weber numbers as follows: where ∆p is the measured pressure step across the drop interface, from equation (6). Simulations were performed as described in section 4.1, above. Data for the deformed radius were extracted and fitted to equation (40) as set-out below, by first computing the centre of gravity as: and from it the radial distance and zenithal angle of inter-facial locations were computed as follows: R(θ * ) = (x * + 0.5) 2 + (z * − < z >) 2 1/2 , A conjugate gradients grid search optimised a fit to equation (40), using as adjustable parameters the undeformed drop radius, R, and amplitudes a 2 , a 3 should, apparently, produce results related to Wb and Re after equation (41). For deformed drop data presented in section 5, the value of the phase field segregation parameter was reduced to β = 0.3, with a commensurate doubling of undeformed drop radius (to R * = 40) and lattice dimensions. This was necessary to stabilize simulations, whilst maintaining the width of the inter-facial region (which varies as 1/β), relative to R, constant compared with data in section 4.1. The surface tension parameter σ = 6.0 × 10 −4 lattice units, Wb 2.0 × 10 −2 and Re < 0.1.

Results and discussion
All data presented and discussed in this section correspond to maximum density contrasts Λ < 10. We return to this matter in our conclusions. The main interest of this article is not the stability of the chromo-dynamic density difference method, rather it is the fundamental accuracy of its hydrodynamics. In this respect, results are encouraging. Consider our first test-bench solution of flow past a spherical, red tethered drop, as discussed in section 4. Figure 4 below shows, side by side, the full velocity field of the shifted analytical solution, computed from spherical polar co-ordinate velocity components in equation (37) (panel (A)) and the corresponding velocity field obtained from simulation (panel (B)). The latter were adjusted for the micro-current as discussed shortly. For these data, Λ = 1. The correspondence between these two solutions over the whole internal domain is excellent.
Panel (A) of figure 5 shows steady-state simulation data for the scaled internal velocity field, after applying the Galilean transformation in equation (39). In these data, the effect of the inter-facial micro-current is very clear and flow is certainly not parallel to the interface (the solid black contour corresponds to the interface center, ρ N = 0). A localized micro-current circulation may be subtracted from the overall flow in (A), to reveal the physical flow in the linear, Stokes' regime [35]. The flow in panel (B) of figure 5, which is not normalised to the same quantity as that in panel (A), corresponds to the micro-current of the 'frozen' colour field, ρ N . This is subtracted, vectorally, to expose the physical flow field in panel (C). The correspondence with theory and, in particular, the compliance of the flow with the kinematic condition of mutual impenetrability are striking in these data. Whilst a MCLBE interface is not sharp, flow is directed tangentially to contours of ρ N = constant everywhere within a shell of finite thickness and we can say that the physical flow lies tangent to the curved interface at all relevant points in the simulated domain. The significance of this result should be stressed. Apparently, in this challenging Re, Ca =0 regime, where received wisdom would argue that the influence of the inter-facial micro-current is proportionally large in comparison with the  (That is, the vortex is predicted to lie at x = r = R √ 2 , ∀Λ ). Moreover, the extent of compliance with the kinematic condition of mutual impenetrability across the range of Λ is also clear in these data. true, hydrodynamic signal, our data, together with the supporting analysis contained in [24], demonstrate that, at steady state, its influence may be removed altogether, simply by computing, then subtracting the flow induced by the frozen phase field, ρ N . Figure 6 assesses the quantitative variation of the velocity field components from figure 5(C), measured along two transects from the lattice co-ordinate closest to the drop centre. To facilitate comparison with the theory of equation (37) we have defined: Panel (A) compares the spherical polar radial velocity, v s r , with its measured value along the equator (θ = 0), where there is no tangential component of motion, note. Panel (B) compares spherical polar v s r with its measured value along the line subtending an angle of π 4 at the positive z axis, intersecting the latter close to the centre of the drop (i.e. θ ≈ π 4 ). Panel (C) : as panel (B) but for spherical polar v s θ . In these data, correspondence between theory and simulation is weakest for r → R as expected, owing to the diffuse nature of the MCLBE interface. In simulation, the interface corresponds to a shell of finite thickness, not a surface. Again, we observe that despite its diffuse nature, both the magnitude and direction of flow conform with a sensible interpretation of the relevant continuum no-traction and kinematic conditions. In figure 7, we show the simulated internal flow field (with the inter-facial micro-current removed, as discussed above), for flow past an effectively tethered red drop computed, for a range of density contrasts, Λ, all in Stokes' regime. (All flows have Re< 2 × 10 −4 .) Note that the buoyancy parameter varies between the data on panels (A)..(D). Insofar as one can judge from these data, the location of the primary vortex does not change, which agrees with the prediction implicit in equation (37). Moreover, the extent of compliance with the kinematic condition of mutual impenetrability, across the range of Λ is again very clear in these data. We defer all further discussion until section 6. We proceed to consider our second test-bench solution of flow past a slightly deformed drop. It will be noted that the range of Λ achieved with the chosen level of resolution, for this more complex interface shape is reduced relative to that in e.g. figure 4. We return to this issue in section 6. Figure 8 shows data obtained for a drop sedimenting vertically downwards, with deformation, viewed from the rest frame of the drop, at steady state. For these data, Λ = 5. Panel (A) shows the pressure field (in lattice units), panel (B) shows Stokes' stream-function for the flow, whilst the red contour corresponds to the ρ N = 0 contour, which is the nominal interface. In figure 9 we see the time development of deformation in a set of interface configurations for a drop sedimenting vertically downwards. Here, the red line corresponds to the initial ρ N = 0 contour, the black contour is a later time and the blue contour corresponds to steady state.
In figure 10 we see points on the lightly deformed drop interface for Wb ≈ 2.0 × 10 −2 , Re ≈ 0.01, for three lightly deformed drops with Λ = 1 2 , 1, 2 (panels (A), (B), (C) respectively). The solid line corresponds to a conjugate gradients grid search optimised fit to the prediction of equations (40)-(42). The measured fit coefficients corresponding to these data are recorded in table 2, below.
The fit to the interface is in reasonable agreement with the predictions of analytical perturbation theory.

Conclusion
Here, we have developed a 2D pseudo Cartesian, chromo-dynamic MCLBE scheme for efficient simulation of multi-component systems with 3D axial symmetry, continuously distributed sources (and/or sinks) of momentum and a density contrast between the simulated fluids. This methodology makes the complex inter-facial validations we report computationally accessible. Without our 2D pseudo Cartesian model, our simulations would have been at least two orders of magnitude more expensive. Our analysis motivates a novel and versatile treatment for velocity-dependent forces in lattice Boltzmann methodology (here, the metric force terms of the cylindrical Navier-Stokes equations), which should facilitate wider application of the method, where such forces are present e.g. magneto-hydrodynamics, porous flow and geophysical flow. All the data we present with our novel scheme correspond to maximum density contrasts of Λ ≈ 10. The investigations of Burgin et al [24], using a similar chromo-dynamic MRT model in 2D only achieve larger Λ in the case of no applied flow. Moreover, the complex simulations with chromo-dynamic MCLBE method [20] have Λ = 2. Apparently, the density contrasts achievable with a given level of resolution appear to be less than those achieved with the MCLBE method of e.g. Innamuro et al [37]. Nevertheless flows of immiscible fluids at The methodological developments we relate here facilitate our principal aim: validation of interfacial hydrodynamic properties of chromo-dynamic MCLBE, with complex interface shapes, by computation of solved problems, namely flow past a spherical liquid drop for Re,Ca → 0 and, also, a lightly deformed drop. A glance at figure 5 illustrates the question. In panel (A), the dominant micro-current produces a strong transverse flow, which clearly violates the kinematic of condition of mutual impenetrability. However, in the case of a spherical drop, simulation data, once corrected for the inter-facial micro-current (see below), actually show excellent agreement with well-accepted analytical theory, over a range of density contrasts and strongly imply the operation of the kinematic condition of mutual impenetrability (e.g. figure 4) and the no-traction condition (e.g. figure 6). Broadly, agreement between data from chromo-dynamic MCLBE variants (including those with sources, herein) and the kinematic condition of mutual impenetrability is very good for flat interfaces [24] and, we now confirm, also for curved interfaces in 3D flows. This result clearly adds considerable support to a number of recent chromo-dynamic MCLBE, MRT schemes [17,19,20], as well as future applications with the method. Further, for lightly deformed drops, simulation data are in reasonable agreement with analytical theory, The excellent agreement obtained in the challenging Re =0, Ca → 0 regime, in which the inter-facial micro-current is a relatively large contribution to the total flow, shows the precaution of phase field freezing, outlined in section 5, to be a very effective means of resolving the physical flow alone. Generally, the success of this procedure supports the view that the interfacial micro-current is a super-posed, hydrodynamic response in the linear regime (which implies that interfacial micro-currents in other MCLBE variants can be removed in a similar fashion). It also provides a way forward for the use of MCLBE throughout suspension rheology, especially in the limit of a high volume fraction of the suspended phase, where microcurrent activity infects -indeed, dominates-hydrodynamic signals over most of the domain.

Appendix A. Multi-relaxation-time scheme for large density contrast immiscible fluids with density sources/sinks
We present a detailed derivation of the Navier-Stokes equations from the multiple-relaxationtime (MRT) lattice Boltzmann equation adapted for multi-component applications with a large density difference between completely immiscible components, where a body force is present. The latter is necessary to carry the Lishchuk interface force. Density gradients associated with a chromo-dynamic or phase field must not affect the macroscopic dynamics, of course. Equally the dynamics of the developed scheme should not affect the physics of the segregation rule. The challenge is to compensate for their presence accurately, whilst retaining algorithmic stability and simplicity. The key advance outlined in this appendix is the consideration of fluid cources and sinks.
In the interest of a compact literature, we retain here the overall structure of the analyses of Guo et al [14], Dellar [25] and Hou et al [26]. For this reason we work in this appendix in two dimensions, using space variables x and y and we denote flow velocity u. We choose to extend the scheme of Dellar because it is efficient (due to a careful choice of non-hydrodynamic modes N, J x and J y [25] with zero equilibria), robust, straightforward to implement and, not least, logical. Note, within this and following appendices, i = At the kinetic scale, the 'forced' MRT LBE for a system subject to an 'external' force term can be expressed in the following form: where the density-difference supporting equilibrium distributes mass away from the rest ( j = 0) link via term φ j , is in the form of: and where the kinetic equation source term, F i , is assumed to have the following properties: where n is a scalar to be determined and tensor C is also to be determined.
In this appendix, we first set-out the basics, then proceed to the Chapman-Enskog analysis to obtain the thermodynamic limit of the kinetic scheme defined in equation (A.1) (i.e. find appropriate expressions for tensor C, which represents the crux of the problem of recovering correct hydrodynamics with the MRT scheme), then we transform to a modal description, and finally, we invert that transformation to obtain an explicit expression for the post-collision distribution function. To maintain parity with the analysis of Guo et al [14] at the outset, we now 'relax' the form of the definition of lattice velocity given in equation (3) as follows: with m being a constant to be determined. Dellar's [25] eigenvalues and left row eigenvectors for the collision matrix A ij can be tabulated as in table 1, where we define for α, β = x, y, and the Π ( p) αβ , p = 0, 1 have the usual meaning, which is also set-out later in this appendix.
As set out in table 1, matrix A ij has the following properties which, it will be seen, are necessary if one is to recover correct hydrodynamics: Here α and β represent either x or y direction in the lattice grid. We also assume that the lattice basis c i and the corresponding weights t i have the following symmetry properties: Let us consider the most rapid behaviour in the model. Equation (9b) of [14] (Guo et al) can be obtained straightforwardly as the following: Taking summation i on both sides of equation (A.8) leads to: Using property (A.3), we therefore obtain: which gives the continuity equation and is the MRT equivalence of equation (10a) in [14]. Before proceeding further, we present the counterpart result, deriving from the kinematics of our model. It is reasonable to interpret equation (11), which is order δ t accurate, as corresponding to the statement DρR Dt = DρB Dt = 0, from which it is straightforward to deduce [24]: That is, on the shortest timescales, the chromo-dynamic field is a material invariant. We will indicate where we appeal to this fact, to eliminate certain t 0 derivatives, as we proceed.
Multiplying every term of equation (A.8) by c ix and taking summation i on both sides, we have: where i F i c ix = nF x , n is a constant to be determined, and we will use property (A.4). Similar as in [26] ( Hou et al), the momentum flux tensor is defined as: where the zeroth-order momentum flux tensor Π (0) αx = (2φ 1 + 4φ 2 )ρδ αx + ρu α u x and δ αx is the Kronecker delta. Equation (A.12) is the MRT equivalence of (10b) in [14] Guo et al. We note that the equivalent result in Guo et al [14] couples n, m and τ (the collision parameter) in the case of the single-relaxation time (SRT) variant. Following Guo et al [14] we recover the appropriate form of the Euler equation by setting: with no constraint on m at O( ).
Proceeding to O( ) 2 in the Chapman-Enskog expansion, equation (9c) of [14] can be rewritten as the following: The source term F i in equation (A.1) may now be conveniently partitioned into: (i) a term responsible for correcting the model dynamics in the presence of density gradients associated with component changes, (ii) a term responsible for the Lishchuk force and (iii) a term resulting from the continuity equation not equalling to zero: where: From equations (20) and (A.22) we identify: The development of an MRT scheme for MCLBE after Dellar et al now reduces to one of recasting the model in terms of a set of linearly-dependent modes (defined in table 1). This process is described in section 3.2. The particular modal equations, determined from equation (1) and table 1 are the following 'forced' modal evolution equations of simple, scalar relaxation: where α ∈ [x, y]. We note the simple form of the relaxation equations for h (6) , · · · , h (8) , i.e. N, J x , J y , which for λ 6 = λ 7 = 1, reduce to N + = J + x = J + y = 0. This is a direct consequence of the choice of equilibria. Note that the fact that the scheme devised contains source and sinks is apparent from the first of the above modal equations.