A Hybrid Particle-mesh Method for Incompressible Active Polar Viscous Gels

a r t i c l e i n f o a b s t r a c t We present a hybrid particle-mesh method for numerically solving the hydrodynamic equations of incompressible active polar viscous gels. These equations model the dynamics of polar active agents, embedded in a viscous medium, in which stresses are induced through constant consumption of energy. The numerical method is based on Lagrangian particles and staggered Cartesian finite-difference meshes. We show that the method is second-order and first-order accurate with respect to grid and time-step sizes, respectively. Using the present method, we simulate the hydrodynamics in rectangular geometries, of a passive liquid crystal, of an active polar film and of active gels with topological defects in polarization. We show the emergence of spontaneous flow due to Fréedericksz transition, and transformation in the nature of topological defects by tuning the activity of the system.


Introduction
Dynamics of many processes is governed by the mechanics of active matter. Active matter is a system of interacting agents that exhibit coordinated motion mediated by active internal stresses induced by consumption of energy [41,58]. Consumption of energy occurs at the level of individual agents keeping these systems persistently out of equilibrium [41,58]. Processes such as the collective motion of robots [60] and the motion of colloidal particles propelled by catalytic activity [50,41] are well described as active matter [41]. More striking examples of active matter are living systems like collections of bacteria [21,41], groups of biological cells [34,41], meshworks of cytoskeletal filaments [37,33,32], and flocks of birds [70,71] that exhibit mechanical motion by dissipating chemical energy. Such active systems constitute a class of physical processes where flow is driven from within by consumption of chemical energy.
Inspired by energy-fueled phenomena such as cortical cytoskeleton flows [46,45,32] during biological morphogenesis, the theory of active polar viscous gels has been developed [37,33]. The theory models the continuum, macroscopic mechanics of a collection of uniaxial active agents, embedded in a viscous bulk medium, in which internal stresses are induced due to dissipation of energy [41,58]. The energy-consuming uniaxial polar agents constituting the gel are modeled as unit vectors.
The average of unit vectors in a small local volume at each point defines the macroscopic directionality of the agents and is described by a polarization field. The polarization field is governed by an equation of motion accounting for energy consumption and for the strain rate in the fluid. The relationship between the strain rate and the stress in the fluid is provided by a constitutive equation that accounts for anisotropic, polar agents and consumption of energy. These equations, along with conservation of momentum, provide a continuum hydrodynamic description modeling active polar viscous gels as an energy consuming, anisotropic, non-Newtonian fluid [37,33,32,41]. The resulting partial differential equations governing the hydrodynamics of active polar viscous gels are, however, in general analytically intractable.
There have been attempts to numerically simulate the physics of active polar viscous gels at different length scales. At the microscopic scale, simulations have been carried out to understand in vitro organization of microtubule filaments in the presence of oligomeric motor proteins [49,47]. In addition, Brownian dynamics simulations have been used to model the collective dynamics of flexible cytoskeletal fibers [48]. However, such microscopic simulations do not allow large systems to be studied and are computationally demanding. At the mesoscopic scale, lattice Boltzmann simulations have been used to simulate the emergence of spontaneous flow in active liquid crystal films [42,43], spontaneous symmetry breaking in active liquid droplets [69] and pattern formation in active fluids [26]. While this increased the accessible length scales, it is challenging to impose constant polarity magnitude due to the limitation of the numerical method [7,42]. We are not aware of any numerical method to simulate at the macroscopic scale the hydrodynamic equations of active polar viscous gels proposed by Kruse et al. (2005) [37], Jülicher et al. (2007) [33] and, Joanny and Prost (2009) [32]. Lack of such a numerical method prohibited macroscopic simulations of complex active processes including coupled mechanochemical processes that have been postulated as the fundamental basis of biological morphogenesis [31,5,45,27,28].
In this paper, we present a hybrid particle-mesh method [29,30,40,9,63,10] to simulate the macroscopic hydrodynamics of incompressible active polar viscous gels on rectangular geometries with arbitrary boundary conditions. The method is based on Lagrangian particles and Cartesian meshes. In Lagrangian particle methods, field quantities are discretized onto a set of computational elements called particles [19,20,12,35,23,62,63]. Every particle occupies a position that evolves according to a local velocity field. The field quantities carried by particles evolve according to the material derivative in a Lagrangian frame of reference. Such a description forms the basis of numerical methods like smoothed particle hydrodynamics [52,8] and vortex methods [22,35,12] that have been used for various flow simulations. In such methods, advection of a field is performed by moving particles according to the velocity field. This enhances numerical stability in flow-dominated simulations [19,20,12,35,23,62,63] since the method is not restricted by the Eulerian Courant-Friedrichs-Lewy (CFL) condition [24]. Remeshing the field quantities from particles to meshes guarantees consistency in discretizing the advection operator [12,11] and allows using mesh-based methods for discretizing other differential operators. Mesh based methods, like finite differences and finite volumes, offer simple and accurate operator discretization schemes, especially in simple geometries. Hybrid particle-mesh methods therefore combine the advantages of Lagrangian particle methods and mesh based methods [29,30,62,63,3,10]. Using these advantages we develop a conservative numerical method where the novelty is the preservation of the constant magnitude of the polarization field. The numerical method uses an analytically derived Lagrange multiplier which is consistent with the physics of active polar viscous gels to conserve the magnitude of the polarization field. The Lagrange multiplier is a function of the velocity and polarization fields, and introduces additional coupling between the dynamics of the polarization and velocity fields. In addition, we perform moment conserving remeshing of the polarization field from particles to an underlying mesh that preserves the constant magnitude of the polarization field. We exactly impose the incompressibility constraint using staggered finite-difference meshes and direct linear solvers to obtain the exact solution of the resulting linear system of equations. Finally, we use particles to perform advection so that the discretization scheme is not necessarily restricted by the CFL condition for advection-dominated problems.
In Section 2, we present the governing equations of the continuum hydrodynamic description of incompressible active polar viscous gels in two-dimensional geometries. We present a hybrid particle-mesh method to numerically solve the governing equations in two-dimensional rectangular geometries in Section 3. In Sections 3.3 and 3.4, we present the accuracy and computational cost of the hybrid particle-mesh method. We validate our numerical method by demonstrating its consistency and convergence in Section 4. In Section 5, we use the hybrid particle-mesh method to simulate passive liquid crystal flows, spontaneous flows in active polar films, and hydrodynamics of active polar viscous gels with topological defects in the polarization field. We provide conclusions and discuss the limitations of the current method in Section 6.

Continuum hydrodynamic description of incompressible active polar viscous gels in two-dimensions
The continuum hydrodynamic description of incompressible active polar viscous gels in two-dimensions [37,33,32] (see Appendix A for more details) is governed by five partial differential equations (PDEs): one for each velocity component, one for the incompressibility condition, and one for each component of the polarization field.
Assuming that the viscous forces in an active gel dominate the inertial forces, the two equations governing the velocity field in the low-Reynolds-number limit are obtained by imposing force balance in each of the two directions: Here, p x (x, y, t) and p y (x, y, t) denote the x and y components of the polarization field p(x, y, t) at coordinates (x, y) at time t such that the magnitude |p(x, y, t)| is constant (or, without loss of generality, |p| = 1) .
denote the x and y components of the velocity field v(x, y, t), g ext x and g ext y are the x and y components of a externally applied force density g ext , is the pressure field, η is the viscosity, ν is coefficient of coupling between the mechanical stress and the polarization field, γ is the rotational viscosity, and ζ is the coefficient coupling the stress to the activity (or energy consumption) measured by the difference in chemical potential μ. The components u αβ of the symmetric strain-rate tensor u, the components σ (e) αβ of the Ericksen stress tensor σ (e) , and the transverse component h ⊥ of the molecular field h (h being the driving force tending to minimize distortions away from a spatially homogeneous polarization field) are defined as where K s and K b are the splay and bend elastic constants respectively. The last six terms on the left-hand side and the last two terms on the right-hand side in both of Eqs. (1) and (2) (1) and (2) are, therefore, a generalization of low-Reynolds-number flow equations for active polar viscous gels [37,33,32]. The third PDE is the incompressibility condition that is enforced by the pressure and is The last two PDEs govern the two components of the polarization field where Dp x Dt are the material (Lagrangian) derivatives of p x and p y respectively, γ is the rotational viscosity, and λ is a coefficient coupling the activity of the system (measured by μ) to the polarization dynamics. In addition, h || is the component of h parallel to the local polarization that enforces constant amplitude of the polarity vector (see Appendix A for more details), and ω αβ are the components of the antisymmetric vorticity tensor ω: Eqs. (10) and (11) reduce to the polarization dynamics of an incompressible passive liquid crystal when μ = 0. They hence are a generalization accounting for consumption of energy in active polar gels [37,33,32]. To solve Eqs. (1), (2), (9), (10), and (11) for the velocity and the polarization fields, we require an initial condition and various boundary conditions. The initial condition is given by the polarization field at time t = 0. The active gel equations can have a variety of boundary conditions. The velocity at the boundary is prescribed by a Robin (mixed) condition or a stress-free condition. Robin boundary conditions for the magnitude v of tangential and normal velocity components, on the one hand, are defined as a v + b (n · ∇)v = constant where a, b are constants, n is a unit vector pointing outwards at the boundary and ∇ is the gradient differential operator. On the other hand, velocity at stress-free boundaries are such that the shear stress at the boundary vanishes (see Eqs. (A.5), (A.6), (A.7) and (A.3) for the expression of the deviatoric stress tensor). Additionally, the polarization field at the boundaries is determined according to a Dirichlet or a Neumann boundary condition. Therefore, every point on the boundary is governed by four boundary conditions to define the two velocity components and the two components of polarization. In summary, Eqs. (1), (2), (9), (10), and (11), along with appropriate boundary conditions for velocity and polarization, govern the hydrodynamics of incompressible active polar viscous gels in two dimensions as a function of material constants, namely, η, ν, γ , ζ , λ, K s , and K b , and of the activity μ. We assume that the material constants are spatially homogeneous whereas μ ≡ μ(x, y, t) is, in general, spatially and temporally variable across the two-dimensional domain.

A hybrid particle-mesh method
We present a hybrid particle-mesh method to numerically solve the continuum hydrodynamic equations of incompressible active polar viscous gels (see Section 2) in two-dimensional rectangular geometries with arbitrary boundary conditions. The rectangular computational domain contains a set of points with coordinates (x, y) such that 0 ≤ x ≤ L x and 0 ≤ y ≤ L y , where L x and L y are the dimensions of the computational domain in x and y directions, respectively. All differential operators are discretized on staggered Cartesian finite-difference meshes, while the advection of the polarization field is represented on particles (for details on particle-function approximation see, for example, Hockney and Eastwood (1988) [30], and Cottet and Koumoutsakos (2000) [12]). We initialize particles carrying the local polarization field at mesh nodes at every time step. We integrate the material derivative and the local velocity field to compute the polarization field carried by the particles, and accordingly update the positions of the particles in the next time step. Subsequently, we interpolate the polarization field carried by the particles back to the mesh nodes. Using staggered finite differences, we then compute the velocity field at the next time step. This procedure combines the superior numerical stability of Lagrangian particle methods in performing advection of field quantities with the accuracy and simplicity of discretizing differential operators on Cartesian meshes to simulate the hydrodynamics of active gels.
Numerically simulating the hydrodynamic equations requires computing the velocity field v(x, y, t) and the polarization field p(x, y, t) across the computational domain as a function of time t. To do so, we first compute the initial velocity field at time t = 0 according to Eqs. (1), (2), and (9) from the initial condition for the polarization field. Subsequently, we integrate the equations of motion of the polarization field (Eqs. (10) and (11)) until time δt to compute the polarization field at t = δt, where δt is the integration time step. We then use the polarization at time t = δt to compute the velocity field at time t = δt according to Eqs. (1), (2), and (9). We use this procedure repeatedly to compute the polarization and velocity fields for all time steps until a final time t f = N t δt, N t being the total number of integration time steps.
The method, therefore, consists of two principal computational modules. One for integrating the polarization field from time t − δt to t according to Eqs. (10) and (11) using point particles and a finite difference mesh. We then make use of the polarization field at time t in the second computational module. The second computational module uses staggered finite difference meshes to compute the velocity field at time t. We make use of these two computational modules in the hybrid particle-mesh method to compute the velocity field and the polarization fields across the computational domain at all times.

Computing the polarization field
The computational module for numerically evaluating the polarization field is conceptually similar to solving for a concentration field governed by an advection-reaction-diffusion equation using particles and meshes. In the case of advectionreaction-diffusion, the diffusion and the reaction parts are evaluated on a mesh, and combined with the advection of particles [29,30,10,3].
To compute p(x, y, t), we use a finite-difference mesh M p and point particles that are initially placed at the nodes of M p . We assume the knowledge of velocity v , strain-rate u, vorticity ω, the polarization p, and the molecular field h at time t − δt on the nodes of M p . The nodes of M p are indexed as ( j, i) where j = 0, . . . , N x − 1 and i = 0, . . . , N y − 1. The x and y coordinates of node ( j, i) are j δx and i δ y, respectively, so that (N x − 1)δx = L x and (N y − 1)δ y = L y . Here δx and δ y are the spacing between nodes in the x and y directions, respectively, and N x and N y are the number of nodes in the x and y directions, respectively. The particle at the node ( j, i) is indexed with φ = i N x + j such that φ = 0, . . . , N x N y − 1. The particles are characterized by their positions x φ and carry the polarization field p φ such that x φ (t − δt) is the coordinate of node ( j, i) and p φ (t − δt) is the polarization at node ( j, i). We then integrate the equations of motion of the polarization field (Eqs. (10) and (11)) from time t − δt to t to compute the polarization on the particles at time t. We perform the integration according to the equations of motions of the particles in the Lagrangian frame of reference [30,44,12]: where g( j, i, t) denotes the value of any function g at the node ( j, i) at time t, and Dp Dt is the material (Lagrangian) derivative of p whose expression is given by the right hand side of Eqs. (10) and (11).
We numerically integrate Eqs. (15) and (16) from t − δt to t to obtain x φ (t) and p φ (t) for all particles. We integrate Eq. (15) from t − δt to t using the explicit Euler integration scheme with the initial condition given by x φ (t − δt), which is the coordinate of node ( j, i), at time t − δt. This advects the particles as illustrated in Fig. 1. The polarization field, however, can be purely oscillatory and using the explicit Euler scheme to integrate Eq. (16) will result in an unstable numerical simulation. As an instance of an integration scheme that can be stable in purely oscillatory regimes, we use the fourth-order Runge-Kutta scheme (RK4) [61,38] with the initial condition given by p φ (t − δt), which is the polarization on node ( j, i) at time t − δt. In order to remain consistent with the explicit Euler integration of Eq. (15), we let the velocity field, strain-rate tensor, and the vorticity tensor at the intermediate stages of RK4 be equal to that at time t − δt. The polarization field, however, changes at the intermediate stages of RK4. This requires recomputation of the material derivative (right hand side of Eq. (16)) at intermediate integration stages. We perform this computation on the nodes of M p by accounting for the boundary condition of the polarization field. This computation requires recomputing h || and h ⊥ at each integration stage.
We compute h || and h ⊥ on the nodes ( j, i) of M p according to Eqs. (13) and (8), respectively. The differential operators for computing h ⊥ are discretized as described in Section 3.1.2.
Once we computed x φ (t) and p φ (t), we interpolate the polarization field back to the nodes of M p . This procedure is generally referred to as remeshing [12] and serves to obtain p( j, i, t) at the nodes of M p for the next time step. We perform remeshing in order to enable discretization of differential operators on the mesh using simple finite difference schemes. The details of the remeshing scheme are described in Section 3.1.1. ) Particle-to-mesh interpolation using the 4,4 function. Illustration of interpolating to a node ( j, i) (denoted by a cross) from particles (hollow circles) that are within 3 mesh nodes from ( j, i). The interpolation weight (Eq. (18)) as a function of particle location with respect to node ( j, i) is color-coded. The figure illustrates interpolation of a field from 4 particles to node ( j, i). The interpolation weight (Eq. (18)) is a Cartesian product of the 4,4 function in each direction, so that the weights along the x-and y-directions can be computed independently and multiplied to obtain the final weight.

Remeshing
Remeshing of polarization p φ (t) carried by particles φ to the nodes ( j, i) of mesh M p is done using the 4,4 interpolation kernel [11]. The 4,4 interpolation kernel is a class C 4 function, conserves up to four moments of the interpolated field, and is fifth-order accurate if p φ (t) is computed analytically [11]. Accounting for the numerical error in explicit Euler time-integration of x φ (t) using v at t − δt (see Eq. (15)), remeshing using the 4,4 interpolation kernel is fourth-order accurate with respect to spatial discretization [11]. In addition, the computational cost of remeshing is linear in the number of target nodes. Due to its accuracy and favorable computational cost, we employ the 4,4 interpolation kernel.
We obtain the polarization field p( j, i, t) on the nodes of mesh M p as: where W φ ( j, i) is the weight (a multiplicative factor) of particle φ that contributes to the polarization field on node ( j, i).
This weight is: where x and ŷ are unit vectors pointing in the x and y directions, respectively, and 4,4 (·) is the interpolation kernel defined in Eq. (B.1). Because of the limited support of the 4,4 kernel (Eq. (B.1)), each particle contributes only to a bounded number of mesh nodes within its neighborhood. We, therefore, loop over all particles, find the mesh nodes that are within the support domain of the particle and add the particle's contribution to the those mesh nodes. Fig. 2 shows W φ ( j, i) and illustrates the interpolation from particles to a mesh node ( j, i) using the 4,4 interpolation kernel.
The 4,4 kernel is two-sided (symmetric) and for nodes that are within two mesh nodes of the boundary, the kernel does not have the necessary support domain. This leads to incorrect interpolation at these nodes. In order to faithfully interpolate polarization from particles to mesh nodes near the boundary, we impose polarization boundary conditions on the mesh M p through particles. For every particle φ at a distance d within three nodes of the boundary, carrying a polarization p φ , we introduce a ghost particle φ g at a distance −d from the boundary. When the polarization at the boundary is set to p b , the ghost particle φ g gets to carry polarization 2p b − p φ , which amounts to a second-order interpolation of the boundary value on the mesh. Similarly, we use second-order interpolation to compute the polarization carried by the ghost particle in the case of a Neumann boundary condition.
In order to preserve the unit magnitude of the polarization field upon remeshing, we also remesh the magnitude |p φ | of the polarization field from particles φ to nodes ( j, i) of mesh M p . The remeshed magnitude is used to renormalize the magnitude of the remeshed polarization field.

Computing the transverse component h ⊥ of the molecular field
When the polarization field p changes, the transverse component h ⊥ of the molecular field h needs to be recomputed. We need to perform this computation at all intermediate stages of the RK4 scheme used for numerically integrating Eq. (16), and to obtain h ⊥ at time t following the computation of p at the nodes ( j, i).
We compute h ⊥ on the nodes ( j, i) of M p according to Eq. (8). Computing h ⊥ involves computing four simple secondorder derivatives ∂ 2 x p x , ∂ 2 y p x , ∂ 2 x p y and ∂ 2 y p y , and two mixed second-order derivatives ∂ x ∂ y p x and ∂ x ∂ y p y (Eq. (8)). We numerically compute these derivatives with a second-order accurate spatial discretization. We compute the simple second derivatives according to the three-point central finite-difference stencil for all (interior) nodes that are not at the boundary.
For boundary nodes, we use a second-order one-sided four-point finite-difference stencil to compute the simple second derivatives. We compute the mixed derivatives by applying the first-order partial derivative in one direction followed by the first-order partial derivative in the other direction. We compute first-order derivatives on the nodes using a two-point second-order central finite-difference stencil for all interior nodes. For the nodes at the boundary, we compute the first-order partial derivatives using a second-order three-point one-sided finite-difference stencil.

Computing the velocity field
The computational module for numerically evaluating the velocity field is conceptually similar to finite volume [51,4,24] or staggered finite difference [51,4,24] approaches to solving time-independent partial differential equations like the Stokes equation.
We compute the velocity field v using mesh M p and three additional staggered finite-difference meshes M u , M v , and M . Mesh M u is shifted by δ y 2 in the y-direction with respect to M p . The nodes of M u are indexed as ( j, i − 1 2 ), where j = 0, . . . , N x − 1 and i = 0, . . . , N y (see Fig. 3 Nodes with i = 0 or i = N y are outside the computational domain and are referred to as ghost nodes (see Fig. 3(b)). All nodes except ghost nodes and boundary nodes are referred to as internal nodes. Mesh M v is shifted by δx Fig. 3(a)). Similar to the nodes of M u , nodes ( j, i − 1 2 ) indexed with i = 0 or i = N y − 1 and 1 ≤ j ≤ N x − 1 are boundary nodes, nodes with i = 0 or i = N y are ghost nodes (see Fig. 3(b)), and the rest are internal nodes. Mesh M is shifted by δx 2 and δ y 2 in the x and y directions respectively. The nodes ( j + 1 2 , i + 1 2 ) of M carry the pressure field , where j = 0, . . . , N x − 2 and i = 0, . . . , N y − 2 (see Fig. 3(a)). This set of staggered meshes M p , M u , M v , and M is used to numerically solve for the velocity v according to Eqs. (1), (2) and (9).
In order to numerically solve for the velocity v across the computational domain at time t, we assume the knowledge of p( j, i, t) and h ⊥ ( j, i, t) on mesh M p . We use these quantities in Eqs. (1) and (2), and along with Eq. (9), compute v x and v y (Eqs. (1) and (2)). We discretize Eq. (1) on the internal nodes of M u and evaluate all terms in the equation using a second-order central finite-difference for all differential operators. On the boundary nodes we evaluate a one-sided secondorder finite-difference scheme to enforce Robin boundary condition for the velocity component normal to the boundary. We use the ghost nodes and central second-order finite differences to enforce the Robin boundary condition for the velocity component tangential to the boundary. We use the ghost nodes to enforce the stress-free boundary condition using central second-order finite differences. Similarly, we discretize Eq. (2) on the nodes of M v . Additionally, we discretize Eq. (9) on the nodes of M , except on one node where we set to an arbitrary scalar. This pressure node acts as the reference pressure node since the gradient of the pressure field and not its absolute value is relevant. In our implementation, we choose to set N x 2 , N y 2 to zero. We solve the resulting unsymmetric, sparse system of N linear equations for the same number of variables using a direct method, the unsymmetric multifrontal method implemented in UMFPACK [15][16][17], to obtain v x on the nodes of M u , v y on the nodes of M v , and on nodes of M . Subsequently, we interpolate v x to the nodes of M p using second-order two-point central difference interpolation in the y-direction. Similarly, we interpolate v y to the nodes of M p using two-point central difference interpolation in the x-direction. We additionally evaluate the components of the strain-rate tensor u and of the vorticity tensor ω with secondorder accuracy on the nodes of M p . Using u, we numerically evaluate h || on the nodes of M p according to Eq. (13).
Combining all steps of the hybrid particle-mesh method, the complete algorithm, given the material constants η, ν, γ , ζ , λ, K s , and K b , the activity μ(x, y, t), the external force density g ext , and the boundary conditions for p and v, is as follows: : compute the polarization field p on the nodes of mesh M p using the initial condition for the polarization field.

6.2.
Compute p( j, i, t) using point particles and a finite-difference mesh: compute the polarization field p on the nodes , material constants γ , λ and ν, and the activity μ( j, i, t − δt), as described in Section 3.1.

6.3.
Compute h ⊥ ( j, i, t): compute h ⊥ on the nodes ( j, i) of M p using p( j, i, t), as described in Section 3.1.2.

Accuracy
The accuracy of the hybrid particle-mesh method is dictated by the numerical error of the time integration scheme used to compute the polarization field on particles, the numerical error in remeshing the polarization field from particles to mesh nodes, and the numerical error of discretizing differential operators on the meshes M p , M u , M v , and M .
The time integration scheme used here is globally first-order accurate due to the explicit Euler integration scheme used to advect particles (see Eq. (15)). In other words, the numerical error due to time integration is O(N −1 t ) where N t is the number of integration time steps. The remeshing from particles to mesh nodes is second-order accurate. This is mediated by the accuracy of imposing polarization boundary conditions on the mesh using ghost particles. The numerical error of remeshing using the 4,4 interpolation kernel is O(N −5 x + N −5 y ) [11]. Accounting for the explicit Euler time-integration of particle positions, the numerical error increases to O(N −4 x + N −4 y ) [11]. In addition, we impose boundary conditions for the remeshing procedure that are second-order accurate and hence the numerical error is O(N −2 x + N −2 y ). Accounting for all these errors, the overall numerical error of the remeshing procedure is O(N −2 x + N −2 y ), limited by the error in imposing boundary conditions. The discretization scheme used to evaluate differential operators on the meshes M p , M u , M v and M is second-order accurate. This is due to the use of one-sided second-order finite differences at boundary nodes and central second-order finite differences at internal nodes. The numerical error of discretizing differential operators is therefore O( In summary, the numerical error of the hybrid particle mesh method is O(N −2 x + N −2 y ) in the number of mesh nodes N x and N y in the x and y-directions, respectively. If the number of mesh nodes in each direction is equal to N m , the numerical error is O(N −2 m ). As a function of the number of integration time steps N t , the numerical error of the present hybrid particle-mesh method is O(N −1 t ).

Computational cost
The computational cost of calculating the polarization and velocity fields at each time point is mediated by the cost of each step in the algorithm. The scaling of the cost for computing v x , v y , and at mesh nodes is dominated by the unsymmetric multifrontal method [15][16][17] to solve for v x , v y and from the unsymmetric, sparse system of linear equations obtained by discretizing Eqs. (1), (2) and (9) on meshes M p , M u , M v , and M . By discretizing Eqs. (1), (2) and (9), we obtain N x (N y + 1) + N y (N x + 1) + (N x − 1)(N y − 1) linear equations involving v x , v y , and at the mesh nodes. We observe that the cost of solving for v x , v y , and at the nodes of M u , M v , and M , respectively using the unsymmetric multifrontal method [15][16][17] is O ((N x N y ) a ) where 1 < a < 2 with a ≈ 3 2 for the system considered here. The cost of computing h ⊥ and h || is O (N x N y ). The cost of the time integration scheme to compute the polarization field on particles is O (N x N y ). The cost of the remeshing procedure is O (N x N y ).
In summary, the computational cost of each time-step of the present implementation of our method is dominated by the cost of the direct method to solve the linear system of equations to compute the velocity components and pressure field, and is O ((N x N y ) 3 2 ). The overall computational cost is linear in the number of time steps, and a hybrid particle-mesh method can potentially use large time-step sizes since they are not restricted by the Eulerian CFL condition for advection.

Validation
We validate the hybrid particle-mesh method by studying its consistency and convergence on two test problems. On a test problem for which the analytical solutions of p(x, y, t) and v(x, y, t) are not known, we demonstrate consistency by showing that the error in the numerical solution decreases with increasing spatial and temporal resolution of the numerical simulation. A consistent numerical method, however, does not guarantee convergence of the numerical solution towards the correct solution of Eqs. (1), (2), (9), (10) and (11) upon increasing spatial and temporal resolution of the simulation. We demonstrate convergence of the numerical solution to the analytical solution at steady state on a second test problem where the steady-state analytical solutions of p(x, y) and v(x, y) can be computed. In addition, we show that convergence of the hybrid particle-mesh method is preserved for time steps corresponding to CFL numbers larger than 1.

Consistency
We show consistency of the hybrid particle-mesh method on a computational test case where the analytical solution is not known. We choose a computational domain with L x = L y = 10. The initial condition of the polarization field is a smooth harmonic map [1] p x (x, y, 0) = sin a, where a = 2 π (cos x 1 − sin y 1 ) (20) with We choose these initial conditions because they have been previously used to test accuracy of numerical methods for simulating passive liquid crystals [2,1]. The boundary conditions are where (x b , y b ) represent the coordinates of points on the boundary. The material constants are η = 1, ν = −0.5, γ = 0.1, ζ = 0.07, λ = 0.1, K s = 1, and K b = 1. The activity μ(x, y, t) = −1 and the external force density g ext = 0. Fig. 4 shows snapshots of the polarization, velocity, and pressure fields across the computational domain at times t = 0, 7.5, and 15 obtained with N x = N y = 17 and δt = 0.003. We observe that the system reaches steady state beyond t = 12 with non-zero steady state flow (Fig. 4(f)). We confirm the constant magnitude of the polarization field in Fig. 5, showing the maximum deviation of the magnitude of the polarization field from unity over the course of the simulation. We observe that the maximum deviation is close to machine epsilon, which is 10 −16 for double-precision floating-point arithmetic.
To measure the scaling of the numerical error due to spatial discretization, we perform simulations with increasing numbers of mesh nodes for a fixed δt. We perform simulations with N x = N y = N m for N m = 9, 17, 33, 65, and 129. We use δt = 2 × 10 −7 in all these simulations. Subsequently, we compute the error of the resulting numerical solutions at t = t f = 2 × 10 −6 against the numerical solution obtained from a more finely resolved computational domain with N x = N y = 257.
. See Section 4.1 for more details.  . This observation demonstrates that the present hybrid-particle mesh method is second-order accurate with respect to spatial discretization.
To measure the scaling of the numerical error with the number of integration time steps, we perform simulations with decreasing δt for a fixed N m . We perform simulations with δt = 4 × 10 −3 , 2 × 10 −3 , 1 × 10 −3 , 5 × 10 −4 , 2.5 × 10 −4 , 1.25 × 10 −4 , 6.25 × 10 −5 , and 3.125 × 10 −5 using a fixed number of mesh nodes in each direction N m = 17. We compute the error of the numerical solutions at t = t f = 1.024 against the numerical solution obtained with a finer integration time step δt = 1.5625 × 10 −5 . Fig. 6(b) shows the error in p x , p y , v x , and v y as a function of the number of integration time steps N t = t f δt . We observe that the error decreases with increasing N t and is O(N −1 t ). This demonstrates that the hybrid particle-mesh method is first-order accurate with respect to the time stepping scheme. In summary, the present hybrid particle-mesh method is consistent. In addition, the measured numerical error scales according to the theoretically expected bound (see Section 3.3) when increasing spatial or temporal resolution.

Convergence
We show convergence of the present hybrid particle-mesh method simulating a two-dimensional film where the steadystate analytical solution can be computed on a computational domain with L x = 10 and L y = 10. The initial condition for the polarization is The boundary conditions are where σ xy (x, L y , t) = 0 imposes a stress-free boundary along y = L y (see Eqs. (A.5), (A.6), (A.7) and (A.3) for the expression of σ xy ). We set η = 1, ν = 0, γ = 1, ζ = 1, λ = 1, K s = 1, K b = 1, and μ(x, y, t) = −1. The external force density is We perform numerical simulations until steady state is reached. We observe that by t = 40 the system reaches steady state. We, therefore, perform simulations with t f = 40 with increasing numbers of mesh nodes in each direction. We run simulations with N x = N y = N m for N m = 9, 17, 33, and 65. We run these simulations with δt = 0.04, 0.01, 0.0025, and 0.000625, respectively. In words, we decrease δt and increase the number of time steps N t = t f δt by a factor of 4 upon increasing the number of mesh cells N m − 1 by a factor of 2. More specifically, N m = 8N t 125 + 1. Fig. 7 shows the polarization and velocity fields as a function of time obtained by running a simulation with N m = 17. We observe that the velocity field is translationally invariant in the x-direction and v y (x, y, t) is on the order of machine epsilon. Fig. 8 shows the error in p x , p y , h ⊥ , and u xy from our numerical simulation at t = t f with respect to the steady-state analytical solution (Eqs. (C.11), (C.12), (C.13), (C.14), and (C.15)). We observe that the error decreases with increasing N m and is O(N −2 m ). Since N m = 8N t 125 + 1, we hence observe that the error is O(N −1 t ). In summary, the present hybrid particle-mesh method is consistent and converges to the analytical solution. Consequently, we conclude that the method is numerically stable. The error in the numerical solution decreases according to the theoretical error bounds (see Section 3.3).

Convergence for CFL numbers larger than 1
For a purely mesh-based approach with explicit numerical time integration, it is necessary that the CFL number |v x |δt δx + |v y |δt δ y ∞ is less than 1 for numerical stability [13,24], where the norm || · || ∞ denotes the maximum across the computational domain. This constraint limits the size of the time step δt. Hybrid particle-mesh methods have the advantage that the value of δt is not restricted by the CFL condition. Instead, stability of the advection operator in any hybrid particle- The magnitude of the pressure field is color coded. The parameters used in the simulation are L x = L y = 10, N x = N y = 17, δt = 0.01, η = 1, ν = 0, γ = 1, ζ = 1, λ = 1, K s = 1, K b = 1, μ = −1 and g ext = 0. See Section 4.2 for more details.  mesh method that uses explicit Euler for integrating particle positions and 4,4 remeshing requires that the Lagrangian CFL number δt ∇ v ∞ is less than 1 [11,9]. The Lagrangian CFL condition ensures that pathlines of particles do not intersect and is potentially less restrictive on δt. Here, we demonstrate the numerical accuracy of our hybrid particle-mesh method for CFL numbers larger than 1.
These time steps correspond to median CFL numbers of 0.32733 and 26.1866 respectively, where the median is computed across time steps. The corresponding median of the Lagrangian CFL numbers are 0.011354 and 0.90838, respectively. We observe that the steady-state polarization field at CFL number much larger than 1 agrees with the polarization field obtained using CFL number smaller than 1 (see Fig. 9). This shows that the hybrid particle-mesh method is accurate at CFL numbers larger than 1, thereby relaxing the constraint on δt. The CFL condition becomes more stringent with increasing number of mesh nodes N x and N y making the hybrid particle-mesh method more attractive for simulations with high spatial resolution.
In summary, the hybrid particle-mesh method potentially allows for time steps larger than those prescribed by the CFL condition. This potentially reduces the number of integration time steps and therefore the overall computational cost of simulating the dynamics of active polar viscous gels. When the CFL condition is fulfilled, or when the polarization dynamics is not advection dominated, the computational cost of our hybrid particle-mesh method becomes identical with that of a stable finite difference scheme.

Computational experiments
We report results from computational experiments on incompressible active polar viscous gels in rectangular geometries using the present hybrid particle-mesh method. The method enables simulations of the low-Reynolds-number hydrodynamics of simple fluids, of passive liquid crystal flows [65,39,1,2,14], and finally, of active polar viscous gels. Here, we present numerical results of the time-evolution of the polarization and velocity fields in a passive liquid crystal [1,2], in an active polar film [33,32], and in active polar viscous gels with topological defects in the polarization field [33,36].
We use a test case previously used to simulate passive liquid crystal flow [1,2]. We choose a computational domain with L x = 20 and L y = 10. The initial condition of the polarization field at interior points is given by Eq. (19). The polarization at the boundary is dictated by the boundary condition: for all times t. The velocities at the boundary follows a no-slip condition given by where (x b , y b ) denotes boundary points. We set η = 1, ν = 1, γ = 1, ζ = 0, λ = 0, K s = 1, K b = 1, μ(x, y, t) = 0, and g ext = 0. For the simulation we choose N x = 21, N y = 31, δt = 0.04, and t f = 500. Fig. 10 shows snapshots of the polarization and velocity fields at t = 0, 10.4, and 500. We observe that the polarization field relaxes to a homogeneous orientation across the computational domain (Figs. 10(a)-10(c)). At t = 500, the system is at steady state and the orientation of the polarization field is approximately 2.55 radians (or 146 • ). This direction of alignment is governed by the initial polarization field. Due to lack of gradients in the steady-state polarization field, we observe that the steady-state velocity field is zero ( Fig. 10(f)).

Active polar film
We show numerical simulation results of an active polar viscous gel of certain thickness with periodic boundary conditions in the other dimension. Such an active polar film is a model for in vitro actomyosin films and in vivo cortical actomyosin meshworks in biological cells [33,32].
We choose a computational domain with L x = L y = 10. The initial condition of the polarization field is The boundary conditions are p(x, 0, t) = p(x, 0, 0), where σ xy (x, L y , t) = 0 imposes a stress-free boundary along y = L y (see Eqs. (A.5), (A.6), (A.7) and (A.3) for the expression of σ xy ). We set η = 1, ν = −5, γ = 1, ζ = −4, λ = −1, K s = 1, K b = 1, μ(x, y, t) = −1, and g ext = 0. We choose N x = N y = 17 and simulate the system until t f = 160 with δt = 0.01. We observe that the system reaches steady state at around t = 90. Fig. 11 shows the polarization and velocity fields at times t = 0, 10, and 90. The polarization field is homogeneous along the x-direction (see Figs. 11(a)-11(c)). The polarization, however, is non-homogeneous along the y-direction and therefore the gradient of the polarization field is non-zero (see Figs. 11(a)-11(c)). The velocity field is also invariant along the x-direction, and the y-component of the velocity field is zero (see Figs. 11(d)-11(f)). Importantly, the velocity field at steady state is non-zero (see Fig. 11(f)). This non-zero steady-state flow is mediated by gradients in the steady-state polarization field, which in turn is imposed by the polarization boundary condition.
We next show a case of steady-state flow in active polar films due to gradients in the polarization field emerging spontaneously only when the ratio of activity μ to a critical activity μ c is greater than 1. When μ μ c < 1, we observe no steady-state gradients in the polarization field. The corresponding steady-state velocity field is, therefore, zero. When μ μ c > 1, we observe a non-zero steady-state velocity field. This transition to a non-zero flow velocity for μ μ c > 1 is referred to as the Fréedericksz transition [65,18,33,72]. To simulate a Fréedericksz transition, we choose a computational domain with L x = L y = 10. The initial polarization field is constant across the computational domain, except for a noisy perturbation: Eq. (27). We set η = 1, ν = −2, γ = 1, ζ = 1, λ = 0.1, K s = 1, K b = 1, and g ext = 0. For these parameters, the theoretical value of μ c = 0.2673 (see Eq. (C.18)). We perform simulations using N x = N y = 17 and δt = 0.05. Fig. 12 shows the polarization and velocity fields at t = 0, 100, and 300 for μ = 0.1 so that μ μ c < 1. The system reaches steady state at approximately t = 280. The polarization angle relaxes to π 2 across the computational domain (see Figs. 12(a)-12(c)). In other words, the polarization gradient at steady state is zero (see Fig. 12(c)). Consequently, the steadystate velocity field is zero (see Fig. 12(f)). We now repeat the simulation with μ = 0.4 so that μ μ c > 1. Fig. 13 shows the polarization and velocity fields at t = 0, 40, and 160. Beyond t = 150, the system is at steady state. The polarization angle at y = L y 2 relaxes to approximately 0.983 radians (or 56.3 • ) at steady state (Figs. 13(a)-13(c)). The gradient in the polarization field is, therefore, non-zero in the y-direction (see Fig. 13(c)). Consequently, the velocity field at steady state is non-zero (see Fig. 13(f)). This shows activity μ larger than the critical value μ c induces gradients in the polarization field leading to non-zero flow velocities.

Aster and spiral defects
We present simulation results of the dynamics of Eqs. (1), (2), (9), (10), and (11) in the presence of topological defects in the polarization field. At the topological defects, the polarization field is not defined (|p| = 0) and the gradient of the kinetic energy density (sum of squares of the two velocity components) is color coded. The kinetic energy density at t = 300 is almost 0, indicating that the magnitude of velocity is almost 0 everywhere. The lengths of the arrows representing the normalized velocity vector, however, are not zero, due to normalization with the maximum magnitude of velocity. The parameters used in the simulation are L x = L y = 10, N x = N y = 17, δt = 0.05, η = 1, ν = −2, γ = 1, ζ = 1, λ = 0.1, K s = 1, K b = 1, μ = 0.1 and g ext = 0. See Section 5.2 for more details. polarization field diverges. Such a point defect could be the center of an aster, a spiral, or a vortex. In an aster, the polarization field at each point is oriented along the line connecting the point with the defect. In a vortex, the polarization field at each point is orthogonal to the line connecting the point with the defect. Defects with any other constant orientation between the polarization field at each point and the line connecting the defect with that point are referred to as a spiral defects [33,36].
We present simulation results with spiral and aster defects in an incompressible active polar viscous gel. We perform a simulation with δt = 0.01 on a computational domain with L x = L y = 10 and N x = N y = 17. The initial polarization field represents a spiral defect (see Fig. 14(a)). The defect is at the center of the computational domain at x = L x 2 , y = L y 2 , where p = 0. In order to avoid a diverging gradient in the polarization field at the defect location, we smooth (regularize) the magnitude of the polarization field so that it increases smoothly to 1 over a distance of approximately 3 units from the defect. At the boundaries, we impose the velocities and the normal derivative of the polarization field to be zero. We set  The parameters used in the simulation are L x = L y = 10, N x = N y = 17, δt = 0.05, η = 1, ν = −2, γ = 1, ζ = 1, λ = 0.1, K s = 1, K b = 1, μ = 0.4 and g ext = 0. See Section 5.2 for more details. μ = 1.5 (Figs. 14(f) and 15(f)). The spiral character of the initial polarization field persists over time when μ = 1.5 (see Figs. 14(a)-14(c)). For the smaller μ = 0.5, the same initial polarization field gets transformed to an aster over time (see Figs. 15(a)-15(c)). This shows that tuning the activity μ can change a spiral defect to an aster defect.
We now switch the initial polarization field to an aster ( Fig. 16(a)) and perform simulations retaining all other settings. to a spiral defect. Such transformations of topological defects, for instance from a spiral to an aster or vice versa, have been predicted using linear perturbation analysis of the hydrodynamic equations [36]. Experimentally, such transformations have been observed in in vitro microtubule dynamics [49,66,67]. In these experiments, modulating μ by changing the concentration of motor proteins induced a change in the overall filament orientations from a spiral to an aster configuration and vice versa [49,66,67]. This is the first time such topological transitions are confirmed in a numerical simulations of the macroscopic governing equations of active polar viscous gels.

Conclusions and discussion
We have presented an algorithm to solve the macroscopic governing hydrodynamic equations of incompressible active polar viscous gels [37,33,32,41] in two-dimensional rectangular domains using a hybrid particle-mesh method. The method is a hybrid of Lagrangian particles and staggered finite-difference meshes. Combining the superior numerical stability of Lagrangian particle methods for flow-dominated problems with the simplicity of discretizing differential operators on staggered Cartesian meshes, we numerically solve the equations of active polar viscous gels. As the main novelty, the numerical method uses an analytically derived Lagrange multiplier to impose constant magnitude of the polarization field that is consistent with the physics of active polar viscous gels. We use particles to transport the polarization field so that the advection is not restricted by the Courant-Friedrichs-Lewy (CFL) condition. Staggered finite differences are used to exactly impose incompressibility, and a direct solver to obtain the exact solution of the resulting linear system of equations. Remeshing of the polarization field from particles to mesh nodes is performed using a moment-conserving interpolation kernel preserving the constant magnitude of the polarization field. Using these ingredients, we developed a novel conservative numerical method to accurately solve the hydrodynamic equations of incompressible active polar viscous gels. We have validated the hybrid particle-mesh method on two test problems. On one test problem, where the analytical solution is not known, we showed that the numerical error of the method decreases with increasing the spatial and temporal resolution of the simulation. On another test problem, where the analytical solution can be computed, we showed that the numerical solution converges to the correct analytical solution with increasing spatial and temporal resolution of the simulation. Specifically, we showed that the numerical error decreases by a factor of four upon doubling the number of mesh nodes in each direction and decreases by a factor of two upon halving the time-step size. These observations demonstrate that the presented hybrid particle-mesh method is consistent and converges to the correct solution of the underlying governing equations. Further, we demonstrated that the presented hybrid particle-mesh method is stable and converges to the correct solution even for CFL numbers much larger than 1 making the method advantageous especially for advection dominated dynamics.
Using the present hybrid particle-mesh method, we have performed computational experiments on a few instances of active polar viscous gels. Setting the energy consumption or the activity of the gel to zero, we simulated the hydrodynamics of passive liquid crystals using the same method. Using active polar films, which are a model system for cortical actomyosin layers in biological cells [33,32,72], we demonstrated the Fréedericksz transition [33,72,42] when the activity is larger than a critical value. At supra-critical activities, we showed the emergence of non-zero gradients in the steady-state polarization field. As a result we observed a non-zero velocity field. In active polar viscous gels with defects in the polarization field, we showed the transition from a spiral to an aster defect, and vice versa, by tuning the activity of the gel as observed in in vitro experiments on microtubule dynamics [49,66,67].
The present hybrid particle-mesh method has several limitations. The method is limited to two-dimensional rectangular domains and has a computational cost proportional to the total number of mesh nodes to the power 3 2 . Even though extending the method to three-dimensional geometries is straightforward, the scaling of the computational cost might prohibit efficient simulations in three dimensions. The computational cost of the method may, however, be reduced by using iterative linear solvers at the expense of numerical accuracy. For complex geometries, discretizing differential operators using simple finite difference schemes is problematic. This limitation can be addressed by using generalized finite-difference schemes, like discretization-corrected particle strength exchange (DC-PSE) operators [64,59]. In addition, the current time integration scheme is explicit and is therefore only conditionally stable. Using an implicit integration scheme can alleviate this limitation at the expense of larger computational cost. Future work will involve relaxing these limitations, leading to a robust hybrid particle-mesh method for simulating active polar viscous gels in complex three dimensional geometries. Additional future extension of the present method include coupling the local activity of active gels to concentration of chemicals that are advected by the flow velocity of the gel, diffuse, and undergo reactions that are deterministic or stochastic depending on the population of the chemicals [37,33,68,3,53,[55][56][57]54].
The theory of active polar viscous gels offers a framework for simulating active matter like migrating cells and cytoskeletal meshworks. The theory, in general, provides a basis for modeling active mechanochemical processes that are crucial to The parameters used in the simulation are L x = L y = 10, N x = N y = 17, δt = 0.01, η = 0.1, ν = 2, γ = 0.1, ζ = −1, λ = 1, K s = 1, K b = 1.5, μ = 1.5 and g ext = 0. See Section 5.3 for more details. See Section 5.3 for more details.
biological morphogenesis [31,5,27]. We envision an important role for computer simulations of such active mechanochemical processes in understanding complex phenomena in developmental biology.
αβ , , Eq. (A.1) is a force balance equation that is equivalent to the equation of conservation of momentum assuming negligible inertial forces (low Reynolds number regime). In Eq. (A.1), is the pressure, σ αβ denotes one of the four components of the deviatoric stress tensor σ , g ext β denotes one of the two components of a external force density vector g ext and ∂ β denotes the spatial derivative along direction β. The pressure enforces the incompressibility condition defined by Eq. (A.2). The deviatoric stress tensor σ with components σ αβ , is a sum of a symmetric stress tensor σ (s) with components σ (s) αβ , an antisymmetric stress tensor σ (a) with components σ (a) αβ , and the Ericksen stress tensor σ (e) with components σ (e) αβ (Eq. (A.5)). The Ericksen stress tensor is an equilibrium stress that generalizes the hydrostatic pressure to anisotropic fluids [18,25,6].
Eq. (A.3) is the constitutive equation for incompressible active polar viscous gels relating the strain-rate tensor u (with components u αβ ; see Eq. (A.10)), and the symmetric stress tensor σ (s) . The material constant η is the viscosity of the gel, h is the molecular field vector with two components h α , ζ is the coefficient coupling the material stress to the activity in the system as measured by the chemical potential difference μ, ν is the coefficient coupling mechanical stress to polarization field, and δ αβ are the components of the Kronecker-delta tensor such that δ αβ = 1 if α = β and 0 otherwise.
Eq. (A.3) models active polar viscous gels as a non-Newtonian fluid, since the stress is non-zero even when the spatial velocity gradients are zero [65]. γ is the rotational viscosity, λ is the coefficient coupling μ to the polarization dynamics, and ω αβ (Eq. (A.11)) are the components of the vorticity tensor ω. Eq. (A.7) defines the components σ (e) αβ of the Ericksen stress σ (e) in terms of a distortion free-energy density (or Franck free-energy density) f of polar nematic liquid crystals [18,65,33]. Eq. (A.6) defines the antisymmetric stress tensor σ (a) in terms of the molecular field h. The molecular field h is defined as the functional derivative of the volume integral of f with respect to the polarization vector p (Eq. (A.8)). The distortion free-energy density f defines the increase in the energy density due to distortions in polar nematic liquid crystals from its uniformly aligned configuration and is given by Eq. (A.9) where αβ are the components of the permutation (Levi-Civita) tensor. The free-energy density f is parametrized by the splay elastic constant K s and the bend elastic constant K b , neglecting the twist elastic constant since it is irrelevant in two dimensions. The free-energy density also includes a contribution from the parallel component h || of the molecular field assuming that fluctuations in polarity orientation dominate the fluctuations in polarity amplitudes. This assumption implies that the amplitude p γ p γ is a constant and can be assumed to be 1 without loss of generality [33]. In other words, h || acts as a Lagrange multiplier to ensure that the amplitude of the polarization field does not fluctuate. Using Eqs. (A.8), (A.9) and imposing p γ p γ = 1, h || = h x p x + h y p y .
(A. 13) The transverse component of h, h ⊥ , creates a torque that tends to align the polarization field. It is given by h ⊥ = p x h y − p y h x (A.14) . We consider a film that is infinitely long along the x-direction and has a thickness L in the y-direction. The surface of the film at y = L is stress free (σ xy (x, L) = 0) and impenetrable (v y (x, L, t) = 0). The surface of the film at y = 0 is also impenetrable (v y (x, L, t) = 0) and has a no-slip boundary condition (v x (x, L, t) = 0). The orientation of polarization at y = L and y = 0 is fixed to θ t and θ b , respectively, so that p x (x, L, t) = cos θ t , p y (x, L, t) = sin θ t and p x (x, 0, t) = cos θ b , p y (x, 0, t) = sin θ b .