extfieldsofstudy:
Computer Science
Mathematics

text:
Primal and mixed finite element formulations for the relaxed micromorphic model

The classical Cauchy continuum theory is suitable to model highly homogeneous materials. However, many materials, such as porous media or metamaterials, exhibit a pronounced microstructure. As a result, the classical continuum theory cannot capture their mechanical behaviour without fully resolving the underlying microstructure. In terms of finite element computations, this can be done by modelling the entire body, including every interior cell. The relaxed micromorphic continuum offers an alternative method by instead enriching the kinematics of the mathematical model. The theory introduces a microdistortion field, encompassing nine extra degrees of freedom for each material point. The corresponding elastic energy functional contains the gradient of the displacement field, the microdistortion field and its Curl (the micro-dislocation). Therefore, the natural spaces of the fields are $[\mathit{H}^1]^3$ for the displacement and $[\mathit{H}(\mathrm{curl})]^3$ for the microdistortion, leading to unusual finite element formulations. In this work we describe the construction of appropriate finite elements using N\'ed\'elec and Raviart-Thomas subspaces, encompassing solutions to the orientation problem and the discrete consistent coupling condition. Further, we explore the numerical behaviour of the relaxed micromorphic model for both a primal and a mixed formulation. The focus of our benchmarks lies in the influence of the characteristic length $L_\mathrm{c}$ and the correlation to the classical Cauchy continuum theory.

Introduction
A common problem in the computation of materials with a pronounced micro-structure is the internal complexity of the geometry. In order to fully capture the kinematics one might resolve the underlying micro-structure. This can be done either with multi-scale finite element methods [21] or by modelling the finite element mesh to fully incorporate the microstructure. In both cases, the computational cost increases, leading to longer computation times and decreasing applicability of such approaches. Alternatively, one may enrich the mathematical model in order to account for the increase in the kinematical complexity. This approach gives rise to generalized continuum theories such as higher order gradient methods [29,36,48] or micromorphic continua [28,44,67]. Micromorphic continuum theories extend the kinematics of the material point with additional degrees of freedom, the choice of which defines the specific theory. Common examples are micropolar Cosserat [27,39], microstretch [60], and microstrain models [23,26]. The latter represent sub-types of the micromorphic continuum. In its most general case, micromorphic continua assume an affine deformable microbody for each material point. As such, this deformation, called here the microdistortion P , is fully captured by three-by-three matrices and introduces nine extra degrees of freedom. The full micromorphic continuum, as introduced by Eringen and Mindlin [22,37], incorporates the gradient of the microdistortion DP into the free energy functional. The resulting hyperstress term is a third order tensor. As such, it is unclear how this term is to be interpreted or applied. Typically, micromorphic continuum theories introduce a characteristic length scale parameter L c , which abstractly relates the dimension of the micro-body to that of the macro-body. In the full micromorphic model, when the characteristic length L c becomes very large, the microdistortion P must become constant in order for the theory to generate finite energies, which may lead to boundary layer problems [66].
The relaxed micromorphic continuum theory [46] takes a different approach by instead incorporating only the Curl of the microdistortion Curl P into the free energy function. The latter term, known as the micro-dislocation, remains a second order tensor and as such, induces a matrix-valued right-hand-side, known as the micro-moment M . Further, large characteristic lengths L c maintain finite energies [47,66]. The theory aims to capture the mechanical behaviour of both highly homogeneous materials and materials with a pronounced microstructure by governing the relation to the classical Cauchy continuum using the characteristic length L c [6] and shows great promise with respect to applications utilizing metamaterials, such as band-gap materials [7,16,34,35] and shielding against elastic waves [54,59]. Furthermore, analytical solutions have been derived for bending [56], torsion [55], shear [57], and extension [58] kinematics. The inclusion of the Curl of microdistortion in the free energy function implies the existence of unique solutions [24,45] in the space X = [H 1 ] 3 ×[H (curl)] 3 , as shown in Section 3.1.1. While the construction of finite elements for the space H 1 is well-known in the field of mechanics, H (curl)-finite elements are commonly used for the Maxwell equations, for example in magnetostatics [64]. For the construction of finite elements for H (curl) one may use Nédélec subspaces [40,41]. The formulation of higher order elements is detailed in [68]. For large characteristic length values L c the computation with the primal formulation (X h ⊂ X = [H 1 ] 3 × [H (curl)] 3 ) becomes unstable [66]. However, it can be re-stabilised by using a mixed formulation. The latter requires the employment of Raviart-Thomas- [53] or Brezzi-Douglas-Marini elements [13] and fully discontinuous finite elements as per the de Rham diagram [5,19,20].
In this work we demonstrate the existence and uniqueness of the primal formulation using the Lax-Milgram theorem, Section 3.1.1. Further, we introduce a mixed formulation, which is stable for large characteristic length L c values, Section 3.1.3. In Section 3.2 we derive the corresponding convergence rates for the discrete spaces. The construction of lower order finite elements for both formulations is explained in Section 4, with focus on a solution to the orientation problem and application of the discrete consistent coupling condition. Section 5 is devoted to numerical benchmarks of the finite element formulations but also features the convergence characteristics of higher order elements using NETGEN/NGSolve [61,63]. Finally, we present our conclusions and outlook in Section 6.

The linear relaxed micromorphic continuum
The linear relaxed micromorphic continuum [43,45,46] is described by its free energy functional, incorporating the gradient of the displacement field, the microdistortion and its Curl I(u, P ) = 1 2 Ω C e sym(Du − P ), sym(Du − P ) + C micro sym P , sym P (2.1) + C c skew(Du − P ), skew(Du − P ) + µ macro L 2 c CurlP 2 dX − Ω u, f + P , M dX , with u : Ω ⊂ R 3 → R 3 and P : Ω ⊂ R 3 → R 3×3 representing the displacement and the non-symmetric microdistortion, respectively. Here, C e and C micro are standard fourth order elasticity tensors and C c is a positive semi-definite coupling tensor for (infinitesimal) rotations. The macroscopic shear modulus is denoted by µ macro and the parameter L c ≥ 0 represents the characteristic length scale motivated by the geometry of the microstructure. The body forces and micro-moments are denoted with f and M , respectively. The differential operators are defined as Du =   u 1,1 u 1,2 u 1,3 u 2,1 u 2,2 u 2,3 u 3,1 u 3,2 u 3,3   , CurlP =    (curl P 11 P 12 P 13 ) T (curl P 21 P 22 P 23 ) T (curl P 31 P 32 For isotropic materials the material tensors have the following structure C e = 2µ e S +λ e 1 ⊗ 1 , C micro = 2µ micro S +λ micro 1 ⊗ 1 , where S : R 3×3 → Sym(3) and A : R 3×3 → so(3) are the fourth order symmetry and anti-symmetry tensors, respectively. Taking variations with respect to the displacement u δ u I = Ω C e sym Dδu, sym(Du − P ) + C c skew Dδu, skew(Du − P ) − δu, f dX , (2.4) and the microdistortion P δ P I = Ω − C e sym δ P , sym(Du − P ) − C c skew δ P , skew(Du − P ) + C micro sym δ P , sym P + µ macro L 2 c Curl δ P , Curl P − δ P , M dX , (2.5) Figure 1: Depiction of the reference domain Ω with an underlying micro-structure, the Dirichlet boundary Γ D according to the consistent coupling condition and the Neumann boundaries Γ u N and Γ P N for the displacement and the microdistortion, respectively. The vector ν is the unit normal vector on the boundary. The unit tangent vectors are illustrated with τ . The body force is given by f and the micro-moment by M . yields the symmetric bilinear form and the linear form for the load The corresponding strong form follows from partial integration, see Appendix A where ν denotes the outer unit normal vector, such that P × ν is the projection to the tangent surface on the boundary. The terms u and P are the prescribed displacement and microdistortion fields on Γ u D and Γ P D , respectively, see Fig. 1.
From a physical point of view, it is impossible to control the micro-movements of the material point on the Dirichlet boundary without also controlling the displacement. Consequently, the relaxed micromorphic theory introduces the so called consistent coupling condition [17] where the prescribed displacement u on the boundary automatically prescribes the tangential component of the microdistortion P on the same boundary, effectively inducing the definition Γ D = Γ P D = Γ u D .
3 Solvability and limit problems

Continuous case
In this section existence and uniqueness of the weak formulation of the relaxed micromorphic continuum for both a primal and a mixed method is discussed. In order to simplify the proof, homogeneous boundary conditions are assumed on the entire boundary, Γ u D = ∂Ω. The proof can be easily adjusted for mixed or inhomogeneous boundary conditions as long as the Dirichlet boundary does not vanish |Γ u D | > 0.
For the following formulations, we define corresponding Hilbert spaces and their particular norms from which we derive the corresponding spaces for our higher dimensional problem where both spaces are to be understood as row-wise matrices of the vectorial spaces.
Remark 3.1. In the following sections we assume a contractible domain Ω for all proofs.

Primal form
The proof of existence and uniqueness has already been given in [46] and subsequently generalized to the dynamic setting in [24]. For the reader's convenience and for later use we present a proof based on the Lax-Milgram theorem and Korn's inequalities for incompatible fields [8, 9, 31-33, 49, 50], similar to [45], where only the case µ c = 0 together with Γ P D = ∂Ω has been considered. Therefore, we define the product space X = [H 1 0 (Ω)] 3 × H (Curl, Ω) with its standard product norm {u, P } 2 X = u 2 H 1 + P 2 H (Curl) and examine the problem: Find {u, P } ∈ X such that where a(·, ·) and l(·) are defined as in (2.6) and (2.7), respectively.
, and C e , C micro be positive definite on Sym(3). Further, assume that C c is positive definite on so(3), or positive semi-definite and Γ P D = ∅. Then Problem 3.3 is uniquely solvable and there holds the stability estimate Proof. We show continuity and coercivity of (2.6). During the proof we denote with c > 0 a generic constant which may change from line to line. Continuity follows immediately with Cauchy-Schwarz inequality For the coercivity we first consider the case of positive definiteness of the tensors First, we consider the symmetric terms. Analogously to [66], with Young's inequality 1 and Korn's inequality 2 [42] we obtain We can choose ε = 1 2 (1 + c1 c1+cmicro ) such that both terms are positive. Next, we estimate the skew-symmetric part With 0 < δ < 1 only the second term is positive. By combining both estimates we conclude by choosing Thus, with Poincarè-Friedrich's inequality 3 we obtain the coercivity In the case of a positive semi-definite C c (even C c = 0 is allowed: absence of rotational coupling) together with Γ P D = ∅ we must use the generalized Korn's inequality for incompatible fields, cf. [9,9,24,32,33,49,50] sym P 2 and estimate Continuity of the right-hand side l is obvious and thus we can apply the Lax-Milgram theorem finishing the proof.
Remark 3.2. The proof above fails if we consider P ∈ [H 1 (Ω)] 3×3 instead of H (Curl, Ω). In fact, a(·, ·) is then no longer coercive as we cannot bound P 2 L 2 + Curl P 2 L 2 ≥ c P 2 H 1 uniformly in P . In [66] and [65] it is demonstrated that already in two spatial dimensions a loss of optimal convergence rates is obtained when using P ∈ [H 1 (Ω)] 2×2 . Remark 3.3. Note that in [51,52] it is shown that for certain sufficiently smooth data, the regularity of P can yet be improved to P ∈ [H 1 (Ω)] 3×3 .
3.1.2 Limit of vanishing characteristic length L c → 0 In the limit L c → 0 the bilinear form (2.6) reduces to (3.10) Therefore, we lose control over the Curl of P yielding a loss of regularity for P from H (Curl, Ω) to [L 2 (Ω)] 3×3 . We emphasise that the proof of Theorem 3.1 can be directly applied with the adapted product space together with the requirement of positive definite C c on the set of skew-symmetric matrices. Otherwise, control over the skew-symmetric part of P is completely lost in the limit leading to unstable results. Note that in this case no prescription of boundary conditions for P is possible.
Eq. (2.8b) can be reformulated as and used to express P algebraically P = (C e + C micro ) −1 C e sym Du + skew Du +M . Applying the latter to Eq. (2.8a) yields − Div(C macro sym Du) = f , The upper definition is derived in [6] and relates the meso-and micro-elasticity tensors to the classical macroelasticity tensor C macro of the Cauchy continuum, allowing to extract the macro material constants in the isotropic case .
In fact, C macro contains the material constants that arise from classical homogenization for large periodic structures. From Eq. (3.12) we observe that in general P is not a gradient field, even if M = 0. Therefore, by setting λ e = λ micro = 0 we obtain that skew P = skew Du and sym P = µe µe+µmicro sym Du. As µ micro > 0 we deduce that P is not a gradient field. This is a significant deviation from the two dimensional relaxed micromorphic model of antiplane shear analyzed in [66], where a gradient field as a right-hand side leads to a gradient field for the microdistortion.

Mixed form
A major aspect of the relaxed micromorphic continuum is the relation to the classical Cauchy continuum theory. This relation is governed by the material constants, where the characteristic length L c plays a significant role.
We are therefore interested in robust computations with respect to L c . To that end, we reformulate the problem as a mixed formulation. The first step consists in introducing the new unknown D = µ macro L 2 c Curl P ∈ H (Div, Ω) , (3.16) reminiscent of the micro-dislocation, and examining its distribution with a test function In fact, D must be solenoidal as and thus the appropriate space for D is We again assume that Dirichlet boundary conditions are prescribed on the whole boundary and thus, from P ∈ H 0 (Curl, Ω) there follows D ∈ H 0 (Div, Ω) as Dν = Curl P ν = 0.
The (bi-)linear forms are now given by where the Lagrange multiplier D has the physical meaning of a hyperstress. Notice we now approximate the hyperstress directly with its own variable, therefore recovering lost precision due to differentiation.
Therefore, we have to prove that a(·, ·) is coercive on the kernel of b(·, ·) The last equality follows from the exact property Curl H 0 (Curl, Ω) = • H 0 (Div, Ω). However, we already know that a({u, P }, {δu, δP }) + Curl P , Curl δP L 2 is coercive, compare the proof of Theorem 3.1. This leaves us with the Ladyzhenskaya-Babuška-Brezzi (LBB) condition to be satisfied We choose u = 0 and P such that The existence of such a potential P for D ∈ • H 0 (Div, Ω) follows directly from the theory of stable Helmholtz decompositions (3.23).
Thus, with Brezzi's theorem [11], there exists a unique solution independent of L c for 1 µmacro L 2 c ≤ 1 fulfilling the stability estimate where the constant c does not depend on L c .

Remark 3.4.
In the proof we made use of the exactness property of the de Rham complex by restricting the space of D to be divergence free. Otherwise we cannot prove the LBB condition in the limit L c → ∞. Numerically we observed that this restriction is necessary, otherwise the resulting matrix becomes singular.
The limit case lim L c → ∞ of Eq. (3.21) is well-defined, resulting in the problem: Consequently, at the limit lim L c → ∞ we have Curl P = 0, i.e. the rows of P are gradient fields. This can be observed by considering that δD only tests functions in range(Curl), where the range is fully given by range(Curl) = • H (Div). As such, by the Helmholtz decomposition only the gradient part of the microdistortion P remains. The latter is also clearly observable through the minimization of the energy function, as P = DΨ is needed for finite energies. Therefore, the consistent coupling condition is crucial, prescribing only gradient fields as Dirichlet boundary conditions for P . From the theory of mixed methods we obtain, analogously to the 2D case in [66], quadratic convergence in L c towards the limit case Remark 3.5. Note that for the limit L c → ∞ one finds P → DΨ, inducing finite energies due to Curl DΨ = 0.
Since L c defines a zoom into the microstructure, the latter can be interpreted as the entire domain being the micro-body. Consequently, setting

30)
and as such, taking the divergence of Eq. (2.8b) results in The divergence of the micro-moment Div M can be interpreted as the micro body-force. The latter implies that the limit L c → ∞ defines again a classical Cauchy continuum theory with a finite stiffness C micro , representing the upper limit of the stiffness for the relaxed micromorphic continuum. Further, due to the consistent coupling condition and Eq. (3.30) there holds DΨ = Du (for a thorough derivation see [6]).

Reformulation of the mixed divergence free constraint
The construction of finite elements which are exactly divergence free, H (Div, Ω), is possible for higher polynomials, [68]. However, at least for the lowest order shape functions, one has to use an additional Lagrange multiplier q ∈ [L 2 0 (Ω)] 3 , compare Remark 3.6 below, forcing Div D = 0. We now show how the mixed formulation has to be adapted leading to a method which can be directly implemented. By defining X = X ×[L 2 0 (Ω)] 3 we introduce the adapted bilinear forms and the problem reads: find ({u, P , q}, D) ∈ X × H 0 (Div, Ω) such that Theorem 3.3. Problem 3.33 is uniquely solvable and the solution solves also Problem 3.21.
Proof. The kernel of the bilinear form b is now given by The last equality follows by choosing on the one hand D = Curl P ∈ H 0 (Div, Ω) yielding Curl P = 0 and on the other hand from the de Rham sequence we find that Div (H 0 (Div, Ω)) = L 2 0 (Ω) such that we can choose D ∈ H 0 (Div, Ω) such that Div D = q. With the same argument as in the proof of Theorem 3.2 the kernel coercivity follows. For the LBB condition we use the Helmholtz decomposition

Now, with
Brezzi's theorem we conclude that Problem 3.33 is uniquely solvable with a stability constant independent of L c . Due to the exactness property of the de Rham complex the additional Lagrange multiplier q enforces that D ∈ • H 0 (Div, Ω). Therefore, taking ({u, P }, D) ∈ X × • H 0 (Div, Ω) from the solution solves also Problem 3.21.
Remark 3.6. In case of prescribed Dirichlet boundary conditions on Γ P D for P ∈ H (Curl, Ω), due to the de Rham complex, compatible Dirichlet conditions have to be used for D ∈ H (Div, Ω). If Γ P D = ∂Ω then from P ∈ H 0 (Curl, Ω) there follows D ∈ H 0 (Div, Ω). Thus, with Div(H 0 (Div, Ω)) = [L 2 0 (Ω)] 3 the space for q has a zero mean. Consequently, another vector valued variable in R 3 is needed for numerics to ensure the zero average over the domain of q, since elements of L 2 cannot be prescribed on the boundary.

Discrete case
We note that existence and uniqueness in the primal formulation is given by the Lax-Milgram theorem. As a result, the use of a commuting diagram for the relation between the displacement u and the microdistortion P is to also be satisfied in the discrete case. Consequently, we rely on the commuting de Rham diagram for the construction of our finite elements. Specifically, for the lower order elements we make use of linear and quadratic tetrahedral Lagrangian elements, linear elements from the first and second Nédélec spaces and lowest-order Raviart-Thomas elements, see Fig. 3. We denote in the following On each tetrahedral element T we define P p (T ) as the set of all polynomials up to order p, P p (T ) = {x i y j z l i + j + l ≤ p, i, j, l ≥ 0} and denote with P p (Ω) the set of piece-wise polynomials of order p on the triangulation of Ω.
Remark 3.7. The Raviart-Thomas element from the space RT 0 can only produce constant normal projections on an element's outer surface. Further, the dimension of the space is smaller than the full linear space, The lowest-order Brezzi-Douglas-Marini elements BDM 1 have linear normal components as the full linear polynomial space is considered. On the one hand better L 2 -estimates are therefore possible, however, on the other hand the range of the divergence operator is not increased. There holds the sequence From Section 3.1.3 we know that the hyperstress field D ∈ H (Div, Ω) is solenoidal, Div D = 0. Therefore, aside from the lowest order case RT 0 , it is desirable to use a basis for RT p (or BDM p ) where all higher-order shape functions are divergence-free. We denote the adapted space by RT p 0 = {d ∈ RT p | d ∈ RT 0 or div d = 0} and the corresponding interpolation operator as Π p d . With such a construction, we can on the one hand save several redundant degrees of freedom, reducing the computational costs for assembly and solution steps without any loss of accuracy, and on the other hand, only a constant correction term q ∈ [P 0 (Ω)] 3 is required to enforce Div D = 0 for the remaining lowest-order RT 0 shape functions in the mixed formulation since In the numerical examples we use NGSolve for the high-order elements where such high-order divergence-free shape functions exist following the construction of [68].
Remark 3.8. The polynomial order of the sequence can be increased by using the Nédélec elements of the first type N 1 I instead of N 1 II , increasing the range of the curl differential operator. We note that both types yield a tangential projection of the same polynomial power on the outer boundaries of the element. Similar to the Raviart-Thomas and Brezzi-Douglas-Marini elements there holds In the Lax-Milgram setting we define X h = V h ×U h . Due to X h ⊂ X the solvability of the discrete problem follows directly from the continuous case [10,15]. We use Cea's lemma for the quasi-best approximation as a basis for the convergence estimate Figure 3: The de Rham complex for linear (upper) and quadratic (lower) sequences. The differential and interpolation operators commute between the continuous and the discrete spaces.
Remark 3.9. When using the canonical interpolation operators in the de Rham complex depicted in Figure 3 additional smoothness of the function spaces has to be assumed, as e.g., point evaluation is not well-defined for functions u ∈ H 1 (Ω), Ω ⊂ R d in dimension d > 1. First commuting, but non-local projections without additional regularity assumptions were constructed in [62] and [14]. Very recently, local L 2 -bounded and commuting projection operators from the function spaces into the finite element spaces have been established in [4]. Therefore, in the theoretical proofs we will make use of these novel projection operators, whereas for the finite element implementations the usual canonical operators are considered.
Proof. During the proof Π I,p−1 c denotes, with an abuse of notation, the projection operator without extra regularity assumptions, see Remark 3.9. With Cea's lemma and the approximation properties of the interpolation operators [11,15] we obtain where | · | H p denotes the standard Sobolev semi-norm.
When using the Nédélec elements of second type N p−1 II instead of N p−1 I , we lose one order of convergence for the microdistortion tensor due to the decreased range of the curl. This should also lead to sub-optimal convergence of the displacement. For example, the quadratic sequence depicted in Figure 3 would only have linear convergence for u although quadratic elements are used. In the extreme case of vanishing characteristic length L c = 0 the optimal convergence rates are achieved since only the L 2 -norm of P is considered, compare Section 3.1.2. On the subspace U 0 = {v ∈ H (curl, Ω) | curl v = 0} the Nédélec elements of first and second kind coincide, N p I = N p II , leading to optimal rates, which corresponds to the limit case L c = ∞. However, in the numerical experiments we continued to observe improved convergence rates for the L 2 -norm of P and H 1 -norm of u for L c ∈ (0, ∞). Under the assumption of an s-regular problem we can prove improved estimates for these norms with an Aubin-Nitsche technique. We call Problem 3.3 s-regular if there holds for the solution {u, P } ∈ X with s ∈ (0, 1] where Proof. We solve the following dual problem, where the differences u − u h and P − P h are used as right-hand side: find {w, Q} ∈ X such that By using the test-function Using Young's inequality on the left side, dividing through u − u h L 2 + P − P h L 2 , and using that {u − u h , P − P h } X ≤ ch p−1 for P h ∈ [N p−1 II ] 3 yields the claim for the L 2 -norm of u and P . For the H 1 -norm we first note that with the Galerkin-orthogonality, adding and subtracting Π p g u, and the specific structure of the bilinear form we obtain Thus, we can estimate ) where we used Cauchy-Schwarz and the explicit structure of the bilinear form. With the triangle inequality and positive definiteness of C e + C c (in combination with the generalized Korn's inequality for incompatible tensor fields if C c = 0) the claim follows.
Remark 3.10. We do not discuss the s-regularity of Problem 3.3 in this work, but mention that for the Laplace as well as for Maxwell's equations there holds s = 1 for convex domains Ω.
Corollary 3.1 also shows that we cannot expect improved convergence for u in the L 2 -norm compared to the H 1 -norm, which is confirmed by the numerics.
Unfortunately, the approximation constant depends a priori on L c . As we are interested in the limit L c → ∞ robust estimates with respect to L c are desirable. We show such an estimate with the help of the equivalent mixed formulation of the problem.
Normally, in contrast to the primal formulation, the kernel coercivity as well as the LBB-condition for the discrete case are not inherited from the continuous proof. However, thanks to the de Rham complex and a discrete Helmholtz decomposition [3], these properties hold also in the discrete case as long as the appropriate conforming finite element spaces are employed [11]. Consequently, Cea's lemma yields Additionally, with {u ∞ , P ∞ , D ∞ , q ∞ } for the (smooth) solution of the limit problem one obtains Proof. The proof of the first claim starts by following the same lines as in the primal case. The estimates for D and q are shown in detail where we used that the exact hyperstress D is divergence free. This estimate, however, would only yield linear convergence for all fields. By noting that on the subspace of divergence free hyperstresses the solution is independent of q we obtain the optimal rates for u h , P h , and D h . For the second inequality we add and subtract the exact solution, use the triangle inequality, and (3.29) and (3.40).
Remark 3.11. Note that only linear convergence is obtained for the correction term q. If one is interested in this quantity, a cheap post-processing step can be applied element-wise to obtain a high-order approximation [30]. For U h = [N p−1 II ] 3 the results of Lemma 3.2 can directly be adapted by following the same lines as in Corollary 3.1. Inequalities (3.41) show that besides the model error in L c , also the approximation error has to be considered for convergence studies in the limit L c → ∞.

Finite element formulations
In this work we employ tetrahedral finite elements. The construction of H 1 , H (div), H (curl), and L 2 -conforming finite elements to obtain the lowest order combination is presented in the following in detail. Higher order elements are used from the open source finite element software NETGEN/NGSolve 4 [61,63] and we refer to [68] for the construction of higher order finite elements. The elements are mapped from the reference element to the physical element using the barycentric base functions

Nédélec base
For the interpolation of the microdistortion tensor P we make use of Nédélec base functions of the first and second types [40,41,68]. The edge degrees of freedom are defined by the functionals where s i is the curve of an edge on the element and q j is a test function. In accordance with the de Rham complex, the appropriate polynomial space for the linear sequence is where P is the space of homogeneous polynomials. The corresponding test function q j = 1 yields the base functions on the reference element (see Fig. 6) representing the base of the lowest order Nédélec element of the first type N 0 I . For the sequence of the quadratic element we employ the linear Nédélec elements of the second type N 1 II . The corresponding polynomial space reads In other words, the space is given by 12 Nédélec functions. Using the ansatz 15) and the degrees of freedom from Eq. (4.11) for each edge on the reference element where µ i is the parameter for the corresponding edge, we find the base functions (see Fig. 7) Remark 4.1. The resulting Nédélec base functions are not unique in the sense that they are determined by the choice of the test functions q j . In our case we chose the linear Lagrangian base on µ, leading to a Nédélec-Lagrangian base.
Remark 4.2. Note that the decrease in polynomial order of the second discrete sequence (see Fig. 2) can be alleviated by employing first order Nédélec elements of the first type N 1 I . This increases the interpolation power of the solenoidal part of the microdistortion field P , but not of the irrotational part.
The base functions are defined on the reference element. In order to preserve their tangential projection on edges in the physical domain to find the covariant Piola transformation [38] θ, ds = θ, J dµ = ϑ, dµ ⇐⇒ θ = J −T ϑ . Further, for the transformation of the curl we find which is the so called contravariant Piola transformation. The formula is derived using the identity 5 ∇ x × J −T = 0. The microdistortion can now be interpolated by Θ P e , (4.22) In other words, the microdistortion is interpolated by 18 Nédélec base functions in the linear case and 36 in the quadratic case on each element.

Raviart-Thomas base
For the mixed formulation we also require Raviart-Thomas [53,68] elements. The lowest order face degrees of freedom of Raviart-Thomas elements are defined as The polynomial space for the construction of base functions is given by where P is the space of homogeneous polynomials. In the lowest order one finds Using the ansatz (4.28) and the degrees of freedom from Eq. (4.25) for each face on the reference tetrahedron where Λ denotes a face on the reference tetrahedron, we find the base functions (see Fig. 8) we consider the transformation of surfaces Consequently, the base functions are mapped using the contravariant Piola transformation Remark 4.3. Note that this the same transformation as in Eq. (4.21) since curl θ ∈ H (div) , θ ∈ H (curl).
Considering the divergence of functions undergoing a contravariant Piola transformation we observe As a result, one finds Thus, the hyperstress field is interpolated by such that, on each element 12 Raviart-Thomas base functions define the hyperstress.

Discontinuous basis
Finally, for the discontinuous elements in L 2 we employ piece-wise constants The entire space has the dimension dim[P 0 ] 3 = 3 and the interpolation reads

The orientation problem
The co-and contravariant Piola transformations do not suffice to assert the consistent orientation of the tangential or normal projections of the Nédélec and Raviart-Thomas base functions, respectively. The transformations control the size of the projections, but not whether these are parallel or anti-parallel with respect to neighbouring elements. Consistent projections is a key requirement in ensuring no jumps occur in the trace of the respective space and as such, there exist various methods for dealing with this so called orientation problem [1, 2,66,68]. In this work we present a solution based on the sequencing of vertices and the separation of orientational data. We define the following rule for the orientation of edges This means each edge starts at the lower vertex index and ends at the higher vertex index. This definition determines the orientation of the edge tangent vector (see Fig. 4) and consequently, the tangential projection of the Nédélec base functions. Analogously, for surfaces we define such that each surface is given by a sequence of increasing vertex indices. The orientation of the surface normal is given according to the left-hand rule. In other words, the direction of the normal is determined by the cross product of the vectors arising from the edges The orientation of the Nédélec-and Raviart-Thomas base functions is according to these rules. Consequently, in order to map each tetrahedron in the mesh to this orientation, we define each element as an increasing vertex-index sequence (see Fig. 9) The latter ensures the consistent orientation of the base functions, since they are all mapped from the same reference domain. However, integration in the reference element is determined by the determinant of the Jacobi matrix which may be negative due to a reflection of the element in the mapping from the reference to the physical domain. We correct for the error by taking only the absolute value of the determinant Consequently, consistency is guaranteed by mapping from a single reference element and the use of correction functions or considerations of neighbouring elements are circumvented.
Remark 4.5. The absolute value of det J is only used for the integration over the element. In all other use-cases, the information of the sign is necessary.

The discrete consistent coupling condition
In order to exactly satisfy the consistent coupling, one may use the degrees of freedom from [20] to set the boundary conditions, ensuring a commuting projection. However, this requires the solution of a Dirichlet type problem for every edge on Γ D . Alternatively, one can first find the Lagrangian interpolation of the boundary condition of the displacement on each edge and derive the boundary condition for the microdistortion directly from it. We observe that on each element's edge, the displacement is given by where n i are the Lagrangian base functions on the reference edge and m is the mid-point of the edge. The consistent coupling condition on a curve transforms according to the covariant Piola transformation where p represents a row of the microdistortion P and Π c p, ς results from the transformations of the tangent vector Eq. (4.19) and the Nédélec base functions Eq. (4.20). The gradient of the Lagrangian interpolation of the displacement Π g u can be reformulated in dependence of the curve parameter µ ∈ [0, 1] We note that in order to maintain an invariant transformation, the tangent vectors on the reference domain are given by and are not necessarily unit vectors, consider edges Fig. 4. Consequently, one finds (Π g u) ,µ D ξ µ, ς = (Π g u) ,µ D ξ µ, ξ ,µ = (Π g u) ,µ . (4.49) As a result, we can simplify the consistent coupling condition to For the linear element, each edge has one Nédélec base function of the first type with a constant tangential projection allowing to easily embed the consistent coupling condition Consequently, on each boundary edge the interpolation of one row of the microdistortion is given by In the quadratic sequence, each edge is associated with two Nédélec base functions of the second type ϑ v1 , ϑ v2 . Their projections on the tangent vector at their respective vertices is one and zero on the other vertex Consequently, it suffices to control the tangential projection of each base function at its respective vertex in order to impose the boundary condition (4.56) The resulting interpolation of one row of the microdistortion at the boundary edge reads For both sequences the above demonstrated methodology satisfies the discrete consistent coupling condition exactly by construction.

Element stiffness matrices
The matrices N , Θ, Φ, and 1 define the interpolation for each element. The resulting finite element stiffness matrices for the primal formulation read The linear element's stiffness matrix has the dimensions K e ∈ R 30×30 and the matrix of the quadratic element the dimensions K e ∈ R 66×66 . The load vector of the primal formulation is given by with the dimensions f e ∈ R 30 for the linear sequence and f e ∈ R 66 for the quadratic sequence, respectively. The mixed formulation induces the following additional stiffness matrices where K uu , K uP , and K P u remain unchanged with respect to the primal formulation. The linear sequence yields the matrix dimensions K e ∈ R 45×45 , whereas the quadratic sequence results in K e ∈ R 81×81 . The load vector changes its dimensions in the mixed formulation, but not its content with dimensions corresponding to the element stiffness matrix.
Remark 4.6. Note that the symmetry and anti-symmetry operators have been dropped. This is because the material tensors are defined as C e , C micro : R 3×3 → Sym(3) and C c : R 3×3 → so (3), implying the former. For materials where this is not guaranteed, the operators can be recovered by modifying the material tensors C e sym P = C e S P = C e P , C micro sym P = C micro S P = C micro P , C c skew P = C c A P = C c P . (4.62)

Numerical examples
In the following we test the convergence rates of the model. The behaviour with respect to the characteristic length L c is crucial, as this parameter governs the impact of the micro-structure on the displacement field u and the relation to the Cauchy continuum.
In order to study convergence we construct artificial analytical solutions for the displacement field u and the microdistortion field P (see Appendix A) and compare them with the numerical approximation of the model. Errors are measured in the L 2 -norm

Convergence for L c → 0
In this benchmark the domain is defined as the axis-symmetric cube Ω = [−1, 1] 3 . For simplicity the material constants in the isotropic case are set to λ e , µ e , µ c , λ micro , µ micro , µ macro = 1 .
The characteristic length is varied between L c ∈ (0, 1] in order to test the stability of the formulations. The prescribed fields are defined as being of higher polynomial order than the quadratic element. The entire boundary is set to ∂Ω = Γ D . By applying the strong form one finds the force and micro-moment Reducing the characteristic length L c or setting it to zero does not disturb the stability of the computation, even though control of the Curl of P is lost, as proven by the convergence rates in Fig. 10. This is a direct result of the derivation in Section 3.1.2, asserting stable computations for P ∈ [L 2 (Ω)] 3×3 .
In an alternative example we take the Cosserat couple modulus µ c = 0 (i.e., C c = 0), resulting in the force 6) and the micro-moment The primal formulation maintains its stability up to small L c -values coupled with small elements, see Fig. 11. For L c = 10 −10 the formulation becomes completely unstable and the matrix singular. The deterioration in the convergence is clearly depicted in Fig. 12. Since µ c = 0 and L c → 0, no term controls the skew-symmetric part of the microdistortion and the formulation only converges for the displacement u, up to L c = 10 −10 , where the stiffness matrix becomes singular.

Robustness in L c
To test for robustness of the formulation in the characteristic length L c we again consider Ω = [−1, 1] 3 , the parameters λ e = µ e = µ c = λ micro = µ micro = µ macro = 1 and µ c = 0, (5.8) and the following exact solution fields  where the first part of P is curl-free. The second term is scaled by L −2 c to obtain a well-defined limit for L c → ∞. The corresponding right-hand sides read We use NGSolve in combination with the direct solver UMFPACK [18] for inverting the arising stiffness matrix for the primal and mixed method. As can be clearly observed in Figure 13 for larger values of L c ≥ 10 6 the primal method becomes unstable, whereas the mixed method still converges optimally.

Convergence for L c → ∞
Next, we consider the convergence to the limit solution L c → ∞. Rewriting the exact solution (5.9) in the mixed formulation (3.28) and taking the limit L c → ∞ yields the following displacement, microdistortion, and hyperstress fields (5.10) Note, that Curl P = 0 and Div D = 0. The hyperstress field D acts as Lagrange multiplier forcing the microdistortion tensor to be curl-free. By taking the limit L c → ∞ we obtain directly the forces We consider five different structured tetrahedral grids consisting of 6, 48, 384, 3072, and 24576 elements and solve the primal and mixed formulation with NGSolve and UMFPACK using different values for L c . As depicted in Figure 14 the expected quadratic convergence rate in L c (3.29) is observed until the discretization error is reached (3.41). For finer grids we can reduce this error. However, the primal method again becomes unstable for large values of L c . The mixed method, where also the limit problem with L c = ∞ is well defined, is again completely stable with respect to L c .
Assuming micro to be always stiffer than macro, we set the macro and micro parameters λ macro = 115.4 , µ macro = 76.9 , λ micro = 10 λ macro = 1154 , µ micro = 10 µ macro = 769 ,  The behaviour of the deviation in both the displacement and the energy in Fig. 15 and Fig. 16 demonstrates the derivation in Eq. (3.14). Further, it is apparent that this characteristic strongly depends on the approximation capacity of the finite subspace. In other words, the use of p-refinement greatly influences the result. Decreasing the element size (h-refinement) can also be used to alleviate the error.

Bounded stiffness
In the example in Section 5.4 one can observe the influence of the characteristic length L c on the reaction of the system. In fact, the characteristic length acts as an interpolation parameter between the micro-stiffness and the macro-stiffness, as shown in Remark 3.5 and Section 3.1.2. In order to demonstrate this phenomenon, we employ the axis-symmetric cubic domain Ω = [−1, 1] 3 and the material parameters from the previous example. The characteristic length L c is varied in [10 −3 , 10 3 ]. The Dirichlet boundary is set to Note that while the problem is reminiscent of simple shear [57], it is different since we are using a finite volume and only the upper and lower parts of the boundary are constrained. The problem is defined entirely as a boundary value problem and no body forces or micro-moments are applied. In order to extract an indicator for the stiffness of the system we measure the Lebesgue integral of the displacement reaction forces in x-direction on the upper Dirichlet boundary where t x represents the equivalent traction in x-direction (see Appendix B). Considering Fig. 17, it is clear that the linear finite element formulation largely overestimates the stiffness of the system by about a factor of 28% in the lower limit for the finer discretization. Further, the effects of mesh refinement are considerable. From the depiction of the displacement field we recognize a transition from a nearly linear displacement field for large L c values to a higher order displacement function for low L c values. Lastly, we observe that the intensity of the microdistortion field (first row of P ) shifts from the centre to the sides of the domain for a decreasing characteristic length L c .
Remark 5.1. The lower limit C macro was estimated using 6000 finite elements of the quadratic sequence in the primal formulation while setting the characteristic length to L c = 0. In order to approximate the upper limit C micro , we make use of the quadratic sequence in the mixed formulation with L c = 10 9 and 6000 elements, since the mixed formulation is stable for very large values of the characteristic length.

Conclusions and outlook
Existence and uniqueness in the primal formulation of the relaxed micromorphic model requires the employment of Nédélec elements, since N p I , N p II ⊂ H (curl). The construction of the lowest order Nédélec elements was demonstrated along with a solution to the arising orientation problem, such that no correction functions are necessary and a single reference element suffices.
While not being a mixed formulation, the use of the exact de Rham sequence is still recommended in order to exactly satisfy the discrete consistent coupling condition, as shown in Section 4.6. This is otherwise not possible for the general case if the exact sequence is not respected.
The example in Section 5.4 depicted the behaviour of the relaxed micromorphic theory with L c → 0 as given by the derivation in Eq. (3.14). Clearly, the characteristic length L c is a crucial component of the theory as it governs the relation to the classical Cauchy continuum. Further, it is demonstrated that the latter comparison is strongly dependent on the approximation capacity of the finite subspace. Due to the relatively large error with respect to the expected result using the linear finite element sequence, u Cauchy − u relaxed L 2 / u Cauchy L 2 ≥ 9.5%, it is recommended to use at least quadratic or higher order finite elements. This conclusion is further supported by the example in Section 5.5, where the linear formulation largely overestimated the expected stiffness of the system for L c → ∞.
While the primal formulation is relatively robust in L c , it may become unstable as L c → ∞. As demonstrated in the example in Section 5.2, the computation can be re-stabilised by using a mixed formulation. As per the de Rham complex, this requires the usage of H (div) elements, such as Raviart-Thomas RT or Brezzi-Douglas-Marini BDM, and completely discontinuous finite elements. Further, for efficient simulations with higher order   Figure 17: (a) Reaction forces with varying characteristic lengths L c , (b) initial geometry with 3072 elements, (c) cutout view of the displacement field of the quadratic element with L c = 10 3 , (d) cutout view of the displacement field of the quadratic element with L c = 1, (d) cutout view of the displacement field of the quadratic element with L c = 10 −3 , (f) microdistortion field of the quadratic element with L c = 10 3 , (g) microdistortion field of the quadratic element with L c = 10 −3 .
elements it is advantageous to split the H (div) elements into source and solenoidal parts as the hyperstress field D is completely solenoidal, meaning the source part of the interpolant should be dropped or must be otherwise compensated. The construction of the lowest order Raviart-Thomas elements along with the lowest order discontinuous elements was demonstrated in Section 4.3 and Section 4.4, respectively. The orientation problem for the latter elements is circumvented by the same methodology as for the Nédélec elements. For further development we note the computation time increases rapidly with more elements and higher order polynomials. The development of an efficient preconditioner is to be considered for future works.

A Derivation of the strong form and analytical solutions
The variation of the free energy functional with respect to the displacement field u reads δ u I = Ω C e sym Dδu, sym(Du − P ) + C c skew Dδu, skew(Du − P ) − δu, f dX ∀ δu ∈ C 1 Γ D (Ω, R 3 ) , (A.1) where C 1 Γ D (Ω, R 3 ) denotes the set of differentiable functions which are zero at the Dirichlet boundary. By using the Green-type identity where P × ν on Γ P D is directly embedded into the space. Consequently, the strong form reads

B Equivalent traction
The traction on the Neumann boundary is defined using Eq. Remark B.1. Note that using the latter, no new computation of the stiffness matrix K, the nodal values vector u d or the force vector f d is required.


extfieldsofstudy: Biology Physics Medicine

text:
Population genetics in microchannels

Significance Many microbial populations proliferate in small channels. In such environments, reproducing cells organize in parallel lanes. Reproducing cells shift these lanes, potentially expelling other cells from the channel. In this paper, we combine theory and experiments to understand how these dynamics affects the diversity of a microbial population. We theoretically predict that genetic diversity is quickly lost along lanes of cells. Our experiments confirm that a population of proliferating Escherichia coli in a microchannel organizes into lanes of genetically identical cells within a few generations. Our findings elucidate the effect of lane formation on populations evolution, with potential applications ranging from microbial ecology in soil to dynamics of epithelial tissues in higher organisms.

Spatial constraints, such as rigid barriers, affect the dynamics of cell populations, potentially altering the course of natural evolution. In this paper, we investigate the population genetics of Escherichia coli proliferating in microchannels with open ends. Our analysis is based on a population model, in which reproducing cells shift entire lanes of cells toward the open ends of the channel. The model predicts that diversity is lost very rapidly within lanes but at a much slower pace among lanes. As a consequence, two mixed, neutral E. coli strains competing in a microchannel must organize into an ordered regular stripe pattern in the course of a few generations. These predictions are in quantitative agreement with our experiments. We also demonstrate that random mutations appearing in the middle of the channel are much more likely to reach fixation than those occurring elsewhere. Our results illustrate fundamental mechanisms of microbial evolution in spatially confined space. spatial population dynamics | bacterial evolution | microfluidics | individual-based models Biological populations can be spatially organized by landscape barriers that constrain individual movement, generating ordered patterns at the population level. For example, populations of rod-shaped bacteria Escherichia coli growing on surfaces tend to organize into domains of aligned cells (1,2). When E. coli grows in confined channels, cell alignment is affected by the geometry of the channel boundaries (3,4). In narrow channels, populations reach a highly ordered structure, with cells organized parallel to each other and to the boundaries (5,6). In wide channels, such alignment is disrupted at large scales by a buckling instability (7).
Once a cell population densely fills a microchannel, dividing cells push others toward the open ends, potentially leading to expulsion of cells. The timescale at which cells are expelled is typically shorter than their lifetime (an estimate of the latter is, e.g., in ref. 8). Therefore, death events can be usually neglected when focusing on a population inside a microchannel. In genetically diverse populations, we expect such competition to reduce diversity at a pace that depends on the channel dimensions relative to the cell size.
Microfluidic devices constitute ideal experimental systems to study population growth in confined geometries (9). Size and shape of the microchannels in these devices can be tailored to mimic microorganism habitats (10,11). Nutrients can be delivered to the residing microorganisms inside the microchannels by controlled flows. Such devices are often used to track population dynamics of microorganisms at the single-cell level over several generations (12,13). A paradigmatic example is the "mother machine"-a microchannel with one open end and small-enough width to accommodate a single lane of cells (8). Microchannels with two open ends have been used to validate a relation between the cell division time distribution and the population growth rate (14), but see (15).
From the theoretical side, competition in confined geometries has been scarcely studied. Common spatial competition models, ranging from Kimura's stepping stone model (16,17) to generalizations in the theory of evolutionary graphs (18)(19)(20)(21), include birth, death, and diffusion events only and do not account for cell to cell mechanical interactions in confined geometries. Single-lane models in which newborn individuals shift their neighbors away have been theoretically studied in the context of cancer progression (22,23). A microbial population dynamics model incorporating shifting has been investigated with computational approaches (14,24). In particular, simulations of a model of competing bacterial strains (24) show a formation of lanes along the channel axial direction. Theoretical predictions for geometries hosting multiple lanes and quantitative experimental validations have been lacking.
In this paper, we study the population genetics of bacterial colonies growing in confined geometries. We combine theory, numerical simulations, and experiments on E. coli populations growing in rectangular microchannels with two open ends. We introduce our work by first presenting experiments in which two fluorescently marked neutral E. coli strains competing in microchannels demix into a stripe pattern. We explain this phenomenon by means of an individual-based population model describing competition between strains inside the channel. This model reveals that the diversity loss within each

Significance
Many microbial populations proliferate in small channels. In such environments, reproducing cells organize in parallel lanes. Reproducing cells shift these lanes, potentially expelling other cells from the channel. In this paper, we combine theory and experiments to understand how these dynamics affects the diversity of a microbial population. We theoretically predict that genetic diversity is quickly lost along lanes of cells. Our experiments confirm that a population of proliferating Escherichia coli in a microchannel organizes into lanes of genetically identical cells within a few generations. Our findings elucidate the effect of lane formation on populations evolution, with potential applications ranging from microbial ecology in soil to dynamics of epithelial tissues in higher organisms. lane is much faster than that predicted by traditional spatial population genetics, in quantitative agreement with our experimental observations. We discuss the consequences of our results and their implications for the evolution of cell populations in confined geometries.

Competing Neutral E. coli Strains Form Stripe Patterns in
Microchannels. As a motivation, we present an experiment that anticipates a main consequence of our theory. We inoculated a mixture of two E. coli strains into microchannels (Fig. 1). The two strains are labeled with green fluorescent protein (GFP) and red fluorescent protein and are otherwise neutral (i.e., have the same fitness). Bacteria reproduce in nutrient-rich conditions inside the microchannels and push each other toward the open ends. As a result, cells are continuously expelled, while the number of cells inside the channels remains nearly constant. Our microchannels have a rectangular cross section, are 30-μm long, 1-μm deep, and have variable width from 1 to 3 μm, unless specified otherwise. For comparison, E. coli cells are 2.1 ± 0.2 μm in length and 0.65 ± 0.04 μm in width, so that the microchannels host monolayers of cells of width ranging from one to four lanes.
In about 8 h, the two strains organize themselves into regular stripe patterns ( Fig. 1 and Movies S1 and S2). The number of stripes and their width depend on the microchannel width and also fluctuate depending on the initial arrangement of inoculated cells. In these experiments, the average cell division time is ∼95 min, meaning that lanes are formed within a few generations. Our focus is on microchannels hosting monolayers of cells, but we also observe stripe formation in deeper (3-μm) microchannels, harboring multiple layers of cells (Movie S3).
Population Model Predicts the Stripe Pattern. We want to understand how the genetic diversity of a microbial population in a microchannel changes with time. We model the microchannel as a lattice of M × N sites. Each site is always occupied by one cell (Fig. 2A). We define clonal populations as groups of cells that originate from a common ancestor in the initial population. The dynamics of the clonal populations permits us to determine the patterns that the population would develop if some of the cells were fluorescently marked or carried a neutral mutation.
The dynamics proceeds as follows. Cells reproduce binarily at a constant rate b. After reproduction, one daughter cell takes the position of its mother. The other occupies one of the adjacent lattice sites and shifts a lane of existing cells toward one open end of the microchannel. As a result, a cell at the open end is expelled from the microchannel. If a reproducing cell is located next to an open end, its daughter can end up outside the microchannel, thereby being immediately expelled.
Our experiments reveal that the probabilities of choosing neighboring sites are not uniform. Specifically, we identify two effects that bias these probabilities. The first effect is related to the mass of the lane of cells to be shifted. We observe that reproduction events that shift shorter lanes of cells are more likely. We quantify this effect via a mass parameter m ≥ 0. Increasing m biases reproduction in the direction closer to an open end. The second effect is the preference of cells to reproduce within the same lane due to their aligned arrangement and the rod shape of E. coli. We introduce an alignment parameter α > 0 equal to the relative probability of a reproduction event within a lane over that of an event involving a change of lane. In the limiting case α = 1 and m = 0, the reproduction probabilities are uniform as in the model numerically studied in ref. 14.
We determine the parameters N, M, b, m, and α from our experiments with a single E. coli strain in channels of different width (Materials and Methods and SI Appendix). We find that N and b are the only parameters that significantly depend on the channel width (Table 1).
We take the number A(t) of clonal populations in the microchannel at time t as our measure of diversity. At the initial time t = 0, we have A(0) = MN . Diversity decreases with time, as progenies of initial individuals are expelled from the microchannel (Fig. 2B). At intermediate times, the surviving strains tend to form stripe patterns that resemble those in Fig. 1 (SI Appendix, Fig. S1 shows a more extensive qualitative comparison). The model further predicts that competition between stripes should lead to fixation of one of the strains at very long times.

Diversity Loss and Fixation Are Exponentially Fast in Sin-
gle-Lane Microchannels. We study diversity loss starting from microchannels with a single lane. In this case, we compute the rate of diversity loss by focusing on the interfaces between clonal populations (Materials and Methods). We find that the average diversity at time t is equal to where we define Eq. 1 shows that diversity loss in single-lane channels is exponentially fast. This result is in stark contrast with classic spatial population models, such as the voter model, where diversity decays as t −1/2 in one dimension (25). The characteristic rate β at which diversity is lost is on the order of the reproduction rate b, apart from a correction term that depends on the mass parameter m.
At long times, one clonal population eventually takes over the entire microchannel. The time at which this event occurs is called the fixation time. In microchannels with a single lane, the average fixation time is equal to where  Table 2).

First Regime of Diversity Loss: Exponentially Fast Fixation
within Each Lane. Our results in the single-lane case suggest that, in microchannels with multiple lanes, competition within each lane should lead to an exponentially fast diversity loss. In contrast, we expect competition among lanes to be less effective at reducing diversity. The alignment of cells favors reproduction events within each lane, further enhancing this difference. Following this idea, we identify two temporal regimes of diversity loss. In the first regime, diversity rapidly decreases from A = MN to A = M , primarily due to competition within lanes. The second regime ranges from A = M down to A = 1 and is characterized by competition among lanes. The first regime is characterized by negligible interaction among lanes. It follows from Eq. 1 that the average number of clonal populations at time t is approximated by We test this prediction in experiments with a single E. coli strain, where we track descendant of each individual in the initial population (Materials and Methods). We find an excellent agreement (Fig. 3A). Eq. 4 also implies that the quantity log[( A(t) − M )/(M (N − 1))] must be a universal linear function of βt (Fig. 3 A, Inset).
We approximate the average duration T MN →M of the first regime as (SI Appendix). This approximation and numerical simulations of the model agree well with our experiments (Table 2).

Second Regime of Diversity Loss: Slow Competition among
Lanes. In the second temporal regime of diversity loss, competition among lanes becomes relevant. This competition is driven by events in which a cell reproduces in a neighboring lane and its progeny eventually colonizes the entire lane. These events occur at a rate that we estimate to be quite small (SI Appendix). Aside from these events, lanes are typically dominated by a single clonal population (SI Appendix).
Following these ideas, we can consider lanes as single units, which invade each other at a certain rate. This process is called the invasion process in the literature (21). We mathematically solve this invasion process and thereby, estimate the probability to observe a given diversity A(t) in the second regime (SI Appendix). Our experimental results agree very well with simulations of the model and qualitatively agree with the results from the invasion process (Fig. 3B). Our theoretical and numerical results predict that the average fixation time for microchannels with multiple lanes is very long and therefore, inaccessible in our experiments ( Table 2).
Exponentially Fast Diversity Loss in the Mother Machine. We extend our theory to quantify the rate of diversity loss in a mother machine. Conceptually, the mother machine is similar to our microchannels with one lane. The main difference is that, in the mother machine, reproductions can occur in only one direction since one end of the microchannel is sealed. We solve our model with one lane under such conditions (M aterials and M ethods). In this case, we do not consider a mass effect, as the reproduction event can only occur in one direction. We find that, for the mother machine, the diversity loss is still given by Eq. 1 and the average fixation time is given by Eq. 3, where we set m = 0 in both expressions. These results show that the change in boundary conditions does not affect the dynamics of diversity loss.

Cells in the Middle of a Microchannel Possess a Positional
Advantage. We expect cells located far from the open ends of the microchannel to benefit from a positional advantage. We quantify this idea by means of the fixation probability P fix i,j , defined as the probability that the clonal population whose initial ancestor has coordinates i , j eventually takes over the microchannel. In the one-lane case and for large N, the fixation probability is approximated by where φ m,N (i ) = exp[−(i − μ) 2 /(2σ 2 )]/ √ 2πσ 2 is a Gaussian distribution with mean μ = (N − 1)/2 and variance σ 2 = (1 − m/2)(N − 1)/4 (SI Appendix). This means that, at increasing the mass effect, mutants that are likely to take over the population are located in a narrower region at the center of the microchannel (SI Appendix, Fig. S3G). In particular, the value of the mass parameter that we estimated (m = 0.6) leads to a 30% reduction in σ 2 , compared with the case m = 0.
In microchannels with multiple lanes, we approximate the fixation probabilities by ( Fig. 4A). The approximation in Eq. 7 is valid in the limit of large N as well. We also require the two regimes of diversity loss to be  well separated (SI Appendix). The fixation probabilities predicted by Eq. 7 are in good agreement with experimental observations (Fig. 4 B and C).

Discussion
In this paper, we studied the population genetics of microbial populations growing in microchannels with open ends. We base our study on a simple model in which cells are placed in lanes that are shifted by reproduction events. This shifting dynamics, combined with the geometry of the microchannels, causes a fast exponential loss of genetic diversity, rapid fixation within each lane, and slower competition among lanes. Our experiments quantitatively confirm these predictions and reveal that population dynamics generates stripes of clonal populations inside the microchannels. This outcome is in sharp contrast to the case of bacterial populations growing on agar plates, where competing populations organize into sectors whose boundaries perform random walks (26). Previous population models considered shifting dynamics in a single lane. Allen and Nowak (23) studied a similar model on a one-dimensional ring, concluding that this dynamics does not significantly affect the selection strength. However, a model of epithelial tissues has led to the conclusion that the shift dynamics suppresses selection (22). In this latter model, cells are arranged in a linear array and can be expelled only on one end, in a similar fashion as in the mother machine (8). Therefore, random mutations of the cell next to the opposite end are very likely to reach fixation, whereas fixation probabilities of mutations occurring elsewhere are very small. Our results show that this imbalance in fixation probabilities is a robust consequence of pushing dynamics in channels with multiple lanes and open ends, which does not require a strict geometric constraint imposing all cells to descend from a mother cell.
Our results can be extended to populations growing in wider and deeper channels. In such populations, the ordered lane structure is disrupted at large scales, potentially leading to jamming (27,28). In fact, crowding is known to affect the outcome of competition among microbial strains even in nonconfined geometries (29,30). Clarifying the connection between these evidences and our results is an interesting venue for future studies.
Our findings are potentially relevant for microbial ecology in soil. Bacteria in soil are spatially organized into relatively isolated and confined microenvironments (31), where fluid flows supply cells with nutrients (32). Most bacteria colonize micropores with diameter of about three times their body size (33) and in any case, smaller than 6 μm since small pores retain water for longer times. Moreover, bacteria residing in a pore of appropriate size are sheltered against larger predators (34). Although bacteria have very large population sizes, it was estimated that each bacterium in soil interacts with about 120 other individuals on average (35). Taken together, these observations support that the size of our microchannels is comparable with that of typical bacterial microenvironments in soil.
Renewing epithelial tissues in multicellular organisms (36) present a similar spatial organization as the one studied in our work. Two main examples are the epidermis (37) and the intestinal crypt (38). The epidermis consists of multiple compartments of cells that originate from a stem cell layer. These stem cells divide and generate differentiated cells that are shifted toward the top of the compartment and can be eliminated from the tissue once they reach the surface. In the intestine, stem cells divide at the bottom of each crypt, move upward, differentiate, and are removed once they reach the top of the villus. This dynamics permits rapid expulsion of cells that have accumulated deleterious mutations, thereby decreasing the risk of cancer, which would otherwise have high chances to arise in rapidly growing epithelial tissues (39,40). Similarly, proliferation of intestinal stem cells is thought to be disregulated in carcinogenesis (38). Our findings can potentially be extended to understand evolutionary dynamics of these tissues. In particular, adapting our model to study cancer dynamics would require introducing nonneutral clonal populations, as fitness differences between cell types are important for cancer progression (41,42). Moreover, cells that carry tumor-promoting mutations can be eliminated by apoptosis or killed by surrounding cells so that more detailed cellular interactions should be considered in this case.
The fact that a constrained geometry has such a drastic impact on population genetics should be taken into account when designing experimental evolutionary studies. Our results open possibilities for constraining evolution by shaping the geometry of a microchannel hosting a microbial population.

Materials and Methods
Bacteria Strains and Maintenance. MG1655, a derivative of the E. coli K-12 wild-type strain, was used in this study. We transformed a plasmid in MG1655 to constitutively express GFP (pUA66 PrpsL-GFP Kan R ) (43). For the red fluorescent strain, the pUA66 plasmid with mCherry open reading frame replacing that of GFP was customly constructed (VectorBuilder). We cultured the MG1655 hosting pUA66 plasmids in Luria-Bertani (LB) broth (Lennox) supplemented with 50 μg mL −1 kanamycin. A detailed description of bacteria culture and plasmid engineering is provided in SI Appendix.

Microfluidic Chip Design and Microfabrication Protocols.
The poly(dimethylsiloxane) (PDMS) device (on top of the cover glass) consisted of the top flow channels for nutrients delivery and bacteria removal and the bottom growth channels for bacteria growth, removal, and monitoring (SI Appendix, Fig. S2A). In order to perfuse fresh nutrients and create flow to remove bacteria, 20 growth channels of different width were intersected with 16 flow channels (Length × Width × Height [L × W × H] = 4500 × 50 × 15 μm). Growth channels (L × H = 30 × 1 μm) with varying widths (1, 1.5, 2.5, and 3 μm) and interspacing of 10 μm were designed such that bacteria could be expelled from both ends, in contrast to the single-end design of the mother machine (8,44). The flow channels were joined by flow-equalizing treelike channels and flow resistors on both ends (45) and connected to an inlet and an outlet.
A silicon mold with the microstructures designed as above was fabricated by multistep lithography with negative photoresist and maskless direct writing. PDMS microfluidic devices (Sylgard 184; Dow Corning) were fabricated by standard soft lithography (46). The PDMS slab was first cut and punched with an inlet and outlet using a puncher; then, it was bonded to a high-precision cover glass using plasma activation.A 2-mm-thick acrylic frame was cut with a CO 2 laser cutter and affixed on top of the PDMS as a reservoir for bacteria seeding, completing the fabrication of the integrated microdevice. Detailed information of design and microfabrication is provided in SI Appendix.
Bacterial Lineage Tracking and Time-Lapse Microscopy. Log-phase MG1655 E. coli harboring plasmids for fluorescent proteins was grown at 37 • C with vigorous shaking in LB broth supplemented with antibiotic kanamycin until the optical density at 600 nm reached 0.2. The bacteria suspension was concentrated 20 times by centrifugation before being injected into the PDMS microdevice, which was pretreated with a passivation solution to reduce bacterial binding to the microdevice surfaces. The PDMS microdevice was mounted in a microscope on-stage incubator preequilibrated at 37 • C on an inverted motorized epifluorescence microscope. The bacteria were allowed to enter the growth channels under static conditions for 2 h before the M9 media were infused using a syringe pump first at a flow rate of 1.6 μL min −1 , doubling every 2 h until it reached 16 μL min −1 .
The fluorescence images were taken with either a 100× oil immersion objective or a 60 × oil immersion objective with 1.5× intermediate magnification at an interval of Δt = 3 min using a high-sensitivity camera (Prime95B; Photometrics) with GFP or mCherry filter cubes. More details are provided in SI Appendix. Image Analysis and Data Processing. We processed the time-lapse recordings of our experiments using ImageJ software (47). We used the MicrobeJ plugin to detect bacteria in each frame (48) and a custom Python program to track all the bacteria in time. The tracking algorithm is based on construction and comparison of local structures for each cell (SI Appendix, Fig. S3A). We then reconstructed spatial lineage trees for each cell in the channel (SI Appendix, Fig. S3B).
In experiments with multiple lanes, microscopy focus drift issues for longterm live cell imaging experiments caused occasional quality loss in the recordings for a few frames. We cropped the recordings when such issues occurred. The durations of our recordings after this operation are 13.4 ± 4.7h, 9.7 ± 2h, and 7.93 ± 1.5h, corresponding to 5.64 ± 2.02, 5.84 ± 1.35, and 4.76 ± 0.95 generations for the experiments with two, three, and four lanes, respectively.
To explore the second regime, we perform a retracking of the experimental data (SI Appendix). In the retracking, we consider an initial condition, in which each lane is occupied by a single clonal population. Retracking is justified by the observation that, during the time at which A(t) = M in our original tracking, at least 70% of each lane is occupied by a single clonal population (SI Appendix).

Reproduction Rates and Model Parameters.
We assign to each cell its coordinates (i, j), with 1 ≤ i ≤ N and 1 ≤ j ≤ M. In the following, we refer to the coordinate i as the "horizontal" or "axial" coordinate. We define the reproduction probability k (i ,j )(i,j) as the probability that the daughter of a cell at position i, j is placed at position i , j .
We start from the case with one lane, M = 1. In this case, the probabilities p(i) = k i−1,i and q(i) = k i+1,i that a cell at position i reproduces to the left and right, respectively, are expressed by where m is the mass parameter. These probabilities depend linearly on i, as observed in experiments (Fig. 5D).
In the case M > 1, we determine the reproduction probabilities by imposing two constraints. First, the probabilities to divide to the left k (i−1,j)(i,j) and to the right k (i+1,j)(i,j) satisfy the condition k (i−1,j )(i,j) /k (i+1,j )(i,j) = p(i)/q(i) for j = (j − 1, j, j + 1). Second, the ratio between the probability of a cell division within a lane and that of a cell division involving a change of lane must be equal to α: k (i ,j)(i,j) /k (i ,j )(i,j) = α with i = (i − 1, i + 1) and j = (j − 1, j + 1).
The reproduction probabilities for cells in the bulk of the population (1 < j < M) satisfying these constraints are expressed by .

[9]
If i = i, the lane of cells j is shifted either to the left or the right with probabilities given by Eq. 8. We impose the same constraints for cells next to the top boundary of the microchannel (j = 1), obtaining . [10] The reproduction probabilities for cells next to the bottom boundary (j = M) can be similarly expressed. In all cases, if the i coordinate of the daughter is equal to zero or N + 1, she is immediately expelled from the microchannel. The model parameters evaluated from experiments are summarized in Table  1. For each microchannel width, we estimated M and N as the average number of lanes and the average number of cells per lane in our experiments, respectively. SI Appendix has details. The estimation of parameters b, m, and α is detailed in Fig. 5.

Dynamics of Interfaces.
We consider the case M = 1 and assign to each cell at position i at time t the position f i (t) of its ancestor at time t = 0. The quantity f i (t) changes every time the cell at position i is replaced by another one having a different initial ancestor. Two neighboring cells having the same value of f i (t) are conspecific (i.e., they belong to the same clonal population). We assign interfaces to neighboring cells that are not conspecific. We encode these interfaces into a vector σ(t) = (σ 1 (t), σ 2 (t), . . . , σ N−1 (t)), whose components are defined by The initial condition is σ(0) = (1, 1, . . . , 1). The vector of interfaces evolves until it reaches the absorbing state (0, 0, . . . , 0) that corresponds to fixation of one clonal population. Each cell division creates a pair of conspecific cells and shifts all cells by one position,either to their right or to their left.This event implies that one interface {σ i } i=1,...,N−1 is set to zero, and all interfaces on one side of it are shifted by one position, depending on the direction of the cell division. As a consequence, an interface located at the open end may be removed from the vector (Fig. 6). To describe the interface dynamics, we introduce inverse shift operators that take into account the presence or absence of an interface at the open end:

[12]
A stateb i r σ orĉ i r σ evolves to a state σ if the ith cell divides to the right. Similarly, a stateb i l σ orĉ i l σ evolves to a state σ if the (i + 1) th cell divides to the left. The master equation for the interface distribution is [q(i) + p(i + 1)]P σ (t)+ δ σ i ,0 [q(i)(Pb i r σ + Pĉi r σ ) + p(i + 1)(Pbi l σ + Pĉi l σ )], [13] where p(i) and q(i) are defined in Eq. 8 and the Kronecker delta takes care of the fact that a reproduction event necessarily creates an interface of value equal to zero. The solution of Eq. 13 reads where β is defined in Eq. 2. The solution given in Eq. 14 can be verified by direct substitution into Eq. 13. Eq. 14 shows that the interfaces {σ i (t)} i=1,...,N−1 are independent, identically distributed random variables with P(σ i (t) = 1) = e −βt for all i. The diversity is related to the number of interfaces by σ i (t) + 1.

[15]
Computing the average of A(t) using Eq. 14 leads to Eq. 1. Fixation Time in the Mother Machine. We apply the interface formalism to a model of the mother machine. Without loss of generality, we assume that the open end is the right one. Therefore, the vector of interfaces evolves according to a master equation including only right inverse shift operators: ( Eqs. 12 and 13). The solution to this master equation also factorizes ].

[17]
It follows that Eq. 3 with β = b holds for the mother machine as well.
Data Availability. All study data are included in the article and/or supporting information.