{'fraction_non_alphanumeric': 0.040350444225074035, 'fraction_numerical': 0.013845014807502469, 'mean_word_length': 4.512311250170045, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dielectrowetting of a thin nematic liquid crystal layer

We consider a mathematical model that describes the flow of a Nematic Liquid Crystal (NLC) film placed on a flat substrate, across which a spatially-varying electric potential is applied. Due to their polar nature, NLC molecules interact with the (nonuniform) electric field generated, leading to instability of a flat film. Implementation of the long wave scaling leads to a partial differential equation that predicts the subsequent time evolution of the thin film. This equation is coupled to a boundary value problem that describes the interaction between the local molecular orientation of the NLC (the director field) and the electric potential. We investigate numerically the behavior of an initially flat film for a range of film heights and surface anchoring conditions.


I. INTRODUCTION
Dielectrowetting, a consequence of dielectrophoresis, involves the use of a nonuniform electric field to control spreading of dielectric fluids (such as nematic liquid crystals (NLCs), among others), on a substrate. A common focus of dielectrowetting experiments is the controlled manipulation of dielectric droplets using the interaction with an electric field to spread them into a thin film (see e.g. Brown et al. [3]); or vice versa: to break down a large parent droplet into smaller droplets (see e.g. McHale et al. [17]), with numerous possible applications [7,[21][22][23][24]27]. The classic dielectrowetting experiment, described in pioneering works such as those of Cheng et al. [7] and Brown et al. [3], involves a dielectric droplet or film spreading over a flat substrate, which contains interdigitated electrodes. A potential difference applied across these electrodes generates a nonuniform electric field, which acts on the polar molecules, driving flow. An excellent and comprehensive recent review of the field is provided by Edwards et al. [10]. A key feature of dielectrowetting that makes it desirable for applications is that it works with any dielectric liquid, in contrast to the more well-known "Electrowetting on Dielectric" (EWOD) technology, which relies on the movement of free charge and is limited to conducting liquids only. Dielectrowetting also offers the capability of spreading a partially-wetting droplet into a fullywetting film (and vice-versa) [16], which opens doors for new thin-film-based devices in the fields of microfluidics and optofluidics [10].
Theoretical investigations into dielectrowetting to date focus mostly on isotropic dielectric liquids [2,5,16,17]. Considering their wide industrial use, it is imperative to develop mathematical models that describe the behavior of anisotropic dielectric liquids, such as NLCs, in a similar setting. NLCs show rather complex behavior compared to Newtonian fluids because they consist of rodlike molecules, which have a dipole moment. In the absence of an electric field the molecules tend to align locally, due to interactions between the dipoles, which imparts elasticity to the material. Molecular alignment is also mediated by surface effects; socalled surface anchoring [6,9,26]. Applying an electric field also affects molecular alignment: the long axis of the NLC molecules will align parallel or perpendicular to the local electric field direction (depending on the dipole moment). Application of an electric field near a bounding surface sets up a competition between the dielectric force on the NLC molecules and the surface anchoring at the boundary. As a result, mathematical models of NLCs in such settings can be complex.
In this paper, we seek to develop a mathematical model that describes the flow of a thin NLC film under the effect of a nonuniform electric field. We consider a setup similar to that of Brown et al. [1,3] where a thin layer of dielectric liquid is placed on a substrate containing parallel interdigitated electrodes, which generate a periodic electric field profile. We make use of the commonly considered long wave scaling to derive a version of the thin film equation governing the evolution of the film height. This equation is coupled to a pair of boundary value problems, one for the electric potential within the film, and one describing the average orientation of the long axis of the NLC molecules (modeled in the theory by the director field n, a unit vector). The modeling is simplified by exploiting a disparity in timescales: both the director field and electric field are assumed to be in instantaneous equilibrium for the current free surface geometry. The resulting model is still complex, and numerical techniques are used to explore how the initial film height and anchoring conditions at the boundaries influence the evolution of the free surface over time.

II. MATHEMATICAL MODEL
In what follows, we develop a mathematical model that describes the flow and free surface evolution of a nematic liquid crystal (NLC) film under the effect of a nonuniform electric field. The basic setup, shown in Figure 1, consists of a thin NLC layer placed on a horizontal substrate. The substrate contains a pair of planar, interdigitated electrodes, leading to an electric potential that is nonuniform in the plane of the substrate. In our model we will assume that the electrodes are of infinite extent (in the y * direction in Fig. 1), so that the electric field generated by the applied potential will vary only in x * , the coordinate in the plane of the substrate perpendicular to the electrodes, and z * , which measures distance perpendicular to the substrate into the NLC. The main dependent variables that govern the dynamics are the velocity and director field of the NLC, together with the electric potential generated by the electrodes. In line with our assumed electrode geometry, throughout this paper we restrict attention to the 2D case in which the NLC velocity field v * = (u * , w * ) and the director field n = (sin θ, cos θ) are confined to the (x * , z * )-plane; so that the electric potential at time t * may be written as Ψ * (x * , z * , t * ); and the director orientation θ, which represents the angle that the director makes with the z * axis, is θ(x * , z * , t * ). Here and below, * superscripts are used to denote dimensional quantities. Quantities without superscripts are dimensionless.

A. Leslie-Ericksen Equations
The flow of NLC may be described by the standard Leslie-Ericksen model [11,13], which comprises the energy, (zero inertia) momentum balance and mass conservation equations as follows: where W * is the bulk free energy density and σ * is the stress tensor; these and the other quantities in Eqs. (1)-(3) are defined below: Here, subscripts i, j, k denote vector indices and the Einstein summation convention is assumed. The bulk free energy density W * consists of elastic and dielectric contributions. Nematic molecules prefer to align with their neighbors locally, a preference modeled by a bulk elastic (Frank) energy, W * e , where in Eq. (5) the widely-used one-constant approximation is used [9]. In addition, NLC molecules respond to an applied electric field E * , which induces a force causing them to align parallel or perpendicular to the field direction, giving rise to a dielectric free energy contribution modeled by W * d in Eq. (6). The constant ε * 0 is the permittivity of free space and ε and ε ⊥ are the relative dielectric permittivities parallel and perpendicular to the long axis of the nematic molecules. The molecular orientation induced by the field depends on the sign of the dielectric anisotropy, ε − ε ⊥ ; in line with the most common situation we mostly assume parallel orientation, associated with ε − ε ⊥ > 0. For our setup the electric field E * = ∇ * Ψ * is generated by applying a potential difference V * 0 across the electrode pair; for modeling simplicity we follow the approach of Brown et al. [3] and approximate the piecewise continuous substrate poten-tial by Ψ * (x * , 0, t * ) = V * 0 cos πx * 2d * , where d * is the electrode spacing. The parameters γ * 1 and γ * 2 in G * are constant viscosities, while e * ij and ω * ij are symmetric and antisymmetric rate-of-strain tensors for the material.
The stress tensor σ * for the NLC can be written as where σ e * , σ d * , σ v * are the elastic, dielectric and viscous contributions, respectively, and are defined as where ε(n) = ε ⊥ I + (ε − ε ⊥ )n ⊗ n is the dielectric tensor, and the α * i are constant viscosities related to γ * i in Eq. (7) by γ * 1 = α * 3 − α * 2 , γ * 2 = α * 6 − α * 5 and to each other by the Onsager relation, α * 2 + α * 3 = α * 6 − α * 5 . In addition to the electric field, NLC molecular orientation is sensitive to interactions with the bounding surface, a phenomenon known as anchoring. At a solid substrate, anchoring is determined by the chemical interactions between the NLC and the substrate. It is a common practice in applications to treat the substrate (chemically or mechanically) to impose strong planar anchoring with respect to the surface, therefore we assume a Dirichlet condition θ(x * , 0, t * ) = π/2 at the lower substrate. At a free surface, the director typically prefers to align normal to the surface (homeotropic anchoring), hence at z * = h * (x * , t * ) we impose weak homeotropic anchoring with associated anchoring strength A * , modeled by a Rapini-Papoular surface energy contribution, (where ν is the unit outward normal to the free surface), to the total energy [20]. In experiments, typical values for weak anchoring vary between 10 −5 − 10 −6 Jm −2 [19]. For the velocity v * , we assume no-slip and nopenetration conditions at the lower boundary, and a kinematic boundary condition together with a balance of normal and tangential stresses at the free surface.

B. Thin Film Model Derivation
We employ standard long wave theory scalings to non-dimensionalize the governing equations: x * = d * x, z * = δd * z, u * = U * u, w * = δU * w, Here, d * is the length scale of the electrode spacing along the x * -axis, U * is the typical flow speed in the x *direction and µ * = α * 4 /2 is the representative viscosity scaling in the pressure.
Experiments by Brown et al. [1,3,4] consist of a dielectric fluid film in which the typical film height h * 0 is much smaller than the electrode spacing (h * 0 ∼ 15µm and d * ∼ 120µm), hence we set δ = h * 0 /d * 1 to be the aspect ratio of the film. Typical values of δ range from 0.05 − 0.5 in the experiments [1,3]. If the free surface in the dimensional variables is given by z * = h * (x * , t * ) then we write h * = h * 0 h and the dimensionless free surface representation is z = h(x, t).

Energetics
Following an approach similar to that of Lin et al. [15] (which does not include electric field effects), provided the inverse Ericksen number N = K * /(µ * U * d * ) and the dielectric parameter D (defined in Eq. (15) below) are both O(1), the director energetics reduce to the Euler-Lagrange equations for minimizing total dimensionless free energy per unit length in the x-direction, given by where W is the total dimensionless bulk energy density and g(θ) is the Rapini-Papoular surface energy density, given in dimensionless form by Here, D represents the relative strength of the dielectric anisotropy and elasticity, γ = δ 2 d * 2 γ * /(K * h * 0 ) is the dimensionless surface tension coefficient defined in terms of the dimensional surface tension γ * and A = δA * d * /K * is the dimensionless anchoring strength at the free surface.
We follow the approach described by Cummings et al. [8] to minimize the total free energy density of the layer with respect to variations in θ and Ψ. The first variations must vanish at an extremum and the signs of the second variations indicate whether an energy minimum is reached. The bulk terms lead to the Euler-Lagrange equations, while the surface terms lead to one term that can be eliminated, giving the weak anchoring condition, The remaining boundary conditions on θ and Ψ are Dirichlet conditions: strong planar anchoring on z = 0, specified surface potential Ψ(x, 0, t), and we assume the free upper surface is electrically insulated (for simplicity), giving the following complete set of boundary conditions for θ and Ψ: Although both θ and Ψ are functions of the three independent variables x, z, t, in much of the following we will suppress the explicit t-dependence due to the quasistatic nature of the boundary value problem these functions satisfy. We will also, where convenient, suppress the independent variables altogether in θ and Ψ (likewise in h).

Momentum Equations
We use the pressure scale µ * U * /(δ 2 d * ) to nondimensionalize the stress tensor σ * . Then, provided that the inverse Ericksen number K * /(µ * U * d * ) = O(1), the momentum equations reduce to the condition that the stress tensor be divergence-free, and the equations governing the fluid flow can be extracted by retaining the leading order terms in where subscripts 1 and 3 refer to the x and z coordinate directions respectively, in our 2D geometry. As a result, the leading-order equations are found to be for the x and z components respectively, where Following the approach taken by Lin et al. [14,15], the normal and tangential stress balance conditions may be derived. The normal stress is balanced by the surface tension contribution in the Young-Laplace condition, which yields where C = δ 3 γ * /(µ * U * ) is an inverse capillary number. The tangential stress is balanced by the surface energy gradients along the free surface, which (after some algebra) reduces to the condition Equation (26) together with the normal stress condition of Eq. (28) gives the pressure p(x, t) throughout the layer as while Eq. (25) together with the tangential stress condition in Eq. (29) gives the velocity gradient across the layer as Imposing the no-slip and no-penetration boundary conditions, u = v = 0 at z = 0, and the kinematic boundary condition at z = h(x, t), together with mass conservation and the integral relation Substituting Eq. (31) into Eq. (32) above, we obtain a partial differential equation (PDE) that governs the evolution of the NLC film height h, where the functions T 1 (x, z), T 2 (x, h) and T 3 (x, h) take the following forms: In order to simulate a thin film of infinite lateral extent overlying a periodic array of electrodes, we impose periodic boundary conditions in the x-direction as follows: where the length of the computational domain L is taken to be an even integer (corresponding to simulating flow over an integer number of electrode pairs). We use L = 4 for all simulations presented in this paper, representing a single period-unit of the interdigitated electrode setup. Solving these equations simultaneously, even for the 2D geometry considered here, poses a significant computational challenge. Since the boundary-value problem (BVP) Eqs. (17)-(22) is quasistatic (reflecting the assumptions on the much slower time scale of the flow relative to the timescales on which NLC molecules and electric field respond to changes in film geometry) we are able to solve it independently for a dense grid of film heights h. To do this, we rewrite the BVP as a vector system of four first-order ordinary differential equations for θ, θ z , Ψ, Ψ z and apply the MATLAB routine bvp4c. An initial guess [θ 0 , Ψ 0 ] is required to start the routine; for this we use the results of the analogous uniform field problem, obtained as outlined by Mema et al. [8,18], for each height h. The procedure builds a "library" that consists of the director configuration θ(x, z) and electric potential Ψ(x, z) for a range of film heights 0.01 ≤ h ≤ 2.00 with ∆h = 0.01, as an electric potential Ψ(x) = cos π 2 x is applied at z = 0 for 0 ≤ x ≤ 4 with ∆x = 0.02 (the value of ∆z to compute the solutions θ, Ψ for each pair (x, h) is selected by the routine).
Once the BVP is solved for all discrete film heights h ∈ [0.01, 2.0] and all discrete x ∈ [0, 4], we solve Eq. (33) numerically. The integrals in Eq. (33) are difficult to evaluate directly; in the interests of retaining maximum tractability we follow Lin et al. [12,15] in using a two-point trapezoidal rule to approximate them and hence obtain the following fourth-order nonlinear partial differential equation for the film thickness h(x, t), The arguments of h are suppressed in Eq. (38) and those for T i are understood to be (x, z); the quantities are then given by Eqs. (34)-(36) respectively. We use a finite difference method based on central differences to discretize the spatial terms in Eq. (38) and use the built-in MATLAB routine ode15s, which is a variable-step variable-order solver [25], to approximate the height of the film at each time step. The procedure used to obtain the numerical results is as follows: given an initial film profile, h 0 (x) we solve Eq. (38) at each time step, extracting the director and electric potential configurations in the process as necessary from our previously built library.

III. RESULTS
In this section we use the above-derived model comprising Eq. (38) coupled with Eqs. (17)- (22) to investigate the flow of an initially flat film in the presence of a nonuniform electric field. We focus on the evolution of the free surface height h(x, t), as well as the director and electric fields within the layer, as we vary the initial film height, h 0 and the anchoring conditions at the upper bounding surface. With a view to later qualitative comparison with the experiments of Brown et al.  Figures 13-14 show film evolution for strong planar substrate anchoring and weak planar free surface anchoring of strength A = 10 and negative dielectric anisotropy D < 0. In all cases, the layer is subjected to a dimensionless electric potential Ψ(x, 0) = cos πx 2 on the substrate, and an insulating boundary condition Ψ(x, h) = 0 at the free surface. We set the inverse capillary number and Ericksen number to C = N = 0.625 and the relative strength of the dielectric anisotropy D = ±10 throughout. Our simulations are stopped either when the film height h falls below 0.1 or when a steady state is reached.
In many of our simulation plots, we superimpose representative snapshots of the director field as short directed line segments. We note that our model is invariant under the transformation θ → −θ (a limitation of our restriction to a 2D thin film model), hence a choice must be made when selecting θ. Although numerically our code selects values continuously in the range θ ∈ [0, π/2], our plots show line segment representations of the director that align visually with the electric potential lines (there is no mathematical distinction between the two possible states in our model).

A. Effect of initial film height on free surface evolution
We consider NLC films subject to strong planar substrate anchoring θ(x, 0) = π/2, and weak homeotropic anchoring at z = h(x, t) with dimensionless anchoring strength A = 10. We begin by demonstrating film stability in the absence of an applied electric field. As can be seen in Fig. 2(a), a slightly perturbed film of initial thickness h 0 (x) = 1.0 + 1 4π cos(πx) flattens in this case, demonstrating the stability of a flat film to perturbations. With no external field the NLC molecules align with their neighbors within the layer, while respecting the anchoring conditions at the boundaries.
We next consider a flat film, of initial dimensionless height h 0 = 1.0 (intermediate thickness), on application of the nonuniform electric field. In this case the film rapidly deforms, as seen in Fig. 2(b). The fluid collects around x = 1 and x = 3 (points of zero substrate potential) and begins to thin around x = 0, 2 and 4 (where the substrate potential takes values +1, −1, +1, respectively). In the presence of the nonuniform potential the polar nature of the NLC molecules induces a force that tends to align them parallel to the local electric field, competing with the forces due to the internal elasticity and the anchoring boundary conditions at z = 0, h respectively. As a result, molecules begin to migrate from the regions directly above the electrodes, x = 0, 2 and x = 4 where Ψ(x, 0) = ±1 and Ψ z (x, 0) is large, towards the inter-electrode regions above x = 1 and x = 3 where the electric potential Ψ(x, 0) and its gradient Ψ z (x, 0) are both small. This may be seen in Fig. 3, which illustrates the electric potential Ψ(x, z) and its level curves within the NLC layer at the same times shown in Fig. 2. Figure 4 shows the corresponding behavior of the director field: the director angle θ(x, z) and its level curves are plotted within the layer. We see that the director gradient |θ z | is large in the regions above the electrodes, x = 0, 2 and x = 4, and smaller in the regions between electrodes, at x = 1 and x = 3.
This behavior is qualitatively consistent with the ex- perimental results reported by Brown et al. [1,3,4] where a similar set up is considered: a layer of an isotropic dielectric fluid (1-decanol oil), rather than NLC, coats a glass substrate on which interdigitated ITO (indium tin oxide) electrodes are patterned. The authors observe that liquid flows from regions where the electric potential gradients are small to those regions where gradients are highest. We will briefly address ways in which the specific nematic nature of the film in our study may impact results later, when we consider how the details of the anchoring conditions affect the surface height evolution. Figure 4. Evolution of the director field θ(x, z) and its level curves for the film of Fig. 2(a).
We now investigate how changing the initial film height affects evolution. We first consider a thinner film of initial dimensionless thickness h 0 = 0.5, then a thicker film with h 0 = 1.5, both initially flat.  patterns, where the film thickens at x = 0.5, 1.5, 2.5 and x = 3.5 (representing the edges of the electrodes in our simple model); and thins at x = 0, 2, 4 (electrode midpoints) and at x = 1, 3 (points mid-way between adjacent electrodes). As time progresses the film continues to thin significantly at the electrode midpoints x = 0, 2, 4 where the potential gradient is large; these are the global film minima. Between the electrodes, where the potential gradient is small, further thinning of the film is suppressed and local minima persist at x = 1 and x = 3. We halt our simulations when the layer becomes too thin at the global minima and the director field no longer responds to the electric field in the interior of the layer. If the simulations are pushed further, solutions for the director configuration become unreliable. The evolution of the electric potential and the director field is qualitatively similar to the intermediate film height case: both director and electric potential gradients are large in the regions above the electrodes (x = 0, 2, 4) and small between them (x = 1 and x = 3). However, the gradients are larger for thinner films.    (Fig. 2), with the fluid collecting in regions where the electric potential gradients are small. As time progresses, however, the thicker film reaches a steady state, with approximately sinusoidal profile, exhibiting local minima of height h ≈ 1.42 at electrode midpoints x = 0, 2 and 4 and local maxima of height h ≈ 1.55 at x = 1 and x = 3, midway between electrodes. Our simulations stop when |h(x, t n ) − h(x, t n−1 )| < 10 −5 , as we deem a steady state to have been reached.

B. Effect of anchoring conditions on the evolution of surface height
We next investigate the influence of different anchoring conditions on the evolution of an initially flat film. Specifically, we first study how, in the model described above, the strength of the weak homeotropic anchoring at the free surface affects the surface height evolution. Then, with a view to conducting simulations that are closer in spirit to the isotropic dielectric case studied by Brown and co-workers [1][2][3], we briefly study the evolution of NLC films with negative dielectric anisotropy, subject to planar anchoring at both boundaries.

Weak homeotropic anchoring strength
We consider the model described above, comprising Eq. (38) coupled with Eqs. (17)- (22), with strong planar anchoring at the lower substrate and weak homeotropic anchoring at the upper free surface. All previous simulations were for dimensionless (Rapini-Papoular) anchoring strength A = 10; we now simulate initially flat films where A = 1 and A = 50 at the free surface. We again compare results for three different initial film thicknesses, h 0 = 0.5, 1.0 and 1.5. Figures 7 and 8 show the evolution of an initially flat film of height h 0 = 1.0 and the corresponding electric potential Ψ(x, z) and director field θ(x, z) at the final computed time for anchoring strengths A = 1 and A = 50 respectively. The evolution of the corresponding film with homeotropic anchoring of dimensionless strength A = 10 at the free surface was shown in Fig. 2(a). Similar to the A = 10 case, our simulations stop when the layer becomes too thin and the solution for the director field becomes unreliable. While the evolution of the film is similar for all three anchoring strengths considered, there are some differences: first, in the weaker anchoring case (dimensionless anchoring strength A = 1), the director at the free surface deviates significantly from the preferred homeotropic orientation (see Fig. 7(b)); this deviation is especially apparent in regions between the electrodes (where the electric potential is close to zero). For (relatively) strong anchoring A = 50, by contrast, the director almost perfectly respects the preferred anchoring conditions at the free surface: see Fig. 8(b). Additionally, we note that the maxima in the free surface height profile at x = 1 and x = 3 become more diffuse as the anchoring strength increases (peaks shown in Fig. 8(a) are wider than those in Fig. 7(a)).  We now discuss the effect of anchoring strength on a thinner film. Figures 9 and 10 illustrate the evolu- tion for a film of initial height h 0 = 0.5, with weak homeotropic free surface anchoring of strengths A = 1 and A = 50 respectively, together with the corresponding electric potential and director field at the final computed time. The corresponding evolution with A = 10 is shown in Fig. 5. As in that figure, we observe here that the film surface exhibits wrinkle formation and develops a distinctly non-sinusoidal profile. The formation of wrinkling patterns appears to depend primarily on the initial film height and not on free surface anchoring strength. We expand on this conjecture in the next section when we discuss how planar anchoring conditions on both boundaries affect the surface height evolution.
There are, however, some qualitative differences in the evolution as anchoring strength varies. First, for weaker anchoring strength and in the regions between electrodes where the electric potential is near zero, the film height slowly increases in time (see Fig. 9(a) at x = 1 and x = 3). By contrast, as anchoring strength increases to A = 10 and A = 50, at the same locations the film height decreases (see Fig. 5(a) and Fig. 10(a)). Moreover, the director orientation at the free surface deviates significantly from the preferred orientation when A = 1; in fact for this weakest anchoring strength the director orientation is planar throughout the film in the regions between electrodes, where the electric field is weakest: θ = π 2 is a stable (constant) steady solution to the governing equations (cf. Ref. [19]). The strong planar substrate anchoring dominates over both the weak homeotropic anchoring and the weak electric field in regions between the electrodes, forcing the director to align parallel to the substrate throughout the layer: θ(1, z) = θ(3, z) = π/2.
As the anchoring strength increases to A = 10 and A = 50 the constant stable steady state ceases to exist, even in the weak-field regions; the director bends across the layer at all x-locations, deviating only slightly from the preferred homeotropic orientation at z = h. In regions above the electrodes (assuming positive dielectric anisotropy), the electric field helps the director align closely with the free surface homeotropic anchoring orientation, independently of the anchoring strength.
Finally, we comment briefly on the effect of anchoring strength for thicker films: Figures 11 and 12 illustrate the evolution of an initially flat film of thickness h 0 = 1.5 with strong planar substrate anchoring and weak homeotropic free surface anchoring, of strengths A = 1 and A = 50 (respectively). The analogous evolution with A = 10 is shown in Fig. 6(a): similar to that case, both films reach a steady-state near-sinusoidal profile (shown in Figs 11(a) and 12(a)). Small differences in the free surface height are observed for the different cases: as anchoring strength increases, the amplitude of the perturbations to the steady-state surface height decreases, and the film height maxima become more diffuse. In all cases the maxima are achieved at the points mid-way between the electrodes, x = 1 and x = 3. For the smallest anchoring strength A = 1 the maximum film height at the final time t f , h(x, t f ) ≈ 1.6. As anchoring strength increases to A = 10 and A = 50, the maximum film height decreases slightly: h(x, t f ) ≈ 1.55. We note also that the director at the free surface deviates slightly from the preferred homeotropic orientation at the weakest anchoring strength A = 1 when compared to the simulations with stronger anchoring strengths A = 10 and A = 50. As expected, however, this difference is less pronounced than for thinner films. field, we now discuss the evolution of an initially flat NLC film with negative dielectric anisotropy, subject to strong planar anchoring at the lower substrate and weak planar anchoring of strength A = 10 at the free surface. With dielectric anisotropy parameter D < 0 the nematic molecules align perpendicular to the electric field, rather than parallel to it, so that all external effects here favor a director field that orients parallel to the substrate, with θ = π/2 throughout and zero bulk distortion (the elastic contribution W e = θ 2 z /2 to the bulk free energy density in Eq. (15) is zero). We present results for three different initial film thicknesses, h 0 = 0.3, 1.0, 1.5 and compare them with the analogous simulations presented above (with strong planar anchoring at the substrate and weak homeotropic anchoring at the free surface, which we refer to as the planar-homeotropic case). Here, we choose the smaller value h 0 = 0.3 as representative of a thin film (rather than h 0 = 0.5 considered earlier) in order to observe free-surface wrinkling patterns, which are not exhibited when h 0 = 0.5. Figure 13 shows the evolution of an initially flat film of intermediate height h 0 = 1.0. As noted above, the director aligns parallel to the substrate throughout the layer (see Fig. 13(b)), which illustrates the director configuration superimposed on the electric field potential at the final computed time. There are similarities to the evolution of the initially flat planar-homeotropic film of the same height shown in Fig. 2(a): the film thickens around the points x = 1 and x = 3 mid-way be-tween the electrodes where electric field gradients are small, and begins to thin directly above the electrode midpoints at x = 0, 2 and 4, where electric field gradients are large. There are some differences, however: the surface height profile in Fig. 13 has well-defined peaks exactly at x = 1 and x = 3; while the surface height profile for the planar-homeotropic case in Fig. 2 has more diffuse maxima. We attribute the diffuseness to the competing forces that the director experiences in the planar-homeotropic case. This competition becomes more intense as free surface anchoring strength is increased, hence the film maxima are even more diffuse at the largest anchoring strength considered (A = 50 in Fig. 8(a)).
We next discuss the thin film of initial height h 0 = 0.3. The evolution, shown in Fig. 14, is compared to that of the flat planar-homeotropic film of height h 0 = 0.5 as shown in Fig. 5(a) (see also Figs. 9(a) and 10(a) for different free surface anchoring strengths). Complex, strongly non-sinusoidal wrinkling patterns emerge in all cases, with the films thinning/thickening in regions of large/small electric field gradients. We here highlight key differences between the cases: first, note that in Fig. 14 the fluid collects in regions midway between the electrodes where the electric potential gradient is small, forming local film maxima at x = 1 and x = 3 (local minima are formed at these points for the planar-homeotropic case in Figs. 5, 9(a) and 10(a)). Second, we observe four small additional peaks in regions of high electric potential gradients, near x = 0, x = 2 (one on either side) and x = 4 that were not observed in the earlier simulations (where surface height decreases over time at these locations). Despite these differences, the re-emergence of wrinkling patterns for these different anchoring conditions reinforces our conjecture that they are ubiquitous for sufficiently thin films (though the details of the anchoring and dielectric anisotropy may modulate the specific wrinkling patterns observed).
Finally, we compare the evolution of a thick film, h 0 = 1.5, to that observed in the planar-homeotropic case shown previously in Fig. 6 (see also Figs. 11 and 12 for different free surface anchoring strengths). Figure 15 shows the evolution of a film with strong planar substrate anchoring, weak planar free surface anchoring of strength A = 10, and negative dielectric anisotropy. As in the planar-homeotropic cases of Figs. 6, 11 and 12 we observe that the film thickens in regions between electrodes where electric field gradients are small and thins directly above electrodes, where electric field gradients are large. In addition, and again similar to the planarhomeotropic cases, a steady state is quickly reached. The surface height profile in Fig. 15 has well-defined maxima at x = 1 and x = 3 while the corresponding planar-homeotropic surface height profile in Fig. 6 exhibits more diffuse maxima at these points. As with the intermediate film thickness h 0 = 1 simulations, we conjecture that this difference is due to the competition between different forces on the NLC molecules in the planar-homeotropic case.
Finally, we note that from a modeling perspective we anticipated that this last-considered case of NLC film evolution with planar anchoring conditions and negative dielectric anisotropy should provide results most analogous to those for an isotropic dielectric liquid, due to the uniform director orientation throughout the layer. However, when comparing our results with the experimental observations reported by Brown et al. [3], we see that a closer qualitative match is obtained with the earlier planar-homeotropic simulations of Figs. 2-12 than with the results of this section. Since a NLC, even with uniform director orientation, is not the same as an isotropic dielectric, we conclude that other factors must account for the observed differences in evolution. We believe that the stronger qualitative similar- ities between the experiments and our simulations for the planar-homeotropic case are at least partly coincidental.

IV. CONCLUSION
We have presented a mathematical model that describes the flow of a thin NLC film in the presence of a nonuniform electric field. Specifically, we consider a thin layer of NLC coated on a substrate z = 0 that contains embedded planar interlaced electrodes, the effect of which we approximate by a periodic electric field profile on z = 0 in the x-coordinate direction. The mathematical model derived consists of a 4th order nonlinear parabolic partial differential equation for the height h(x, t) of the film, coupled to a boundary value system for the electric potential Ψ(x, z, t) and the director field n = (sin θ, cos θ), characterized in our 2D setup by a single angle θ(x, z, t).
Numerical techniques were used to investigate the temporal evolution of an initially flat film, with initial dimensionless height representative of thick (h 0 = 1.5), intermediate (h 0 = 1.0) or thin (h 0 = 0.5 or h 0 = 0.3) films. For thick films, the surface height evolves to an undulating steady-state profile in all cases considered: fluid collects in the regions between electrodes where the electric potential is near zero and electric field gradients are small, and thins at the electrode midpoints, where the electric and director field gradients are large. For thin and intermediate thickness films, fluid accumulates in regions of high electric potential and director field gradients, migrating from regions of small gradients. Unlike the thicker films, thin and intermediate films thin significantly in these high-gradient regions; no steady state is found and our simulations are halted when accuracy is lost in the solution for the director field. For the thinnest films, complex wrinkling patterns form, with the film thickening at the edges of the electrodes and further thinning observed at the electrode midpoints where the electric potential gradients are large, leading to global film minima at these points. Between the electrodes, film thinning is suppressed and local minima form above x = 1 and x = 3.
Additionally, we investigated the effect of different anchoring conditions on NLC film evolution; specifically we considered the effect of free surface homeotropic anchoring strength A on the surface height evolution by considering three different values: A = 1, A = 10 and A = 50. Changing A primarily affects the director configuration and its gradients, which in turn affect the surface height evolution. We observed that as A is increased the maxima of the film surface profile become more diffuse, possibly due to increased competition of different forces within the film. For sufficiently weak anchoring, in the thinnest film simulated (h 0 = 0.5) we also observe regions in the film where the director field is uniform and planar (Fig. 9), reflecting the fact that θ = π/2 is a stable steady solution to the governing equations under the local film conditions. Finally, with the goal of gaining insight into the behavior of film evolution when there are no elastic forces within the film layer, we considered an NLC film with negative dielectric anisotropy, subject to strong planar anchoring at the lower substrate and weak planar free surface anchoring. Negative dielectric anisotropy means that the nematic molecules align perpendicular to the direction of the electric field instead of parallel to it; the director in this case is uniform throughout the layer, and the elastic energy is negligible. While surface height evolution was generally similar to that of analogous films with planar-homeotropic anchoring, the films with planar anchoring and negative dielectric anisotropy had well-defined maxima, which we attribute to the lack of competition between opposing forces in the film. The thinnest films in this case exhibit complex non-sinusoidal wrinkling profiles, behavior that appears to depend primarily on the initial film height.
We note that our numerical results are in qualitative agreement with the experiments of Brown et al. [1] who investigate dielectrophoresis of an isotropic dielectric fluid over interlaced electrodes. They observe that fluid collects above electrodes, where electric potential gradients are highest, and is removed from regions between electrodes where electric potential gradients are small. A similar experiment carried out by the same group [3] demonstrated that as the initial height of the dielectric film decreases, wrinkled, strongly non-sinusoidal profiles emerge, similar to those in our thinnest film simulations. Though we anticipated that the NLC film with planar free surface anchoring and negative dielectric anisotropy might be more analogous to the isotropic dielectric liquid (having zero net elastic energy), this expectation was not borne out: our simulations of planarhomeotropic flat films provide the best agreement with the experimental results of Brown et al. [1]. We hope that our new model and computational results will inspire future experimental investigations that may provide further insight and lead to improvement of our mathematical model.

V. ACKNOWLEDGMENTS
This work was supported by the NSF under Grant No. DMS 1615719.