Viscous fingering in liquid crystals: Anisotropy and morphological transitions

We show that a minimal model for viscous fingering with a nematic liquid crystal in which anisotropy is considered to enter through two different viscosities in two perpendicular directions can be mapped to a two-fold anisotropy in the surface tension. We numerically integrate the dynamics of the resulting problem with the phase-field approach to find and characterize a transition between tip-splitting and side-branching as a function of both anisotropy and dimensionless surface tension. This anisotropy dependence could explain the experimentally observed (reentrant) transition as temperature and applied pressure are varied. Our observations are also consistent with previous experimental evidence in viscous fingering within an etched cell and simulations of solidification.


I. INTRODUCTION
Interfacial instabilities arise in a wide variety of contexts, often of applied interest, such as dendritic growth, directional solidification, flows in porous media, flame propagation, electrodeposition, or bacterial growth [1]. Notwithstanding this disparity, there has been a search for unifying common features concerning the more fundamental problem of their underlying nonequilibrium dynamics. One such feature seems to be the role of anisotropy in determining the observed morphology.
Thus, the finding that anisotropy is necessary for the needle crystal to solve the steady state solidification problem (see e.g. [2]) and indeed a critical amount of it to stabilize its tip and (possibly) generate side-branches in related local models [3] motivated the inclusion of anisotropy in viscous fingering experiments, either by engraving the plates [4], or by using a liquid crystal [5][6][7][8]. In turn, these experiments in anisotropic viscous fingering confirmed the existence of a tip-splitting / side-branching transition controlled by anisotropy and driving force. Again, theoretical work on the solidification problem [9] and numerical integration of its (nonlocal) dynamics [10] showed the anisotropy in the surface tension to control the transition in this system together with the dimensionless undercooling. To our knowledge, no analytical or numerical work on anisotropic viscous fingering [11] has focused on this transition so far.
In this way, a picture that some kind of anisotropy can control the tip-splitting or sidebranching behavior of different systems emerged. However, it is still not clear how different these anisotropies and systems can be. On the one hand, not all kinds of anisotropy seem to control the transition. The channel walls in a viscous fingering experiment, for instance, are known to play the same role as surface tension anisotropy in free dendritic growth as far as the existence of a single finger steady state solution with surface tension is concerned. Moreover, even with isotropic surface tension, this steady finger is stable up to a certain critical amount of noise [12], but no side-branching is observed, in contrast to free dendritic growth. Above this threshold the tip splits, but again no side-branching is observed, in spite of the anisotropy due to the channel walls. It is necessary to introduce some other type of anisotropy to observe the transition to side-branching. On the other hand, some types of anisotropy not directly acting on the surface tension [13] seem to actually control the transition in some systems, as is the case of liquid crystal viscous fingering experiments, which varied mainly the anisotropy in the viscosity, as well as etched cells, where the exact effect of the grooves on the free boundary equations is unclear. Therefore a connection between surface tension anisotropy (seen to control the transition in the solidification problem both experimentally and in simulations) and other types of anisotropy clearly lacks. Here we present such a connection for the case of a simple model of a liquid crystal.
Specifically, we show that two different viscosities in two perpendicular directions (in addition to some already anisotropic surface tension) can be mapped to a two-fold surface tension anisotropy (times the rescaled original anisotropic surface tension) through a convenient axis rescaling. Moreover, we integrate the resulting problem to confirm the existence of the morphological transition also for the viscous fingering equations. The numerics use a previously developed [14] and thoroughly tested [15] phase-field model for viscous fingering.
We do find such a transition as a function of the amount of anisotropy and the value of the dimensionless surface tension itself (i.e., the driving force for fixed surface tension).
The results are consistent with experimental results in viscous fingering with a liquid crystal [5][6][7][8] and etched cells [4], and also with theory and simulations for the solidification problem [9,10].
The layout of the rest of the paper is as follows: in Sec. II we refer to the special features of liquid crystals concerning viscous fingering and present a simple model for them. We then map this model onto the basic Saffman-Taylor problem with a two-fold anisotropy in the surface tension. In Sec. III we briefly describe the phase-field model used and present the numerical results. Finally, in Sec. IV we discuss their consequences for viscous fingering with a liquid crystal and consistency with related problems.

II. MODEL
In the nematic phase of a liquid crystal its molecules are locally oriented, giving rise to anisotropy in the viscosity and surface tension. The degree of orientation depends on the proximity of the other phase(s), namely the isotropic (and for some liquid crystals the smectic), i.e., it still depends on temperature, and so does the anisotropy, mainly in the viscosity [7]. Therefore one should be able to explain the tip-splitting / side-branching transition as a function of temperature in the nematics by means of the anisotropy in the viscosity alone.
In a viscous fingering experiment, the director forms a small angle with the velocity field, except maybe for the neighborhood of the interface, where it might follow its normal direction. So, as a first approximation, one can consider that there is flow alignment, and, therefore, a velocity dependent viscosity, which would make the flow nonlaplacian. However, in the vicinity of a finger tip we can approximate the direction of the flow for that of the finger, so that we can make a minimal model with only two different viscosities: one in the direction parallel to the finger and one in the perpendicular one. More details can be found in Ref. [7].
Let us now review the formulation of the Saffman-Taylor equations to account for those two different viscosities in two perpendicular directions x and y. We will do it for the channel geometry, although the result also applies to the circular cell used in the experiments of Ref. [7] with minor changes. (i.e. two different viscosities in the radial and tangential directions would also map to the standard viscous fingering equations and the same functional dependence for the surface tension anisotropy). For the sake of generality, we consider both the displaced (1) and the injected (2) fluid to have a certain distinct viscosity (µ 1 ,µ 2 ).
As in the usual Saffman-Taylor problem, in each bulk we assume the flow to be incompressible, (where u is the fluid velocity in the reference frame moving with the mean interface at V ∞ , the injection velocity), and also Darcy's law to hold, but now for two different viscosities µ x,i , µ y,i , where i = 1, 2 stand for each fluid, u x , u y are the x, y components of u, p is the pressure, η x,i = (12/b 2 )µ x,i is an inverse mobility in the x direction, η y,i = (12/b 2 )µ y,i , in the y direction, ρ i , the density, and g ef f , the effective gravity in the plane of the channel. Also as in the usual Saffman-Taylor problem, on the interface the normal velocity is continuous and equals that of the interface,r (where r is a coordinate perpendicular to it increasing towards fluid 1 and v n its normal velocity), and the pressure has a jump given by Laplace' law, with σ(φ) the (anisotropic) surface tension and κ the interface curvature. Due to Eqs. (2.1), and (2.3), the flow can be described by a scalar field ψ, the stream function, defined even on the interface by u x = ∂ y ψ, u y = −∂ x ψ (see e.g. Refs. [14,16]). However, because of the different viscosities in the x and y directions, the problem is nonlaplacian (there is vorticity) in the bulk: To circumvent this, we rescale the x and y axis by a different factor. We also take advantage to adimensionalize the resulting equations in the same way as in Refs. [14,16], so that they can be compared to those in these references, and especially to Ref. [14] in order to generalize the phase-field model described there to the case of anisotropic viscosity. Thus, we perform the following change of variables: where tildes denote new variables, a x , a y have units of length, U * is a velocity, and W , the channel width. We find so that the flow is still incompressible and we can define a new stream functionψ = (W/U * ) [ψ/(a x a y )], which will be laplacian in the bulk if and only if [see Eq. (2.5)] the velocity field is potential in each fluid, which is now the case as long as we choose a x , a y to be such that On the interface, Eq. (2.3) will be formally unchanged as long as the choice of a x , a y is the same at both sides. Note that, according to Eq. (2.9), this implies that the ratio m ≡ η x /η y must be the same for both fluids. In an air-liquid crystal experiment, this is obviously not the case, but, in the limit in which the viscosity of the air is negligible compared to that of the liquid crystal, the anisotropic character of the air viscosity in our model becomes irrelevant. In terms of the stream function, Eq. (2.3) for the new variables then reads where s is the arclength along the interface and such thatŝ ×r =x ×ŷ.
As for ∂rψ, the boundary condition for it will be given by that forũs ≡s ·˜ u. Indeed, it will have a jump on the interface due to the fact that˜ u is not potential on the very interface [see Eq. (2.8] because of the jump inη i , which gives rise to a singular vorticity on it: and therefore, making use of Eq. (2.4), where c ≡ (η 1 −η 2 )/(η 1 +η 2 ). Now choosing a y /a 2 x = 1/W , Eq. (2.9) yields m = a y /W , and defining U * ≡ cV ∞ /m + [mg(ρ 1 − ρ 2 )] /(η 1 +η 2 ) we recover the usual result for viscous fingering in a channel (see Refs. [14,16]), , except for the m factors in the definition of U * -and therefore in B(φ)-and the clue fact that W κ and σ(φ) are still in the old variables and must be rescaled: 15) whereφ is the angle fromx tor.
To summarize, we recover the usual viscous fingering equations, including Eq. (2.13), which finally reads 16) but now with an anisotropic dimensionless surface tension of the form 17) whereB 0 is the dimensionless surface tension of isotropic viscous fingering (2.18) except for the m factors, with σ(φ) ≡ σ 0 Σ(φ). This means that, even with an originally isotropic surface tension, that in the rescaled problem has a two-fold anisotropy with a very specific form given by the last (third) factor at the r.h.s. of Eq. (2.17). On the other hand, the possible original anisotropy in the surface tension will change its functional form in the rescaled problem according tõ 2.19) and the rescaled problem will have the two-fold anisotropy of the mentioned last factor superimposed to the transformed anisotropy of Eq. (2.19) [second factor at the r.h.s. of Eq. (2.17)]. A similar result was found in a different context, namely for the nematicsmectic B transition, where two different heat diffusivities in two perpendicular directions could be mapped to the same type of anisotropy in the surface tension and the same type of transformation in the original anisotropy [17]. However, note that here the assumption is that the growth is in the direction of lowest viscosity (because of flow alignment of the director), which results in growth in the direction of largest surface tension (φ = π/2), whereas in Ref. [17] the situation is just the opposite: growth was found to be in the direction of lowest diffusivity because that is the direction of lowest capillary length. Also for isotropic diffusivities it is known that steady needle crystals can only grow in the direction of minimal capillary length [2], although there the anisotropy is assumed to be four-fold. Finally, for the described minimal model the original anisotropy in the surface tension would be two-fold, e.g.
so that the transformed anisotropy would read (2.21)

III. NUMERICAL INTEGRATION
We now integrate the rescaled problem, namely Laplace equation for the stream function with the boundary conditions Eqs. (2.10) and (2.16). In principle, given an initial condition, we should rescale it, evolve it using the rescaled dynamics, and translate the resulting interface back to the original variables, but we will not perform any rescaling, since the initial condition is free, and the tip-splitting or side-branching character of the result is unaffected by the final translation into the original variables. Instead, we will consider the rescaled problem in its own, and simulate it by means of the following phase-field model, where θ is the phase field, ǫ,ǫ are model parameters which must be small to recover the sharp-interface equations of the rescaled problem, and we have dropped the tildes of the rescaled variables. We have defined f (θ) ≡ θ(1 − θ 2 ), and γ(θ) 2 ≡ŝ(θ) · ∇B(θ)κ(θ) +ŷ , κ(θ) ≡ − ∇ ·r(θ), withr(θ) ≡ ∇θ | ∇θ| andŝ(θ) ≡r(θ) ×ẑ. This model was introduced for isotropic viscous fingering in Ref. [14] and extensively tested in Ref. [15]. From this work we know that it will yield converged results for the steady fingers (and, in particular, for their widths The only change to be made for the anisotropic case is to set B(θ) not merely equal to a constant, but to that given by Eq. (2.17) taking φ = φ(θ) = arccosx ·r(θ). This gives B(θ) = B(φ) + O(ǫ 3 ), which not only satisfies the desired sharp-interface limit, but also ensures that the introduction of anisotropy will not result in any extra first-order correction to the free boundary problem, so that the above conditions on ǫ,ǫ still hold.
The same phase-field equations could be used for the circular geometry reinterpreting the parameters, since an analogue rescaling yields formally the same result. However, the boundary conditions would change. For instance, injection at the center of the cell should also be considered. Anyhow, we choose to simulate the well controlled situation in which an (unstable) steady finger propagates in a channel of width W . This representation is of course exact for single fingers in experiments carried out in the linear Hele-Shaw cell [8], but only a (good) approximation for the vicinity of a finger in a multifinger configuration (with many fingers) in the circular geometry [7].
We investigate the transition between the tip-splitting and side-branching behaviors as both the dimensionless surface tension B 0 and its anisotropy m − 1 are varied. We are interested mainly in the effect of the anisotropy coming from the viscosity (m − 1), so we drop that in the original surface tension (σ(φ) = σ 0 ). The runs use equal viscosities in both fluids, c = 0, for reasons of numerical efficiency (ǫ = 0.2), but we do not expect the viscosity contrast to affect the stability of the tip for similar reasons for which it does not play any role in the (linear) stability of a flat interface. To check this conjecture we ran simulations with B 0 = 10 −3 both for c = 0 and c = 0.8, two values of the viscosity contrast for which a dramatic change in the competition dynamics was seen using the same phase-field model [15], and we found that the transition lies in both cases between m = 2 and m = 2.25. These c = 0.8 runs are indeed the ones shown in Figs. 1 and 2.
We use ǫ = 4 × 10 −3 , so that we can simulate accurately with values of the dimensionless surface tension down to B 0 = 4×10 −4 . We first run a steady finger with a large, isotropic dimensionless surface tension B(φ) = B 0 = 10 −2 , large enough for the finger not to destabilize for the amount of (numerical) noise we have and thus reach its steady width and velocity. Once this is achieved -see inner interface in Figs. 1(a) and 1(b)-, we perform a "quench" in surface tension, i.e., we instantly reduce it to some lower value. Simultaneously, we also introduce some amount of anisotropy m − 1.
The subsequent interface evolution for B 0 = 10 −3 (and c = 0.8,ǫ = 0.08) is also shown within the reference frame moving with the mean interface in Figs. 1(a) (m = 2) and 1(b) (m = 2.25) in the form of snapshots at time intervals 0.11. (Simulations used only half of the channel and reflecting boundary conditions at its center x = 0). The corresponding y position of the interface at the center of the channel (also in the frame of the mean interface) is plotted against time in Fig. 2.
For this value of B 0 the finger clearly destabilizes: First its tip widens and flattens (see Fig. 1) and therefore slows down (see Fig. 2) for any value of the anisotropy. (Note that for t < 0 the tip position would be a straight line in time, since the finger was steady, and, in particular, its velocity). Then, for m = 2 the tip continues to flatten and slow down until its curvature [ Fig. 1(a)] and eventually its velocity in the frame of the mean interface (lower curve in Fig. 2) reverse their signs. Finally, the velocity of the interface at the center of the channel seems to reach some negative constant value (again, in the frame moving with the mean interface) corresponding to the growth of two parallel fingers at each side. We identify this reversal of the curvature sign at the center of the finger and this always convex tip position vs. time plot with the tip-splitting morphology.
In contrast, for m = 2.25 the reversal of the curvature sign takes place at some distance of the center of the channel, while at the center the curvature increases again [ Fig. 1(b)] and makes it possible for the tip to speed up again as well, giving rise to a change of concavity in the tip position vs. time plot (upper curve in Fig. 2). We identify this reversal of the curvature sign at a distance from the center of the channel and this change of concavity in the tip position vs. time plot with the side-branching morphology.
In this way we systematically explore values of the dimensionless surface tension B 0 ranging from B = 10 −2 down to B = 4 × 10 −4 . For each value of B 0 we simulate with several values of the anisotropy m − 1, and we find that there is a relatively sharp transition between the tip-splitting and side-branching morphologies. In Fig. 3 we show for each value of B 0 (x axis) the two closest values of m−1 (y axis) for which the two different morphologies are observed, namely tip-splitting (circles) and side-branching (triangles). Thus we know that the transition line must lie somewhere between the circles and the triangles, and that above (larger values of m − 1) and left (lower values of B 0 ) of that transition line the morphology is side-branching, and below and right of it, tip-splitting. This means that the critical anisotropy m − 1 above which side-branching replaces tip-splitting decreases with decreasing dimensionless surface tension B 0 .
In fact, this critical anisotropy vanishes at B 0 ∼ 5 × 10 −4 , and below this value only side-branching is observed, even if one uses negative anisotropies down to m − 1 = −0.9, which correspond to a viscosity larger in the direction of growth of the finger than in the perpendicular one, and which is not the case of the liquid crystal experiments that motivated this study. (m − 1 > −1 to keep the two viscosities and therefore B(φ = 0) finite and positive). Of course, the specific value of B 0 for which the critical anisotropy vanishes could be affected by the fact that a residual (four-fold) grid anisotropy remains, but it seems unavoidable that there is such a (finite) value of B 0 , since the transition line curves down as B 0 is decreased, and for large enough values of B 0 or the anisotropy m − 1 the grid spacing ∆x = ǫ = 4 × 10 −3 is far too fine to affect the effective anisotropy.
On the other hand, for B 0 ≥ 1.4 × 10 −3 and for the time elapsed in our runs no clear side-branching is actually observed above the transition line extrapolated from lower values of B 0 , whereas tip-splitting still occurs below that line. For even larger values of B 0 , namely B 0 ∼ 2 × 10 −3 , not even tip-splitting is observed again within the time elapsed, although the steady finger still destabilizes through the widening and flattening of its tip. Finally, for B ≥ 10 −2 the finger is completely stable for the amount of noise we have, as was pointed out before.

IV. DISCUSSION AND CONCLUSIONS
We have shown that viscous fingering with two different viscosities in two perpendicular directions maps to the standard viscous fingering equations (i.e., for isotropic viscosity) with an extra two-fold anisotropy in the surface tension, which is such that, together with the hypothesis of flow alignment of the director, leads to growth in the direction of maximal surface tension. We have simulated the resulting problem using a previously developed phase-field model [14,15], and we have found that there is a transition from the tip-splitting to the side-branching morphology as either the anisotropy in the surface tension is increased or the dimensionless surface tension is decreased. We now draw the connection with the liquid crystal experiments of Ref. [7].
The observed anisotropy dependence is consistent with the experimental finding that there is a transition from tip-splitting to side-branching and back to tip-splitting with temperature in the nematic phase, since close to the other phases the director alignment, and consequently the anisotropy, weaken [7]. The transition is found to be also reentrant with injection pressure [7], which is explained there with the hypothesis that too low pressures do not achieve flow alignment, whereas too large ones break down the Hele-Shaw approximation because of the importance of inertial terms in the hydrodynamic equations, which then destroy the flow alignment again. This anisotropy dependence is also consistent with simulations of the boundary layer model [3] and the full solidification problem [10], as well as with analytical approaches to solidification [9].
As for the dimensionless surface tension B 0 dependence, one first needs to relate the values of B 0 used in the channel simulations to the experimental parameters in the circular geometry. To do this we consider a virtual channel whose walls are placed at half the distance between a finger and its nearest neighbors. The channel width W is given by this distance between adjacent finger tips, whereas the effective injection velocity V ∞ turns out to be the ratio between the injection pressure and R, the mean distance between a tip and the injection point. Then, the following dynamic picture of a typical experiment in the circular cell emerges: Initially some fingers develop. If the anisotropy, m − 1, is strong enough, their tips are stable (which corresponds to a point above the transition line in Fig. 3, where the observed morphology is side-branching). As these fingers grow radially, W increases as R, whereas the effective driving force, (η 1 − η 2 )V ∞ , decreases as 1/R, so that the dimensionless surface tension B 0 they experience is found to decrease as 1/R. Thus, the corresponding point in Fig. 3 moves to the left (lower values of B 0 ), and the side-branching behavior is preserved.
In contrast, if the anisotropy m−1 is not strong enough, the tips split (the corresponding point is below the transition line in Fig. 3, where the observed morphology is tip-splitting). As a result the number of fingers increases, which then compensates for the growth of the distance between finger tips as R in such a way that the effective dimensionless surface tension B 0 keeps roughly steady during the pattern development, so that the corresponding point in Fig. 3 basically does not move. Thus the transition line is not crossed and the tip-splitting behavior is also preserved. In this case, we can estimate B 0 in the experiments to be of order 10 −3 , whereas m is said to be around 2 at the transition [7]. This is indeed very close to our transition point B = 10 −3 , 2.05 ≤ m ≤ 2.1 in Fig. 3, but it should be taken into account that the value of B 0 below which fingers destabilize is known to depend on the amount of noise present [12]. Accordingly, one expects the whole transition curve to be shifted to different values of B 0 for a different amount of noise. The latter, however, is uncontrolled both in the experiments and the simulations.
On the other hand, the dependence of the value of the anisotropy at the transition on the driving force has been seen in viscous fingering with an etched cell [4]. However, we find the critical amount of anisotropy for side-branching to replace tip-splitting to vanish at a finite value of B 0 . Below this value we only observe side-branching. A more realistic model of viscous fingering in a liquid crystal should include a velocity dependent viscosity. In general the resulting nonlaplacian character of the problem could not be avoided, but in principle it would still be possible to simulate the dynamics by means of a phase-field model. Fig. 1 Destabilization of the tip of a (stationary) finger after instantly decreasing B 0 to B 0 = 10 −3 at the time of the first interface shown. Successive interfaces are shown in the reference frame moving with the mean interface at time intervals 0.11, for c = 0.8,ǫ = 0.08. The latest interface is represented in bold. (a)Tip-splitting for m = 2. (b) Side-branching for m = 2.25. Fig. 2 y coordinate of the interface at the center of the channel (x=0) in the reference frame of the mean interface as a function of time corresponding to Figs. 1(a) (lower curve) and 1(b) (upper curve). t = 0(0.88) corresponds to the first (last) interface shown there. Fig. 3 Transition between tip-splitting (circles) and side-branching (triangles) as a function of the surface tension anisotropy m − 1 and the dimensionless surface tension B 0 .