Mingfeng Qiu
email: mingfeng.qiu@ens.fr

James J Feng
email: james.feng@ubc.ca

Jean-Christophe Loudet
email: jean-christophe.loudet@u-bordeaux.fr

Phase-field model for elastocapillary flows of liquid crystals

We propose a phase-field model to study interfacial flows of nematic liquid crystals that couple the capillary forces on the interface with the elastic stresses in the nematic phase. The theoretical model has two key ingredients, a tensor order parameter that provides a consistent description of the molecular and distortional elasticity, and a phase-field formalism that accurately represents the interfacial tension and the nematic anchoring stress by approximating a sharp-interface limit. Using this model, we carry out finite-element simulations of drop retraction in a surrounding fluid, with either component being nematic. The results are summarized by eight representative steady-state solutions in planar and axisymmetric geometries, each featuring a distinct configuration for the drop and the defects. The dynamics is dominated by the competition between the interfacial tension and the distortional elasticity in the nematic phase, mediated by the anchoring condition on the drop surface. As consequences of this competition, the steady-state drop deformation and the clearance between the defects and the drop surface both depend linearly on the elastocapillary number.

I. INTRODUCTION

We consider two-phase flows of a nematic liquid crystal (LC) and an immiscible isotropic fluid under coupled elastic stresses of the LC and capillary forces of the interface. Understanding such flows is important for a number of fundamental and applied scenarios, such as drop dynamics with a nematic drop or host [e.g., [START_REF] Doane | Polymer dispersed liquid crystals for display application[END_REF][START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF][START_REF] Mushenheim | Straining soft colloids in aqueous nematic liquid crystals[END_REF][START_REF] Scheff | Tuning shape and internal structure of protein droplets via biopolymer filaments[END_REF], selfassembly of soft colloids in a nematic phase [e.g., [START_REF] Poulin | Inverted and multiple nematic emulsions[END_REF][START_REF] Muševič | Two-dimensional nematic colloidal crystals self-assembled by topological defects[END_REF], elastocapillary interaction of particles on LC interfaces [e.g., [START_REF] Smalyukh | Ordered droplet structures at the liquid crystal surface and elastic-capillary colloidal interactions[END_REF][START_REF] Oettel | Effective interactions of colloids on nematic films[END_REF][START_REF] Jeridi | Capillary-induced giant elastic dipoles in thin nematic films[END_REF], and dynamics of biological matter and cells [e.g., [START_REF] Rofouie | Morphology of elastic nematic liquid crystal membranes[END_REF][START_REF] Metselaar | Topology and morphology of self-deforming active shells[END_REF][START_REF] Leoni | Defect driven shapes in nematic droplets: analogies with cell division[END_REF]. The presence of an interface leads to the interplay among surface tension, the anchoring of LC molecules, and the bulk elasticity, and gives rise to novel dynamics. For instance, Rey developed a general continuum theory for nematic-isotropic interfaces [START_REF] Rey | Marangoni flow in liquid crystal interfaces[END_REF][START_REF] Rey | Viscoelastic theory for nematic interfaces[END_REF][START_REF] Rey | Irreversible thermodynamics of liquid crystal interfaces[END_REF][START_REF] Rey | Dynamical phenomena in liquid-crystalline materials[END_REF], and predicted a Marangoni flow along the surface gradient of the anchoring energy as well as a normal force on the interface from LC molecular alignment deviating from easy directions.

To model the behaviour of elastocapillary phenomena involving a nematic LC, there have been mainly three types of approaches. First, many existing studies prescribe various strong simplifications, such as a fixed interfacial shape [START_REF] Jeridi | Colloidal particles in thin nematic wetting films[END_REF], a class of candidate shapes [START_REF] Prinsen | Shape and directorfield transformation of tactoids[END_REF] or the absence of defects in the bulk of the LC [START_REF] Van Bijnen | Texture and shape of two-dimensional domains of nematic liquid crystals[END_REF]. It is clear, however, that these assumptions rely on decoupling elasticity and capillarity and apply only in limiting situations. Second, some researchers adopt a sharp-interface description of the moving front, and couple the bulk flow equations with the constitutive relations characterizing the nematic-isotropic interface. Care et al. [START_REF] Care | Generalized lattice boltzmann algorithm for the flow of a nematic liquid crystal with variable order parameter[END_REF] proposed a lattice Boltzmann method (LBM) based on a continuum theory of the interface similar to that of Rey's, and applied it to study the shape of a drop submerged in a nematic LC [START_REF] Lishchuk | Shape of an isotropic droplet in a nematic liquid crystal: The role of surfactant[END_REF][START_REF] Lishchuk | A lattice boltzmann scheme for a nematic-isotropic interface[END_REF]. However, because of the complexity of the formulation needed for a sharp nematic-isotropic interface coupled with dynamic flow, these models are difficult to derive and apply in time-dependent phenomena. In particular, the interfacial stress in [START_REF] Care | Generalized lattice boltzmann algorithm for the flow of a nematic liquid crystal with variable order parameter[END_REF] lacks a viscous component, thereby only allowing stationary cases to be considered.

The third type of approach uses a diffuse-interface description based on phase-field (PF) models, the energetic formulations of which constitute a natural advantage in modeling complex fluids such as the LC whose relaxation is driven by the minimization of a free energy. The diffuse interface makes it straightforward to incorporate interfacial dynamics such that there is no need for sophisticated interfacial equations. Yue et al. [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF] proposed a PF method coupled with a regularized Leslie-Ericksen (LE) theory. This model has been used to examine the drop retraction and deformation with either the drop or host being nematic [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF], the dynamics of a Newtonian drop rising in a nematic LC [START_REF] Zhou | The rise of newtonian drops in a nematic liquid crystal[END_REF], and the self-assembly of drops in a nematic phase [START_REF] Zhou | Dynamic simulation of droplet interaction and self-assembly in a nematic liquid crystal[END_REF]. However, the regularization scheme for topological defects employed in those studies does not admit defects of half-integer winding numbers in the bulk of LC [START_REF] Bethuel | Ginzburg-Landau Vortices[END_REF]. As consequences, the ring defects were predicted to be always attached on the interface, and drop interaction was mediated by defects with integer winding numbers in 2D, which were known to be energetically unfavourable. Metselaar et al. [START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF] formulated a PF model to probe the shape of tactoids in a lyotropic LC, and recently generalized the model to investigate the dynamics of microtubules cross-linked and aligned by active molecular motors residing in the mixing region of two fluids-the so-called "active nematic shell" [START_REF] Metselaar | Topology and morphology of self-deforming active shells[END_REF]. Nevertheless, the details of the surface tension model in [START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF] are not explained and it is unclear how the model approaches a sharp-interface limit and represents the interfacial forces.

Generally, to account for fully-coupled interfacial flows of LC, one needs to overcome two challenges: (i) to describe the LC microstructure with potential topological defects accurately, and (ii) to describe the moving interface and related interfacial boundary conditions accurately. From the above survey, existing models typically have employed components which make it challenging to satisfy one of the two requirements. Hence, it is the objective of the present paper to construct a model able to meet both requirements simultaneously.

First, to describe the nematic order and its hydrodynamics, we use a tensor order parameter similar to [START_REF] Care | Generalized lattice boltzmann algorithm for the flow of a nematic liquid crystal with variable order parameter[END_REF][START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF] that stays smooth regardless of topological defects and keeps the elastic energy density bounded. Thus, artificial regularization is no longer needed. Second, to represent the moving interface, we choose the PF method used in [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF][START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF][START_REF] Zhou | The rise of newtonian drops in a nematic liquid crystal[END_REF][START_REF] Zhou | Dynamic simulation of droplet interaction and self-assembly in a nematic liquid crystal[END_REF][START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF]. This model approximates a sharp-interface limit and can thus recover macroscopic surface tension accurately. In the present context, this property also makes it possible to capture the nematic anchoring stress. Our model is numerically solved with a finite element method, making it easier to accommodate complex geometries compared with earlier LBM-based methods. The model is presented in Sec. II. To demonstrate the capability of our method, we study the dynamics of drop retraction involving a nematic phase in Sec. III. We provide a comprehensive map of drop and defect configurations with various combinations of materials, anchoring conditions, and spatial dimensions. Our results show that the competition between the bulk elastic distortion of the LC and surface tension dictates the process and consequence of drop retraction.

II. THEORETICAL MODEL AND NUMERICAL METHOD

A. Nematic order and free energies

We use the tensor order parameter Q to describe the nematic order [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]. It is defined from the secondorder moment of the distribution function for individual nematic molecule orientations u: Q = uu -I/3 , where I is the unit tensor and the operator • denotes average over an ensemble of molecules. Following literature [e.g., [START_REF] Larson | The Structure and Rheology of Complex Fluids[END_REF][START_REF] Feng | Orientational defects near colloidal particles in a nematic liquid crystal[END_REF], we define a scalar order parameter:

q = 3/2 Q F = 3 2 Q : Q 1/2
, where • F is the Frobenius norm. When the nematic phase is in a uniaxial state, we recover, as expected, the relationship Q = q (nn -I/3), where n is the nematic director. Note that some researchers adopt a different definition by requiring q to be proportional to the largest-in-magnitude eigenvalue of Q [e.g., [START_REF] Andrienko | Computer simulation of topological defects around a colloidal particle or droplet dispersed in a nematic host[END_REF][START_REF] Guzmán | Defect structure around two colloids in a liquid crystal[END_REF][START_REF] Hung | Anisotropic nanoparticles immersed in a nematic liquid crystal: defect structures and potentials of mean force[END_REF].

We introduce a phase variable φ to mark the different fluids. φ = -1 denotes the nematic LC and φ = 1 denotes the isotropic phase. These two fluids mix in the thin diffuse interfacial region in which φ and other variables transition smoothly. The total free energy of the fluids includes the mixing energy of the interface, the bulk elastic energy of the nematic phase, and the anchoring energy on the interface. The mixing energy and the bulk elastic energy are standard from the PF model [e.g., 24] and the nematic order theory [e.g., 37], respectively. The mixing energy is described by the Landau-Ginzburg density

f m = λ 2 |∇φ| 2 + λ 4 2 φ 2 -1 2 , ( 1
)
where λ is the mixing energy strength with the dimension of force, and is the constant capillary width governing the thickness of the diffuse interface. The bulk elastic energy in the nematic phase can be described phenomenologically as [START_REF] Ravnik | Landau-de gennes modelling of nematic liquid crystal colloids[END_REF]

f b = f LdG + f e = A 2 Q ij Q ij + B 3 Q ij Q jk Q ki + C 4 (Q ij Q ij ) 2 + L 1 2 (∂ i Q jk )(∂ i Q jk ). (2)
Here A, B, C are material property coefficients (A < 0, B < 0), and L 1 is the bulk elastic constant. The first three terms constitute the celebrated Landau-de-Gennes (LdG) free energy expansion f LdG for local molecular relaxation (microelasticity). Since we are only concerned with thermotropic LCs at a fixed temperature in this paper, A, B, C are constants. Minimizing the total LdG energy F LdG = f LdG dV under the uniaxial assumption gives the equilibrium scalar order parameter q e . We consider a nematic phase far away from phase transitions with strong molecular interaction. In this case,

q e = 3 4   - B 3C + B 3C 2 - 8A 3C   . (3)
The last term in Eq. [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF],

f e = (L 1 /2)(∂ i Q jk )(∂ i Q jk )
, describes the long-range orientational distortions (distortional elasticity). For simplicity, we have applied the common one-elastic-constant approximation [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF] and further omitted all mixed derivatives of the Q components.

The anchoring energy enforces a finite-strength anchoring condition on the interface by penalizing LC molecules deviating from the easy direction. We propose the following diffuse-interface form of the Rapini-Papoular an-choring energy density [START_REF] Rapini | Distorsion d'une lamelle nématique sous champ magnétique conditions d'ancrage aux parois[END_REF][START_REF] Nobili | Disorientation-induced disordering at a nematic-liquid-crystal-solid interface[END_REF]

f a = W 2 |∇φ| 4 Q -q e ee - 1 3 I 2 F , (4)
where W is the constant anchoring strength, and e is a unit vector along the easy direction on the interface. More will be discussed on the value of W in Sec. II C. The factor |∇φ| 4 avoids singularities in the expression of ee (see Eq. ( 5) and Eq. ( 6)). Note that one may consider alternative forms of anchoring energy, such as a LdG style energy expansion truncated at low orders of |∇φ|, which warrant future investigations. We consider two common choices for the easy direction:

(i) Homeotropic anchoring. The easy direction is the same as the unit normal vector m = ∇φ/|∇φ| of the interface. Hence

f h a = W 2 |∇φ| 2 Q -q e ∇φ∇φ - 1 3 |∇φ| 2 I 2 F . (5)
(ii) Planar anchoring in 2D planar and axisymmetric geometries. The easy direction is along a tangential direction to the interface. In 2D planar geometry, e can be written as e = (m 2 , -m 1 , 0), where m 1 and m 2 are components of m. It follows that

f p a = W 2 |∇φ| 2 Q + q e ∇φ∇φ - 1 3 |∇φ| 2 J 2 F , (6)
where the tensor J can be written as

J =   2 0 0 0 2 0 0 0 -1   , (7)
which is rotationally symmetric in the 1-2 coordinate plane. Therefore, as expected, the above expression of J does not rely on a specific coordinate frame. In the 2D axisymmetric geometry, we only consider the case where e stays in the r-z plane in a cylindrical frame, i.e., a "monostable" planar anchoring [START_REF] Jerome | Surface effects and anchoring in liquid crystals[END_REF]. Then, one needs to swap the second and third components in the diagonal of Eq. ( 7).

Hence, the total free energy of the two-phase mixture in a domain Ω is

F = Ω f (φ, Q, ∇φ, ∇Q) dV = Ω f m + 1 -φ 2 f b + f a dV. (8)
The factor (1 -φ)/2 gives the concentration of the nematic phase.

B. Governing equations

The governing equations include an evolution equation for the phase field φ (Cahn-Hilliard), an evolution equation for the tensor Q, and the equations of motion for the fluid flow.

The Cahn-Hilliard equation. The classical Cahn-Hilliard (CH) equation describes the evolution of φ [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF] ∂φ ∂t

+ v • ∇φ = γ∇ 2 µ, (9)
where v is the velocity vector, and the constant γ is the mobility parameter (unit: m 2 /(Pa•s)). The chemical potential µ is

µ = δF m δφ = λ -∇ 2 φ + φ(φ 2 -1) 2 , (10)
where F m = Ω f m dV is the total mixing energy. Thus µ only includes the relaxation of the mixing energy, without the contribution from the bulk elastic and anchoring energies. The same simplification has been adopted in previous works [e.g., [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF][START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF]. It allows for convenient recovery of interfacial forces, which will be elaborated in Sec. II C. The evolution equation for nematic order. There have been a number of tensor-based theories for nematic hydrodynamics [e.g., [START_REF] Olmsted | Theory of the nonequilibrium phase transition for nematic liquid crystals under shear flow[END_REF][START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals Part 1. Model formulation using the Hamiltonian (poisson bracket) formulation[END_REF][START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals: Part 2. Non-homogeneous systems[END_REF][START_REF] Beris | Thermodynamics of Flowing Systems with Internal Microstructure[END_REF][START_REF] Qian | Generalized hydrodynamic equations for nematic liquid crystals[END_REF], including kinetic theories that do not start from the same elastic energy f b as in this paper [e.g. , 46]. Each gives a slightly different evolution equation for Q. Yet, all of them predict qualitatively similar results [e.g., [START_REF] Feng | Orientational defects near colloidal particles in a nematic liquid crystal[END_REF][START_REF] Ravnik | Landau-de gennes modelling of nematic liquid crystal colloids[END_REF][START_REF] Tóth | Hydrodynamics of topological defects in nematic liquid crystals[END_REF][START_REF] Svenšek | Hydrodynamics of pairannihilating disclination lines in nematic liquid crystals[END_REF]. So far, experimental evidence has been inadequate to discriminate among these models. Therefore, we choose the Beris-Edwards (BE) theory [START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals Part 1. Model formulation using the Hamiltonian (poisson bracket) formulation[END_REF][START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals: Part 2. Non-homogeneous systems[END_REF][START_REF] Denniston | Lattice Boltzmann simulations of liquid crystal hydrodynamics[END_REF] based on its wide usage. In our PF formulation, the evolution equation for Q has the same general expression as in the BE theory in the bulk of LC [START_REF] Denniston | Lattice Boltzmann simulations of liquid crystal hydrodynamics[END_REF] ∂Q ∂t

+ v • ∇Q = S + ΓH, ( 11
)
where Γ is the collective rotational diffusion coefficient and is assumed to be a constant in this paper. The molecular field H is still given by H = -δF/δQ + tr (δF/δQ) I/3. Considering Eq. ( 8), we write H as

H = 1 -φ 2 M -G + ∇ • 1 -φ 2 L 1 ∇Q , (12)
where the first and third terms are the same as in the BE theory except for the LC concentration as an additional coefficient. M is a polynomial in Q that stems from the LdG energy:

M(Q) = -AQ-BQ•Q-(B/3)(Q : Q )I- C(Q : Q )Q.
The tensor G is due to the anchoring energy. In cases of homeotropic anchoring (Eq. ( 5)) and planar anchoring (Eq. ( 6)), respectively,

G h = W |∇φ| 2 |∇φ| 2 Q -q e ∇φ∇φ - 1 3 |∇φ| 2 I , (13)
G p = W |∇φ| 2 |∇φ| 2 Q + q e ∇φ∇φ - 1 3 |∇φ| 2 J . ( 14
)
The PF version of the corotation tensor S is defined as

S(∇v, Q) = (1 -φ)ξ 2 D + Ω • Q + 1 3 I + Q + 1 3 I • (1 -φ)ξ 2 D -Ω -(1 -φ)ξ (Q : (∇v) ) Q + 1 3 I , (15)
where D = ((∇v) + ∇v) /2 is the rate of deformation tensor, and Ω = ((∇v) -∇v) /2 is the vorticity tensor. ξ is a material parameter that determines, together with A, B, C, if the LC material is a flow-aligning or tumbling nematic [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Denniston | Lattice Boltzmann simulations of liquid crystal hydrodynamics[END_REF]. Note that in arriving at Eq. ( 15) from the original S in the BE theory, we have postulated the phenomenology of interpolating ξ with the nematic concentration (1 -φ)/2. It ensures that this equation vanishes in the isotropic phase when φ = 1. It avoids the complexity of designing Poisson brackets, similar to what was done in the original BE theory [START_REF] Beris | Thermodynamics of Flowing Systems with Internal Microstructure[END_REF], with the addition of the phase field. Provided that Q vanishes in the isotropic fluid and the interface is thin, this phenomenology does not alter the dynamics of Q in the nematic phase. The modified H and S in the PF formulation preserve the symmetry of the tensors in the original BE theory, as expected.

The equations of motion for the fluids and the LC stress tensor. The continuity and Navier-Stokes (NS) equations for an incompressible fluid govern the velocity and pressure fields:

∇ • v = 0, ρ ∂v ∂t + ρv • ∇v = ρg + ∇ • τ + µ∇φ. ( 16
)
Here ρ is the interpolated density of the two-fluid mixture

ρ = 1 + φ 2 ρ i + 1 -φ 2 ρ n , (17)
where ρ i and ρ n are the densities of the isotropic and nematic phases. g is the gravitational acceleration. The last term is a body-force form of surface tension from the PF formulation [START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF]. The stress tensor τ can be written in the following form:

τ = -pI + (1 + φ)η i D + (1 -φ)η n D + τ n , (18)
where p is pressure, and the second term is the viscous stress from the isotropic phase. The third term is a viscous stress of the nematic phase with a constant viscosity η n , independent of the molecular configuration Q [START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals Part 1. Model formulation using the Hamiltonian (poisson bracket) formulation[END_REF][START_REF] Edwards | Generalized constitutive equation for polymeric liquid crystals: Part 2. Non-homogeneous systems[END_REF][START_REF] Beris | Thermodynamics of Flowing Systems with Internal Microstructure[END_REF].

The last term is the nematic stress tensor τ n , which will be modified from the BE theory to take into account the anchoring condition on the moving interface. Instead of re-deriving this stress tensor using Poisson brackets, we take a simple but equivalent shortcut. Following [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF], we decompose the effect of anchoring into two parts, each of which can be determined from a distinct virtual operation. (i) From a virtual rotation and deformation of the molecular distribution without motion of the fluid. This part is already accounted for in the molecular field (Eq. ( 12)). (ii) From a virtual strain of the fluid with the Q configuration frozen in each fluid particle. This variational procedure gives rise to an additional stress tensor associated with ∇φ, termed the anchoring Ericksen stress τ ae . τ ae can be derived by applying the virtual work principle to the total free energy (Eq. ( 8)) [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF]. For homeotropic (Eq. ( 5)) and planar (Eq. ( 6)) anchoring, respectively, we have

τ h ae ij = - ∂f h a ∂(∂ i φ) (∂ j φ) = -W |∇φ| 2 Q kl -q e (∂ k φ)(∂ l φ) - 1 3 |∇φ| 2 δ kl • 2Q kl (∂ i φ) -q e δ ik (∂ l φ) + (∂ k φ)δ li - 2 3 (∂ i φ)δ kl (∂ j φ), (19)
(τ p ae ) ij = - ∂f p a ∂(∂ i φ) (∂ j φ) = -W |∇φ| 2 Q kl + q e (∂ k φ)(∂ l φ) - 1 3 |∇φ| 2 J kl • 2Q kl (∂ i φ) + q e J ik (∂ l φ) + (∂ k φ)J li - 2 3 (∂ i φ)J kl (∂ j φ). (20)
As expected, τ ae is not symmetric in general. Overall, the PF version of the nematic stress tensor is

(τ n ) ij = - 1 -φ 2 ξH ik Q kj + 1 3 δ kj + ξ Q ik + 1 3 δ ik H kj -2ξ Q ij + 1 3 δ ij (Q kl H kl ) + L 1 (∂ i Q kl ) (∂ j Q kl ) + H ik Q kj -Q ik H kj + (τ ae ) ij . ( 21
)
Boundary conditions. For the evolution of the order parameter Q, there are four types of BCs. (i) On boundaries Γ QN where the gradient of Q vanishes, the homogeneous Neumann BC m • ∇Q| Γ QN = 0 applies. (ii) On boundaries Γ QD with infinite-strength anchoring, we use the Dirichlet BC Q| Γ QD = Q D where Q D is the nematic configuration prescribed by the anchoring condition. (iii) On boundaries Γ QR with finite-strength anchoring, we assume a Rapini-Papoular anchoring energy and equilibrium on Γ QR . This leads to the Robin con-

dition ( W (Q -Q) + L 1 m • ∇Q) Γ QR = 0
, where W and Q are the anchoring strength (surface energy density) and the easy nematic configuration on Γ QR , respectively. (iv) On a symmetry boundary Γ S , one of the three symmetry planes of the ellipsoid described by Q + I/3 needs to be tangential to Γ S . Assuming that Γ S is parallel to a coordinate plane, we can thus impose the mixed BC:

m • ∇ (Q ii )| Γ S = 0, Q ij | Γ S = 0 (i = j)
, where ii does not imply the Einstein summation.

For the continuity and NS equations for the fluid flow, standard BCs can be used. For example, on boundaries that mimic infinity or on Γ S , we can use no penetration and free slip conditions. The fourth-order CH equation needs two BCs. First, one often requires that there be no diffusive flux across all boundaries: m • ∇µ| ∂Ω = 0. The second BC typically deals with φ or the gradient of φ. For instance, on a boundary in contact with a single phase or with a pinned contact line, we impose fixed values of the phase variable φ. The current paper does not consider dynamic contact lines or contact angles, which require different types of BCs and will be the subject for future investigation.

C. Interfacial forces

Since our computations involve a nematic-isotropic interface, the interfacial forces include both surface tension and the anchoring stress. These quantities are clearly defined in a sharp-interface model, but need to be transformed into our diffuse-interface framework.

For this purpose, we consider a flat diffuse interface at local equilibrium, i.e. µ = 0. Let x be the local coordinate normal to the interface, with x = 0 at φ = 0. Using Eq. ( 10) one obtains

φ e (x) = tanh x √ 2 . ( 22
)
This equation also holds asymptotically ( → 0) for curved interfaces in general. Requiring that the mixing energy of the diffuse interface be equal to the surface energy of a sharp interface, one arrives at the matching condition [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF]

σ = 2 √ 2 3 λ , ( 23
)
where σ is the surface tension coefficient. The same local interfacial profile also allows us to specify the PF expression of the anchoring strength W (Eq. ( 4)). Consider a sharp interface S W that is represented in our PF model by an interfacial region Ω W . In the sharp-interface formulation, the Rapini-Papoular anchoring energy is written as the following surface energy density

f s = W s 2 Q -q e ee - 1 3 I 2 F , (24)
where W s is the anchoring strength with the dimension of force/length. Equating the anchoring energies in the sharp-interface and diffuse-interface models yields,

S W W s 2 Q -q e ee - I 3 F dS = Ω W W 2 |∇φ| 4 Q -q e ee - I 3 F dx. (25)
We are free to make S W and Ω W arbitrarily small, and can thus assume that Q -q e (ee -I/3) F is a constant. Using the equilibrium φ profile in Eq. ( 22), we obtain

W = W s ∞ -∞ (φ e ) 4 dx = 35 8 √ 2 W s 3 , (26)
where the prime denotes differentiation with respect to x, and W bears the expected dimension of energy×length (see Eq. ( 4)).

We include only the mixing energy in µ (Eq. ( 10)) in order to derive the local interfacial profile in Eq. ( 22), which in turn allows us to establish the correspondence of interfacial forces between the sharp and diffuse interface descriptions (Eq. ( 23) and Eq. ( 26)). The values of λ and W would have depended on Q otherwise, resulting in a more complicated surface tension model. Compared to the PF method used in [START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF], our model is a simplification in the sense that the former specifies a more complete chemical potential including all free energies.

With the simpler CH dynamics of Eq. ( 10), free from coupling with Q, we are able to exploit the wellestablished separation of time scales in the nonlinear CH diffusion [START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF]. As → 0, on the time scale of O( 4 /(γλ)), the interface stabilizes locally and develops the hyperbolic tangent profile (Eq. ( 22)). On the longer time scale of O( r 3 /(γλ)), where r is a characteristic length of interfacial shape such as the radius of curvature, CH diffusion drives an evolution of φ in addition to flow advection, known as Mullins-Sekerka flow in material science [START_REF] Mullins | Morphological stability of a particle growing by diffusion or heat flow[END_REF]. Therefore, we require to be small and γ to be judiciously chosen such that the flow time scale of interest falls between the two time scales of the CH dynamics mentioned above. In this way, the interface motion is governed by flow advection under the correct interfacial forcing. Practical guidelines on mitigating nonphysical long-time effects have been documented in [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF][START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF][START_REF] Yue | Spontaneous shrinkage of drops and mass conservation in phase-field simulations[END_REF].

D. Numerical approximation

We employ a Galerkin finite element method to solve the governing equations. We use a mixed formulation for the CH equation as is done in [START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF][START_REF] Liu | Mixed methods for the approximation of liquid crystal flows[END_REF], as well as for the NS equation by keeping τ n as an auxiliary variable. The weak form of the governing equations with applicable BCs is discretized on an unstructured mesh of triangles, using piecewise quadratic Lagrange elements (P2) for v, φ, µ, Q, τ n , and piecewise linear Lagrange elements (P1) for p. We employ non-uniform meshes to balance resolving defects and the thin interface and saving computational cost. (See Sec. III and Fig. 1 for an example.)

The model is implemented using the software package COMSOL Multiphysics ® [START_REF]COMSOL Multiphysics ® v. 5.4, www.comsol.com[END_REF] to make use of a se- lection of well-developed numerical solvers. For timestepping, we use a backward-difference-formula (BDF) scheme with adaptive time step sizes and adaptive orders of accuracy between first and second [START_REF] Brenan | Numerical Solution of Initial-value Problems in Differentialalgebraic Equations[END_REF]. LC hydrodynamics displays a range of physical time scales including relaxations of microelasticity, distortional elasticity and the flow field, in addition to the time scale at which the φ field establishes local equilibrium. Thus the adaptive time step size is crucial for the performance of our computation. In each time step, all equations are coupled and solved simultaneously. The nonlinear system is solved by a damped Newton's method. To solve each Newton iteration, we choose the multifrontal massively parallel sparse direct solver (MUMPS) [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF][START_REF]MUMPS: A parallel sparse direct solver[END_REF], an efficient direct linear solver. Since we only deal with 2D problems in this paper, a direct solver is advantageous in its robustness and accuracy without incurring excessive computational cost.

We note two practical considerations in the numerical scheme. First, on symmetry boundaries Γ S , τ n satisfies the same conditions as

Q: m • ∇ ((τ n ) ii )| Γ S = 0, (τ n ) ij Γ S
= 0 (i = j) (also see Sec. II B). This condition needs to be taken into account when choosing the approximation spaces for τ n in the weak form. Second, due to the relaxation of the total mixing energy in a finite domain, after the interfacial profile is stabilized, φ will not be exactly 1 or -1 in the bulk of the two fluids [START_REF] Yue | Spontaneous shrinkage of drops and mass conservation in phase-field simulations[END_REF]. In the current calculation, anti-diffusion could arise from the term ∇ • 1-φ 2 L 1 ∇Q (Eq. ( 12)) when φ shifts slightly above 1, producing spurious oscillations in the Q field. To mitigate this issue, we introduce a thresholding scheme. Let h s (x) denote a regularized Heaviside function that transitions smoothly from 0 to 1 within x ∈ [-d, d]. The advection term in Eq. ( 11) and the concentration factor (1 -φ)/2 in all the equations are multiplied by h s (-φ + φ c ). This ensures that the LC dynamics vanishes where φ > φ c + d. We choose the cutoff φ c = 0.9 and d = 0.05. Provided that we are close to the sharp-interface limit, specific φ c values have little impact on the flow field. We have tested smaller φ c values and the overall dynamics is hardly changed.

To validate our numerical model, we have performed extensive tests on single-phase nematic flows, including formation of Saturn ring and polar ring defects near a spherical solid particle, a flow-aligning nematic in simple shear flow in both infinite space or bounded by walls, and annihilation of defects of opposite winding numbers with and without flow coupling. We have also simulated a deformed Newtonian drop retracting in a Newtonian medium to validate our PF calculation. Quantitative agreement with literature has been achieved. The details are omitted here, but can be found in [START_REF] Qiu | Computational studies on interfacial dynamics in complex fluids[END_REF].

III. DROP RETRACTION WITH A NEMATIC-ISOTROPIC INTERFACE

We consider the problem of a freely suspended microsized liquid drop retracting in a quiescent medium, with one of the fluids being a nematic LC. Besides being relevant for a number of applications [START_REF] Doane | Polymer dispersed liquid crystals for display application[END_REF][START_REF] Mushenheim | Straining soft colloids in aqueous nematic liquid crystals[END_REF][START_REF] Poulin | Novel colloidal interactions in anisotropic fluids[END_REF], drop retraction provides a relatively simple ground for understanding the interplay among surface tension, the bulk elastic energy and the anchoring energy, under the constraint of geometry and topology. The notion of a dynamic surface tension, dependent on transient factors such as interfacial molecular alignment, has been demonstrated experimentally [START_REF] Lee | The deformation and retraction of thermotropic LCP droplets in a flexible polymer matrix[END_REF][START_REF] Wu | A study on interfacial tension between flexible polymer and liquid crystal[END_REF] and computationally [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF] in drop retraction with a nematic-isotropic interface. The results in [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF], however, suffered from limitations due to the regularization scheme in the modified LE model, as aforesaid (Sec. I). Previous modeling work also only covered a few scattered cases. Different combinations of drop-matrix materials, anchoring conditions and spatial dimensions remain to be explored. Here, using our new model, we perform a systematic numerical study of the drop retraction problem with a nematic-isotropic interface. On one hand, it validates our model by reproducing known features. On the other hand, we obtain new insight on elastocapillary flows.

The problem is set up in either 2D planar or axisymmetric geometries, as shown in Fig. 1(a). In all cases considered here, the geometry and nematic alignment have quadrupolar symmetry, so we only need to simulate one quarter of the domain. On the outer boundaries, we use non-penetrating, free-slip walls and specify the phase to be nematic or isotropic. For all cases where the nematic phase is outside the drop, we require the far-field LC molecules to be oriented in the vertical direction (ydirection in 2D planar geometry, and z-direction in 2D axisymmetric geometry) at equilibrium order (Eq. ( 3)). Initially the drop is always elongated in the vertical direction. In most cases, whether the nematic is outside or inside the drop, the LC molecules are initially oriented along the vertical direction at equilibrium order. Exceptions will be explained in specific cases.

We prescribe a sub-domain with highly refined, nearuniform mesh, covering the area through which the interface and topological defects move. Outside this subdomain, the mesh gradually coarsens. Fig. 1(b) shows the mesh structure in a typical simulation. The interfacial region and the defect core are both well resolved during the entire simulation.

Unless explicitly stated otherwise, the results have been obtained with the base parameters listed in the Appendix. The nematic and isotropic phases have a matched density and we set η i = η n as we are focused on the interaction of elasticity and capillarity. In presenting the results, we use dimensionless variables marked by an ϕ = -0.9 ~ 0.9 In the inset, the solid black lines show levelsets of φ = 0.9, φ = 0 (nominal interface) and φ = -0.9, and the shaded area (bottom left) represents the defect core (q ≤ 0.5).

asterisk. We scale length by the nominal drop radius R, time by η i R/σ, energy by L 1 R and energy (volume) density by L 1 /R 2 . Among the numerous dimensionless groups, the following three dimensionless lengths are important in governing the defect configuration and drop shape:

l * c = 1 R L 1 -A , l * e = L 1 W s R , l * ec = L 1 σR . (27)
The nematic coherence length l * c reflects the competition between distortional elasticity and microelasticity, and governs the defect core size. The anchoring extrapolation length l * e describes the competition between the bulk distortion of molecular order and the anchoring of LC molecules on the fluid-fluid interface. The elastocapillary length l * ec denotes the competition between the bulk elastic distortion and surface tension. Note that with our base parameters l * e = 0.01 1. Hence, we are in a regime where the anchoring is strong and topological defects are expected.

In the following sections, we first provide an overview of the different cases considered in this study, including the defect structures and drop deformation, in comparison with previous literature. Next, we explain how the drop deformation originates from the competition between bulk elasticity and surface tension. Lastly, we show the retraction speed of the drops and the interfacial Marangoni flow induced by anchoring energy gradients, in agreement with existing results.

A. Overview of different cases

With isotropic fluids, the drop retraction is solely driven by surface tension, or equivalently, by the reduction of the mixing energy. In the presence of an LC, however, depending on the initial bulk orientation of LC molecules and the prescribed anchoring conditions at the interface, various topological defects may appear and influence the retraction process. Distortions in the nematic alignment during retraction are likely to increase the elastic energies.

We first consider an illustrative case exemplifying the retraction process and energy tradeoff. Fig. 2 shows snapshots of a nematic drop retracting in an isotropic fluid with planar anchoring in the 2D planar geometry. The grey scale indicates the concentration-weighted scalar order parameter 1-φ 2 q. With our parameters, the equilibrium order is at q e = 0.81 (Eq. ( 3)). Consistent with the initial alignment, Fig. 2(a) shows a uniform molecular configuration without any defect at the beginning. While the drop retracts, the LC molecules near the interface reorient to reduce anchoring energy, and a defect nucleates near the pole (Fig. 2b,c), at the expense of both molecular and distortional elastic energies. Later, due to strong anchoring, the defect detaches from the interface and moves inside the drop, transitioning from a surface defect to a point defect of winding number +1/2 in the bulk (Fig. 2c,d ). Eventually the drop settles into a steady state and attains an elongated shape in the ydirection.

We summarize the steady-state defect and drop configurations observed in our study into eight different cases (Fig. 3). Case 1 corresponds to the scenario described above. Cases 2 -4 are drop retractions in 2D planar geometry as well, with different material combinations and anchoring conditions. For convenience of comparison, in Case 3, we set the initial molecular alignment to be along the x-direction, which leads to a final drop elongation in the y-direction. The same cases are calculated for 3D drops with axisymmetry in Cases 5 -8. Note that for Cases 1 -4 and 8, the defects are detached from the interface (as shown in the diagrams) when the surface tension or anchoring strength is high, otherwise it would be energetically favourable for the defects to stay bound to the interface. See Sec. III B for more details.

To our knowledge, the defect structures computed in Cases 1 -2 have not been reported before. Even though Yue et al. [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF] and Liu et al. [START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF] considered the same setups as in Cases 1 -3 here, they obtained different configurations because the formation of defects with half-integer winding numbers in the bulk of LC, like those in Cases 1 -3, were not permitted in their model (see Sec. I). In particular, for Cases 1 -2, the defects could detach from the interface in our model, while those in [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF] bound to the drop surface. In Case 3, [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF] predicted a point defect at the centre of a circular drop, while our calculation shows a pair of defects of half-integer winding number, deforming the drop. Notice the similarity in the defect structure in our Cases 3 -4 with those in Fig. 4 of Mackay and Denniston [START_REF] Mackay | Deformable vesicles interacting in a nematic liquid crystal[END_REF], who computed the configuration of an elastic vesicle shell suspended in a continuous phase of nematic LC in 2D. In [START_REF] Mackay | Deformable vesicles interacting in a nematic liquid crystal[END_REF], the homeotropic anchoring induced defects both inside and outside the shell, which have been reproduced separately in our Cases 3 -4. It also reflected well image charge arguments in the strong anchoring limit. For Case 4, our computed defect agrees with earlier studies on 2D particles or drops submerged in the bulk of a nematic LC [START_REF] Care | Generalized lattice boltzmann algorithm for the flow of a nematic liquid crystal with variable order parameter[END_REF][START_REF] Fukuda | Director configuration and dynamics of a nematic liquid crystal around a twodimensional spherical particle: Numerical analysis using adaptive grids[END_REF].

In the axisymmetric geometry, Cases 5 -6 feature surface-bound defects, the so-called "boojums", a consequence of the topological constraint of the spherical surface in 3D [START_REF] Kamien | The geometry of soft materials: a primer[END_REF][START_REF] Lavrentovich | Topological defects in dispersed words and worlds around liquid crystals, or liquid crystal drops[END_REF][START_REF] Serra | Curvature and defects in nematic liquid crystals[END_REF]. This is in contrast with the detached defects in Cases 1 -2. In a 3D drop with homeotropic anchoring in Case 7, a central point defect opens up spontaneously into a small ring. This is a known feature of the LdG free energy with the one-elastic-constant approximation [START_REF] Ravnik | Landau-de gennes modelling of nematic liquid crystal colloids[END_REF][START_REF] Fukuda | Defect structure of a nematic liquid crystal around a spherical particle: adaptive mesh refinement approach[END_REF]. In Case 8 with an isotropic drop in a nematic host with homeotropic anchoring, we have found a Saturn ring in the bulk of LC or a surface ring attached FIG. 3. Illustrated final defect and drop configurations after drop retraction. Defects are marked by green points or curves. The thin solid lines illustrate streamlines of the director n. The drop deformation may be exaggerated for clarity. P = planar anchoring. H = homeotropic anchoring. These anchoring conditions are illustrated in the boxes for Cases 1 and 3. For all cases with a nematic host [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Mushenheim | Straining soft colloids in aqueous nematic liquid crystals[END_REF][START_REF] Poulin | Inverted and multiple nematic emulsions[END_REF][START_REF] Smalyukh | Ordered droplet structures at the liquid crystal surface and elastic-capillary colloidal interactions[END_REF], in the far field the molecular alignment is along the vertical (y-or z-) direction at equilibrium order. For Cases 1 -4 and 8, the defects may be attached on the interface or not depending on parameters.

on the drop, depending on the parameters. Consistent with experimental observations [START_REF] Khullar | Dynamic evolution of topological defects around drops and bubbles rising in a nematic liquid crystal[END_REF], our result is more physically realistic than that in [START_REF] Zhou | The rise of newtonian drops in a nematic liquid crystal[END_REF], where the ring defect was confined to the drop surface. Again it is because of the regularized LE theory used in the latter study. In addition, in Case 8, it is well known that a dipolar configuration, which consists of a companion hyperbolic point defect located in the close vicinity of the particle, may arise as well [START_REF] Ravnik | Landau-de gennes modelling of nematic liquid crystal colloids[END_REF]. Actually, because of our assumption of isotropic elasticity in the LdG free energy, the point defect will spontaneously open up into a small ring, as in Case 7. However, with our material parameters A, B, C, L 1 (see the Appendix), the dipolar pattern is not stable with respect to the Saturn ring [START_REF] Fukuda | Nematic liquid crystal around a spherical particle: Investigation of the defect structure and its stability using adaptive mesh refinement[END_REF][START_REF] Wang | Topological defects in an unconfined nematic fluid induced by single and double spherical colloidal particles[END_REF].

The non-spherical shapes attained at the end of the retraction in 2D planar geometry (Cases 1 -4) agree with results reported in various earlier studies [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF][START_REF] Van Bijnen | Texture and shape of two-dimensional domains of nematic liquid crystals[END_REF][START_REF] Lishchuk | A lattice boltzmann scheme for a nematic-isotropic interface[END_REF][START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF]. The shapes of axisymmetric drops (Cases 5 -8), although reported less frequently, are not surprising considering their 2D planar counterparts. Here we suggest a heuristic argument that can rationalize intuitively the equilibrium drop shapes in all cases by considering the interaction between the director field (thin lines in Fig. 3) and the interface. In Case 3, for instance, we first consider a circular drop with a director field consistent with the homeotropic anchoring on the interface, illustrated in Fig. 4(a). Because of the anchoring, bend and splay distortions arise in the director field, which exert torques on the drop interface, as shown by the semi-circles with arrows in Fig. 4(a). These torques will tend to rotate the interface and elongate the drop in the y-direction, as shown in Fig. 4(b). This argument is consistent with the well-known fact that defects of equal winding numbers repel each other elastically. In Case 7, the defect is far from the interface so that the alignment is approximately radial at the interface, thereby producing a nearly spherical drop. In Case 8, by a similar argument to that in Case 3, nematic distortions near the interface tend to "flatten" the drop, resulting in an oblate shape. Our explanation agrees with the argument put forward in [START_REF] Mackay | Deformable vesicles interacting in a nematic liquid crystal[END_REF], although in the latter study, the deformation is additionally resisted by the bending rigidity of the vesicle shell. We have not attempted a more quantitative comparison on drop shapes with prior literature. Such a comparison is hampered by strong assumptions on drop shape classes or defect locations required in analytical solutions [e.g., [START_REF] Prinsen | Shape and directorfield transformation of tactoids[END_REF][START_REF] Van Bijnen | Texture and shape of two-dimensional domains of nematic liquid crystals[END_REF], or difficulties in mapping parameters among different models for nematic order and interfacial stresses used in computational studies [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF][START_REF] Lishchuk | A lattice boltzmann scheme for a nematic-isotropic interface[END_REF][START_REF] Metselaar | Electric-field-induced shape transition of nematic tactoids[END_REF].

B. Drop deformation and competition of free energies

Now we examine the steady-state drop configurations in more detail. The qualitative analysis in the last paragraph of the above section suggests that the reduction in the bulk distortional elastic energy at the expense of the mixing energy (surface tension) is the dominant factor determining the drop deformation at steady state. This motivates a scaling argument for the final drop deformation, which has not been reported before. Let a s , b s denote the steady-state lengths of a and b, the semi-axes of the drop (Fig. 1a). The competition between the distortional elastic energy and the mixing en-ergy gives rise to the (dimensional) elastocapillary length l ec = l * ec R = L 1 /σ. Dimensional analysis suggests that the drop deformation, measured by b s -a s , be proportional to l ec . This scaling relation may be developed in more detail in the limit of small deformations. Consider a circular drop in the 2D planar geometry (Cases 1 -4) and a fully relaxed Q-field with a pair of defects. In the nematic phase, the distortional elastic energy density f e ≈ L 1 /ζ 2 , where ζ is a characteristic length of longrange Q distortions. For all 2D planar cases (and also Cases 5 and 6), one can choose ζ to be the distance between the point defects. Typically, for all cases reported in Fig. 3, ζ is on the order of R. We examine the change of free energies when the drop elongates by δl R in one direction. Without loss of generality, we assume that b s > a s , i.e., b s = R + δl. To conserve the drop area, to leading order, a s = R -δl. The distortional elastic energy density becomes approximately L 1 /(ζ + 2δl) 2 as the defects move apart by about 2δl. We further assume that the characteristic area of elastic distortion is A e , which stays constant under small deformations. Thus the change in total distortional elastic energy δF e ≈ (-4A e L 1 /ζ 3 )δl. Meanwhile, the change in the mixing energy δF s ≈ (3πσ/R)δl 2 . Equating their magnitudes gives an equilibrium δl corresponding to the lowest total of the two energies: b s -a s = 2δl ≈ (8A e RL 1 )/(3πζ 3 σ). Scaling all lengths by R, we obtain

b * s -a * s = 2δl * ≈ 8A * e 3πζ * 3 L 1 σR = 8A * e 3πζ * 3 l * ec . (28)
In 3D with axisymmetry, similar arguments give b * s -a * s ≈ (24V * e /πζ * 3 )l * ec , where V * e is the dimensionless characteristic volume of long-range distortion, scaled by R 3 .

Numerical data confirms this proportional relation. In Fig. 5 we show the dimensionless drop elongation |b * s -a * s | as a function of the dimensionless elastocapillary length l * ec with constant l * e and l * c , for Cases 1 -4 in 2D planar geometry (Fig. 5a) and for Cases 5 -8 in 2D axisymmetric geometry (Fig. 5b). In all cases, the deformation |b * s -a * s | is indeed approximately proportional to l * ec . In Case 7 the relation is degenerate with negligible drop deformation regardless of l * ec . This proportionality persists to moderately large magnitudes of |b * s -a * s |, around 0.2 in many cases and over 0.3 in Case 6.

In certain scenarios, the trade-off between the bulk elastic energy and the mixing energy can be demonstrated directly from the evolution of free energies in time. In Case 8, for instance, we calculate different components of the free energy with our PF model, in dimensionless forms:

F * m = 1 L 1 R f m dV, F * LdG = 1 L 1 R 1 -φ 2 f LdG dV, F * e = 1 L 1 R 1 -φ 2 f e dV, F * a = 1 L 1 R f a dV, (29)
where the anchoring energy. Let F * be the total energy, i.e., F * = F * m + F * LdG + F * e + F * a . Fig. 6 shows the evolution of these energies of a simulation in Case 8. The initial transients are dominated by the retraction of the drop and establishment of the defects. As we are most interested in the later stage of drop deformation, we have chosen matching ranges for the various energies to facilitate their comparison in the later stage. The mixing energy decreases during the initial contraction in the axial (z-) direction, and rises again, although very moderately, in the subsequent expansion in the radial (r-) direction (Fig. 6a, also see Case 8 in Fig. 3). The decreasing distortional elastic energy (Fig. 6c) drives the decay in the total free energy (Fig. 6e), while being resisted to an appreciable extent only by the mixing energy. Their antagonism indeed dominates the changes in the micro-elastic (LdG) (Fig. 6b) and the anchoring (Fig. 6d ) energies.

F
In addition, as the drop shape is dictated by the competition between the bulk distortion and the mixing energy, the defect position relative to the interface also depends on l * ec whereas the anchoring strength is kept constant relative to distortional elasticity (constant l * e ). We use d * s to denote the dimensionless steady-state distance from the defect core to the nominal interface (φ = 0 contour) in the drop elongation direction. Fig. 7 shows The anchoring energy. (e) The total energy. We have used a smaller surface tension in this run (l * ec = 0.033) than the base value, so as to magnify the drop deformation for ease of visualization and analysis. Steady state is reached near t * ≈ 200. For clarity here we only show the range up to t * = 100.

d *

s as a function of l * ec for all 2D planar cases as examples. When surface tension decreases (increasing l * ec ), the defect gradually approaches the drop interface and d * s decreases roughly linearly, again consistent with suggestions by dimensional analysis. In this process, although the anchoring energy may increase due to the disordered defect core being closer to the interface and disturbing molecular order, the total energy decreases thanks to relaxed bulk distortions. In Cases 1 and 3, the dependence of d * s on l * ec is weak, and the data points overlap, similar to those in Fig. 5(a). This is because of the confinement of the defects inside the drop. The dependence is stronger in Cases 2 and 4, with the nematic outside the drop. Note that at high values of l * ec with the nematic outside, the defects become attached to the interface, as shown by the insets (Fig. 7).

In the limit of vanishing anchoring (l * e → ∞), the bulk distortion and the mixing energy will be decoupled as there will be no nematic distortions, defects or drop deformation at steady state. Thus the anchoring energy mediates the exchange between the bulk distortion and the mixing energy by providing a topological constraint on the directors in the interfacial region. As examples, we show the effect of anchoring strength in Fig. 8 for Cases 1 and 2, with l * ec = 0.05. There are three features that s between the defect and the interface as a function of l * ec for 2D planar cases. The insets show the distribution of 1-φ 2 q at steady state for two scenarios in Case 2. The grey scale is the same as in Fig. 2. The dark region outside the drop marks the defect. warrant attention. First, when the anchoring strength decreases (increasing l * e ), the defects shift closer to the interface (Fig. 8a), as expected intuitively and in agreement with literature [e.g., 71]. The defects then undergo an abrupt transition and become attached to the interface (Fig. 8a and the insets therein). This transition is also observed by Lishchuk et al. [START_REF] Lishchuk | A lattice boltzmann scheme for a nematic-isotropic interface[END_REF] in their LBM calculation of a static drop shape similar to our Case 2. Second, as the defect becomes attached, the drop elongation jumps to a higher value, reducing the bulk distortion but increasing the mixing energy (Fig. 8b). This is another manifestation that the drop deformation is mainly driven by the reduction of the distortional elastic energy. Third, as anchoring further weakens, the coupling between the bulk distortion and the mixing energy is gradually lost, and the drop deformation decays, in agreement with our reasoning about l * e → ∞. Note that the first two features have counterparts in Fig. 7 when the surface tension is lowered. In the latter case, though, the transition is much smoother.

To summarize the analysis in this section, we see that the trade-off between the bulk elastic distortion and the mixing energy determines the drop deformation. The anchoring energy plays the role of a mediator and moderates this coupling. s | as functions of l * e for Cases 1 and 2. The insets show the distribution of 1-φ 2 q at steady state for two scenarios in Case 2. The grey scale is the same as in Fig. 2. The dark region outside the drop signals the defect.

C. Retraction dynamics

In Fig. 9, we analyze the retraction speed by monitoring the time evolution of the drop deformation parameter D = (b-a)/(b+a) (see Fig. 1 for definition of a, b) for all the cases considered in Fig. 3, in semi-logarithmic scale. Additionally, we perform two simulations of a Newtonian drop retracting in a Newtonian medium in 2D planar and axisymmetric geometries, respectively, with the viscosity and density matched to the isotropic fluid in our nematic-isotropic cases (Appendix). Designated by "N/N", these results serve as baselines for comparison. As known from previous studies [START_REF] Guido | Measurement of interfacial tension by drop retraction analysis[END_REF], the N/N retraction features an exponential decay of D toward 0. Yet with a nematic-isotropic interface, the decay of D deviates gradually from exponential, particularly in the late stage. In general, compared with N/N, the acceleration or deceleration depends on whether retraction facilitates or hinders the reduction of the bulk distortion. The drops in Cases 1, 2, 3, 5 and 6 retract slower than in their corresponding Newtonian cases, because a more spherical drop causes more bending and splay of the directors, which resists further retraction (also see the last paragraph in Sec. III A). Note that the trend in Cases 1 and 2 agrees well with the findings in [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF]. In contrast, those in Cases 4 and 8 retract faster than in N/N as the bulk distortions keep on decreasing while the drop first retracts in the y-or zdirection and then expands in the x-or r-direction. The drop in Case 7 retracts slightly slower than in N/N. In this configuration, the nematic distortions do not play a strong role in the drop retraction. The retraction speed can be faster or slower than in the N/N case, depending on initial conditions and parameter values. Finally, the gradient in the anchoring energy drives a Marangoni flow along the interface from locations of low anchoring energy to those of high anchoring energy, as predicted by Rey [START_REF] Rey | Marangoni flow in liquid crystal interfaces[END_REF]. It has been captured in our sim-ulations, as demonstrated in Fig. 10

IV. CONCLUSION

This paper reports direct numerical simulations of elastocapillary flows with liquid crystals (LCs). We have developed a computational framework that integrates a tensorial description of molecular order in nematic LCs with a phase-field (PF) formalism that approximates a sharp-interface limit of two-phase flows. This method allows for a physically consistent description of topological defect structures in the LC, while capturing interfacial forces including surface tension and the LC anchoring stress.

We demonstrate the capability of our method by applying it to the drop retraction problem with a moving nematic-isotropic interface. We have provided a comprehensive map of drop-defect configurations in a variety of cases, including a nematic drop in an isotropic medium and the inverse, with planar and homeotropic anchoring conditions at the interface, in 2D planar and 2D axisymmetric geometries.

One of our key findings is the interplay between the bulk elastic distortions and surface tension in dictating the drop shape. In particular, the dimensionless steadystate drop deformation is proportional to the elastocapillary length l * ec = L 1 /(σR), which controls the relative importance of bulk distortional elasticity and surface tension. This coupling is mediated by the anchoring energy, and gradually diminishes as anchoring becomes weaker. Our examination of this model problem highlights the interplay between the elastic energies of the LC and surface tension.

In addition, we have reproduced many features observed in prior studies, including defect configurations, elongated drop shapes, retraction speed, and the Marangoni flow induced by anchoring energy gradients. Compared to certain earlier computational studies on the same problem, our predicted defect structures are more physically realistic, admitting defects of half-integer winding numbers in the bulk of LC.

Even though our model offers several advantages over previous descriptions in [START_REF] Yue | Interfacial forces and marangoni flow on a nematic drop retracting in an isotropic fluid[END_REF][START_REF] Liu | Variational approach in two-phase flows of complex fluids: transport and induced elastic stress[END_REF][START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF][START_REF] Zhou | The rise of newtonian drops in a nematic liquid crystal[END_REF][START_REF] Zhou | Dynamic simulation of droplet interaction and self-assembly in a nematic liquid crystal[END_REF], it still has the limitation of high computational cost associated with the fine mesh for the thin interfacial region, and the need to solve the fourth-order CH equation. New schemes proposed for general gradient dynamics problems, including CH diffusion, have provided promising prospects for acceleration of our calculation [e.g., [START_REF] Guillén-González | Linear unconditional energy-stable splitting schemes for a phase-field model for nematic-isotropic flows with anchoring effects[END_REF][START_REF] Shen | A new class of efficient and robust energy stable schemes for gradient flows[END_REF].

Lastly, our model can be extended to treat moving contact lines thanks to the energetic nature of the PF method (see Sec. II B). The finite-element approximation allows for complex geometries and one can further equip the model with adaptive mesh refinement or moving mesh capabilities. Our model is well-suited for further development to tackle more complex problems such as particle interactions under coupled elastocapillary forces at the interface of an LC [e.g., [START_REF] Oettel | Effective interactions of colloids on nematic films[END_REF][START_REF] Jeridi | Colloidal particles in thin nematic wetting films[END_REF][START_REF] Liu | Elastocapillary interactions on nematic films[END_REF].

FIG. 1 .

 1 FIG. 1. (a)Setup of the drop retraction problem. In 2D planar geometry it is in the x-y plane, while in 2D axisymmetric geometry it is in the r-z plane. (b) Mesh used in a typical simulation where a nematic drop retracts in an isotropic medium in a 2D planar geometry. In the inset, the solid black lines show levelsets of φ = 0.9, φ = 0 (nominal interface) and φ = -0.9, and the shaded area (bottom left) represents the defect core (q ≤ 0.5).

FIG. 2 .

 2 FIG.2. Snapshots of a nematic drop retracting in an isotropic host with planar anchoring in 2D planar geometry (Case 1 in Fig.3). The grey scale shows the concentration-weighted scalar order parameter, 1-φ 2 q, with qe = 0.81. The nominal interface φ = 0 is marked with a solid black curve. The dark region inside the drop signals the defect. (a) Well-aligned state at initial condition. (b) As the drop retracts, a defect starts to nucleate at the pole of the drop. (c) The defect core is completely established. (d ) The defect detaches from the interface and moves inside the drop. The drop shape and defect location come to a steady state.

FIG. 4 .

 4 FIG.[START_REF] Mushenheim | Straining soft colloids in aqueous nematic liquid crystals[END_REF]. Conceptual director patterns in Case 3. The thin solid curves represent streamlines of the directors, and the solid dots represent defects. In (a), the torques exerted on the drop interface by the nematic distortions are represented by semi-circles with arrows. These torques tend to rotate the interface to relieve the elastic distortion resulting in a drop elongated in the y-direction (b). Deformation is exaggerated for clarity.

FIG. 5 .

 5 FIG. 5. Steady-state drop deformation |b * s -a * s | as a function of the elastocapillary length l * ec for all the cases reported in Fig. 3. (a) 2D planar (Cases 1 -4), and (b) 2D axisymmetric (Cases 5 -8). The lines are least-square fitted to the data points.

FIG. 6 .

 6 FIG. 6. Temporal evolution of various free energies in a representative simulation of Case 8. (a) The mixing energy. (b) The LdG energy. (c) The distortional elastic energy. (d )The anchoring energy. (e) The total energy. We have used a smaller surface tension in this run (l * ec = 0.033) than the base value, so as to magnify the drop deformation for ease of visualization and analysis. Steady state is reached near t * ≈ 200. For clarity here we only show the range up to t * = 100.

FIG. 7 .

 7 FIG. 7. The dimensionless distance d *s between the defect and the interface as a function of l * ec for 2D planar cases. The insets show the distribution of 1-φ 2 q at steady state for two scenarios in Case 2. The grey scale is the same as in Fig.2. The dark region outside the drop marks the defect.

FIG. 8 .

 8 FIG. 8. (a) Defect location d * s and (b) drop deformation |b * sa *s | as functions of l * e for Cases 1 and 2. The insets show the distribution of 1-φ 2 q at steady state for two scenarios in Case 2. The grey scale is the same as in Fig.2. The dark region outside the drop signals the defect.

FIG. 9 .

 9 FIG. 9. Time evolution of the drop deformation factor D compared with a Newtonian drop retracting in a Newtonian host (N/N). (a) 2D planar (Cases 1 -4). (b) 2D axisymmetric (Cases 5 -8). We only show dimensionless time up to t * = 10 to highlight the interfacial motion. The late-time evolution is dominated by defect movement while the drop shape changes little.

FIG. 10 .

 10 FIG. 10. Anchoring energy density and velocity field for Case 8 near the end of the retraction (t * = 183). The grey scale shows the dimensionless anchoring energy density f * a . The interface and the defect core (q ≤ 0.5) are marked with solid curves. The arrows show the velocity field, displaying the Marangoni flow driven by the anchoring energy gradient.

ACKNOWLEDGEMENTS

The authors thank Julia Yeomans (Oxford, UK), Colin Denniston (Western U, Canada), Géza Tóth (UPV/EHU, Spain), Timm Treskatis (TU Dortmund, Germany), Pengtao Yue (Virginia Tech, USA), Francis Corson (ENS Paris, France), Anjishnu Choudhury (IIT Hyderabad, India), as well as all the members of the Feng-Keshet research group at UBC for valuable scientific discussions. We are grateful for financial support from the EU Marie-Curie fellowship "CoPEC" under Grant No. 794837-H2020-MSCA-IF-2017, the IdEx program of the Université de Bordeaux entitled "Développement des carrières -Volet personnel de recherche", the NSERC Discovery Grant under No. 2019-04162, and the UBC Four-Year Doctoral Fellowship. We acknowledge CMC Microsystems for software licensing.

Appendix: Parameters for drop retraction

The base parameters used for drop retraction are close to those for 5CB [START_REF] Hung | Anisotropic nanoparticles immersed in a nematic liquid crystal: defect structures and potentials of mean force[END_REF][START_REF] Ravnik | Landau-de gennes modelling of nematic liquid crystal colloids[END_REF][START_REF] Kim | Interfacial tension of a nematic liquid crystal/water interface with homeotropic surface alignment[END_REF][START_REF] Rai | Interfacial tension of liquid crystalline droplets[END_REF]. The numerical parameters for the CH dynamics (capillary width , mobility γ) are chosen based on considerations in Sec. II C and [START_REF] Yue | A diffuse-interface method for simulating two-phase flows of complex fluids[END_REF][START_REF] Yue | Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing[END_REF][START_REF] Yue | Spontaneous shrinkage of drops and mass conservation in phase-field simulations[END_REF]. They are, in dimensional form, as follows.