Generalized Berreman's model of the elastic surface free energy of a nematic liquid crystal on a sawtoothed substrate

PHYSICAL REVIEW E 86, 041706 (2012) Generalized Berreman’s model of the elastic surface free energy of a nematic liquid crystal on a sawtoothed substrate O. A. Rojas-G´ omez and J. M. Romero-Enrique Departamento de F´ ısica At´ omica, Molecular y Nuclear, Area de F´ ısica Te´ orica Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla, Spain (Received 18 May 2012; revised manuscript received 5 October 2012; published 26 October 2012) In this paper we present a generalization of Berreman’s model for the elastic contribution to the surface free-energy density of a nematic liquid crystal in presence of a sawtooth substrate which favors homeotropic anchoring as a function of the wave number of the surface structure q, the tilt angle α, and the surface anchoring strength w. In addition to the previously reported nonanalytic contribution proportional to −qln q, due to the nucleation of disclination lines at the wedge bottoms and apexes of the substrate, the next-to-leading contribution is proportional to qfor a given substrate roughness, in agreement with Berreman’s predictions. We characterize this term, finding that it has two contributions: the deviations of the nematic director field with respect to a reference field corresponding to the isolated disclination lines and their associated core free energies. Comparison with the results obtained from the Landau-de Gennes model shows that our model is quite accurate in the limit wL > 1, when strong anchoring conditions are effectively achieved. DOI: 10.1103/PhysRevE.86.041706 PACS number(s): 61.30.Dk, 61.30.Hn, 61.30.Jf I. INTRODUCTION The behavior of nematic liquid crystals in the presence of microstructured substrates has been the subject of intensive research in recent times [1–3]. This problem has practical applications such as the design of zenithally bistable devices [4–9] or the trapping of colloidal particles in specified sites [10,11]. It is well known that the nematic director field, in the presence of a structured substrate, may be distorted, leading to an elastic contribution to the free energy. Since the seminal work by Berreman [12,13], this problem has been extensively studied and generalized in the literature [4,14–27]. Wetting and filling transitions by nematic on these grooved surfaces [28–31] have also been studied. When the substrate presents cusps, topological defects nucleate on them [14–16,26,27], and Berreman’s expression of the elastic contribution to the freeenergy density, which is proportional to the wave number of the substrate structure q, breaks down because of the emergence of a nonanalytical contribution proportional to −qln qassociated to the nucleated defects [27]. This result contrasts with the phenomenology observed in smooth substrates, as sinusoidal substrates, in the absence of topological defects. In this situation a suitable generalization of Berreman’s model may be applied [24,31]. In this paper we will complete the characterization of the elastic contribution to the surface free-energy density for sawtoothed substrates [27]. Beyond the −qln qterm previously mentioned, we find that the next-to-leading contribution follows Berreman’s scaling with the wave number q.The origin of this term is twofold: (i) the deviations of the nematic director field with respect to the distortions imposed by the presence of the nucleated disclination lines and (ii) the defect core contributions. We estimate both contributions, finding a fairly good agreement with the reported values in Ref. [27]. The paper is organized as follows. The problem is set up in Sec. II, where we identify the different contributions to the elastic contribution to the surface free-energy density. Secs. III and IV are devoted to the estimation of these contributions, and the obtained results will be discussed in Sec. IV.The conclusions are presented in Sec. V. II. THE MODEL We consider a nematic liquid in contact with a sawtooth substrate characterized by the angle αand the length side L (see Fig. 1). The substrate favors homeotropic anchoring of the molecules. We assume translational symmetry along the z axis and a periodic distribution of wedges and cusps along the xaxis. Under these conditions, the nematic director field n(r) shows only in-plane distortions [27], so it can be parametrized by using the angle θbetween the local director and the y axis, yielding n(r)=(−sin θ(r),cos θ(r),0). Far from the substrate, the bulk nematic phase orients homogeneously along either the xaxis (perpendicular texture N⊥)orthey axis (parallel texture N). The nematic order of the system can be locally represented by a traceless symmetric orderparameter second-rank tensor Q, with Cartesian components Qij =3 2S[ninj−1 3δij ]+1 2B[lilj−mimj], where Sis the nematic order parameter, which measures the orientational ordering along the nematic director, and Bthe biaxiality parameter, which measures the ordering of the molecules on the orientations perpendicular to n, characterized by the eigenvectors land m. As in previous works [27,29,30], the system will be described within the Landau-de Gennes (LdG) framework. The order-parameter tensor profile is obtained by minimizing the surface free-energy density functional FLdG =V (φbulk +φel)dV +A φsurf ds, (1) where φbulk is the bulk free-energy density, φel is the elastic free-energy density, and φsurf is the surface free energy, defined as Ref. [13] φbulk =aTr Q2−bTr Q3+c[Tr Q2]2,(2) φel =L1 2∂kQij ∂kQij +L2 2∂jQij ∂kQik,(3) φsurf =−2 3¯wTr[Q·Qsurf],(4) 041706-1 1539-3755/2012/86(4)/041706(10) ©2012 American Physical Society O. A. ROJAS-G ´ OMEZ AND J. M. ROMERO-ENRIQUE PHYSICAL REVIEW E 86, 041706 (2012) α Rc Rc c R II Substrate III III II II L I Nematic η C1 x y FIG. 1. Schematic picture of the geometry of the system, characterized by the side length Land the angle α. The different regions I, II, and III are outlined. See text for explanation. where adepends linearly on the temperature, band c are positive constants, and L1and L2are positive parameters related to the elastic constants. This model shows a bulk nematic-isotropic transition for a=b2/24c(i.e., T=TNI), with a nematic order parameter Sb(T)=b/8c+ (b/8c)1−64ac/3b2below the transition. Note that the nematic order parameter in the nematic phase in coexistence with the isotropic phase is Sb(TNI)=b/6c.Inthe nematic phase, a natural length is given by the nematic correlation length ξalong the local nematic director axis, which can be expressed in terms of the LdG model parameters as ξ=√(3L1/2+L2)/(9bSb/4−6a). Finally, the free energy is reduced by ξS2 b(3L1+2L2)/√2. Regarding the surface contribution, homeotropic alignment of the nematic is favored by setting Qsurf =Sref surf(3ν⊗ν−1)/2, with ν= ((x/|x|)sinα, −cos α,0) being the outwards normal vector to the substrate and Sref surf the favored nematic order parameter on the substrate. Note that Eq. (4) corresponds to the linearization on Qof the Nobili-Durand surface free energy [32]. This model favors a perfect nematic ordering on the surface (i.e., S=1), unlike the original Nobili-Durand surface contribution, which favors S=Sref surf. However, if wand Sbare small enough, the results obtained with both surface free energies are qualitatively similar. The nematic-substrate coupling Eq. (4) has been used to study the orientational wetting phenomena on flat [33,34] and rough substrates [29–31,35]. As we are interested in surface phenomena, instead of minimizing directly the LdG functional Eq. (1), we subtract the contribution of an unperturbed bulk nematic phase to get the excess free energy and divide the result by the projection area of the substrate in the x-zplane. In this way, and measuring all the quantities by their natural units described above, we get the following free-energy density functional: F=1 λλ/2 −λ/2 dx ∞ |x|tan α dy3−2τ 4τ−32Tr ˜ Q2 3−1 −2 4τ−34 3Tr ˜ Q3−1+τ 4τ−34 9[Tr ˜ Q2]2−1 +1 3+2κ[∂k˜ Qij ∂k˜ Qij +κ∂j˜ Qij ∂k˜ Qik] −1 λλ/2 −λ/2 dx2 3wTr[ ˜ Q·˜ Qsurf]|y=|x|tan α,(5) where λ≡2Lcos αis the substrate periodicity wavelength, ˜ Q=Q/Sb, where Sb=Sb(T)isthebulkvalueofthe nematic order parameter at the temperature Tand τ= Sb(T)/Sb(TNI) is the ratio between the nematic order parameter at the temperature Tand the nematic-isotropic transition temperature TNI. The positions and lengths are measured in units of √2ξ. Finally, κ=L2/L1(κ>−3/2) and w=√8¯wξ[Sref surf/Sb(T)]/(3L1+2L2) the reduced anchoring strength and ˜ Qsurf =(3ν⊗ν−1)/2. Hereafter, we consider every quantity in its corresponding natural units. The global minimum of the functional Eq. (5) yields the mean-field equilibrium surface free-energy density, f,in reduced units. The contribution due to the elastic deformations induced by the substrate structure, felastic, can be obtained as felastic =f−rσNW(w), where the surface roughness is r= 1/cos αand σNW(w) is the nematic-flat substrate interfacial tension. In a similar way as it was obtained at nematic isotropic coexistence in Ref. [31], the interfacial tension for the LdG model has the expression σNW =√2(τ−1)(˜ S(0) +1)2+˜ S(0)2 6τ2√4τ−3[2τ2˜ S(0)2−τ(2τ+1) ˜ S(0) −4τ2+7τ−3] −√2√4τ−3(−4τ2+6τ−3) 6τ2√4τ−3 −√21−3τ+2τ2 2τ5/2√4τ−3ln τ˜ S(0) +τ−1+√τ(τ−1)(˜ S(0) +1)2+˜ S(0)2 2τ−1+√τ√4τ−3−w˜ S(0),(6) where ˜ S(0) =S(0)/Sband S(0) is the nematic order parameter at the substrate. The value of ˜ S(0) is obtained as the largest solution of the equation (3 −2τ)˜ S2(0) −2˜ S3(0) +τ˜ S4(0) =1−τ+4τ−3 2w2.(7) At coexistence τ=1 and Eq. (6) reduces to σNW =√2(2˜ S(0) +1)( ˜ S(0) −1)2 6−w˜ S(0),(8) where ˜ S(0) =(1 +1+2√2w)/2. A systematic study of this system via full minimization of the LdG functional was done in Ref. [27]. For this purporse, the functional Eq. (5) was numerically minimized by using a conjugate-gradient method. The numerical discretization of the continuum problem is performed with a finite-element method [36] combined with adaptive meshing in order to resolve the different length scales that may emerge in the problem [37]. It was found that the N⊥texture has lower free energy if α<π/4 due to lesser distortion. Conversely, the Ntexture has lower free energy for α>π/4, in agreement with earlier predictions [14–16,26]. For large wL, strong anchoring conditions are effectively achieved, leading to the 041706-2 GENERALIZED BERREMAN’s MODEL OF THE ELASTIC ... PHYSICAL REVIEW E 86, 041706 (2012) nucleation of disclination lines characterized by non-half- integer winding numbers along the ridges and wedges of the substrate [14–16,26,27]. As a consequence, the elastic contribution to the surface free-energy density has the following scaling [27]: felastic ≈−K(α) 2πqln qcos α π+q 2πB(α,w),(9) where q=2π/λ =π/Lcos αis the substrate periodicity wave number and K(α) is defined as K(α)=⎧ ⎨ ⎩ Kπα2 (π 2)2−α2N⊥texture Kπ π 2−α π 2+αNtexture (10) with K=(9/2)(2 +κ)/(3 +2κ) being the reduced bulk elastic constant associated to bend and splay distortions. From the numerical results, the function Bis found to depend on the substrate roughness (i.e., α) and nematic texture, as well as the anchoring w, but asymptotically not on Lfor large wL. However, as Lincreases, the complete minimization becomes very time-consuming. On the other hand, we do not get information about the origin of B(α,w). A closer inspection of the typical nematic textures obtained from minimization of the LdG model for large wL allows us to identify three different spacial regions with distinctive behavior in the nematic tensor field. Far from the substrates, the nematic order parameter Stakes the bulk value everywhere, but the orientational field θdeforms smoothly in order to satisfy the preferred anchoring conditions on the walls. Second, close to the walls but far away from apexes and wedges, there is a layer in which the orientation field is nearly homogeneous but Schanges in a scale of order of the nematic correlation length. In fact, the obtained nematic order parameter profile is basically equal to that obtained for a planar wall with the same anchoring conditions. Finally, around the apexes and wedges, both Sand θare inhomogeneous. This explains the scaling of the surface free-energy density reported in Ref. [27]interms of the substrate wavelength λ. Thus, the leading contribution to fis σNW due to the fact that the anchoring is locally homeotropic on most of the substrate and that Sbehaves as if the substrate were flat. On the other hand, the elastic deformations far away from the substrate can be described by the Frank-Oseen functional [27], which leads to a contribution of order of λ−1to fin absence of topological defects and of order of ln λ/λ when disclination lines are present [29]. In both cases this contribution cannot overcome the leading order due to the surfaces if Lis large enough. Finally, we expect that the regions around the apexes and wedges will generate a line tension term, which contributes as λ−1to f. Inspired by these observations, in this paper we will introduce an alternative way to obtain felastic from the functional Eq. (5). We divide the minimization domain into three regions (see Fig. 1): most of the domain will correspond to region I, formed by the points which are far enough from the substrate. The neighborhood of the substrate will be split into two regions: region II, formed by the union of the circular sections of radii ξ<R cLcentered at each wedge and apex, and region III, which are the points which are at a distance smaller than η∼ξto the substrate, but at distances larger than Rcfrom any substrate ridge. Our hypothesis is that, for large wL,the minimization of the surface free-energy functional restricted to each region (subject to appropriated boundary conditions) gives a good approximation to the complete minimization of F. On the other hand, we anticipate that this analysis will give us some insight in the different contributions to B(α,w). However, the success of our approach can be assessed only when comparing their predictions with the results obtained from the full minimization of the LdG model, as will be discussed in Sec. V. We start with the minimization of region III. As argued in Ref. [27], large wL leads to strong anchoring conditions. This condition is equivalent to impose that Lis larger than the extrapolation length, which is of order of K/w in reduced units. So, in order to minimize the surface free-energy density, we impose that the angle field θ(r) is constant along its boundary, so the nematic director field is homogeneous and equal to the normal to the substrate. Consequently, the minimization of the free-energy functional in this region will lead to a homogeneous director field normal to the substrate, although the nematic order parameter Sat each point will depend on its distance to the substrate. We impose the following fixed boundary conditions for S:S=Sbat the boundary between regions I and III and the equilibrium nematic order parameter profile for the flat wall case at the boundary between regions II and III. Assuming that ηis large enough, this situation is completely equivalent to the flat case, so the minimum value of the surface free-energy density in this region, fIII, will be σNW(1 −2Rc/L)/cos α, with corrections of order of exp(−η/ξ). The next sections will be devoted to the evaluation of the minimum values of the surface free-energy functionals at the remaining regions, fIand fII. III. EVALUATION OF fI The variations of Sare restricted to the neighborhood of the substrate of a width typically of order ξand inside the defect cores. So, in region I, Stakes the bulk value Sbeverywhere [27]. Thus, the surface free-energy functional to minimize in region I reduces to a Frank-Oseen functional, FI≈K 2λI dxdy|∇θ|2,(11) where the integration is restricted to region I, Kis the reduced elastic constant, and θis the orientation field. The minimization of Eq. (11) yields to the Laplace equation for θ, ∇2θ(r)=0. In the far field, we impose Dirichlet boundary conditions limy→∞ θ(r)=α∞, where α∞=0fortheN⊥ texture and α∞=π/2fortheNtexture. Along the contours x=±λ/2, we should impose periodic boundary conditions. However, we impose instead Dirichlet boundary conditions θ=α∞along these contours, as we know from the full LdG model minimization that these are the conditions satisfied by the mean-field solution [27]. Finally, we assume strong anchoring conditions along the boundary between regions I and III: θ(x,y =x2tan α/|x|,z)=α∞+(x/|x|)(α−α∞). As discussed above, this condition will be accurate if wL 1. Their contribution to the free-energy density, fI, comes from a contour integration of the mean-field solution via [27], fI=K(α−α∞) λC1 ν·∇θds,(12) 041706-3 O. A. ROJAS-G ´ OMEZ AND J. M. ROMERO-ENRIQUE PHYSICAL REVIEW E 86, 041706 (2012) where C1is the contour parallel to the boundary between regions I and III between a wedge and apex (see Fig. 1). As argued in Refs. [14–16,27], the presence of cusps in the substrate induces the formation of disclination lines, which are responsible for the non-Berreman scaling of the elastic contribution to the surface free-energy density. In this section, we will complete that analysis, evaluating the next-to-leading contribution. A. Singular contribution The solution θ(r) to the Laplace equation subject to the boundary conditions mentioned above can be split into two terms: a singular contribution θs(r), due to the periodic array of disclination lines nucleated at the ridges of the substrate, and a nonsingular contribution, θns (r). A representation of the singular contribution for each texture is given by Ref. [27] θ⊥ s=−α π 2−α−arctan tan qx 2coth qy 2+arctan tan qx 2 +α π 2+α−arctan tan qx 2tanh q(y−Lsin α) 2+arctan tan qx 2 (13) θ s=π 2+−arctan tan qx 2coth qy 2+arctan tan qx 2 − π 2−α π 2+α−arctan tan qx 2tanh q(y−Lsin α) 2+arctan tan qx 2.(14) Their contribution to the free-energy density, fs I, comes from a contour integration of these solutions, fs I=K(α−α∞) λC1 ν·∇θsds.(15) In Ref. [27] the large-Lbehavior of Eq. (15) was estimated, leading to the non-Berreman term K(α)qln(L/Rc). However, after some algebra, it is possible to obtain explicitly fs Ifrom the solutions Eqs. (13) and (14) via Eq. (15) as fs I=K(α)q 2π−ln qRc+ln 2 cosh π 2tan α−αtan α, (16) where K(α) depends on the texture and substrate geometry as Eq. (10) and we neglected terms of order q3(Rc)2and qη. As the dependence on the nematic texture comes from K(α), fs Iwill be minimum for the N⊥(respectively, N)texturefor α<π/4 (respectively, α>π/4). Two remarks are pertinent at this point. First, we note that although fs Imay depend on the explicit representation of the singular solution, the leading non-Berreman contribution is independent on this representation. The reason for this is that this leading contribution arises from the behavior close to the wedges and apexes of θs, which must asymptotically approach to the corresponding to an isolated disclination line [27]. Second, the next-to-leading contribution gives a first contribution to B(α,w), which we will denote as BI,s(α), given by the expression BI,s(α)=K(α)ln 2 πcosh π 2tan αcos α−αtan α. (17) B. Nonsingular contribution The nonsingular part of the mean-field solution, θns,isa solution of the Laplace equation ∇2θns =0, subject to the boundary conditions θns =0 in the far field, i.e., y→∞ and along the boundaries x=±λ/2. On the other hand, θns =α−θsalong C1. We do not have an explicit expression for θns (however, there is an implicit expression via a Schwarz- Christoffel transformation [14], see below), so we have to resort to numerical methods. We have used two different techniques: a finite-element method, analogous to the method outlined in Sec. II to solve the LdG model but minimizing, instead, the functional Eq. (11) subject to the boundary conditions for θns mentioned above, and, as an alternative, the boundary-element method [38,39]. In this technique the solution θns inside region I can be written as θns(r)=∂I ds([ν(s)·∇sθns(s)]G(s,r) −θns(s)[ν(s)·∇sG(s,r)]),(18) where the contour integral over the boundary ∂I of region I is counterclockwise, ν(s) is the outwards normal to the boundary at s, and G(s,r) is the fundamental solution of the Laplace equation G(s,r)=−ln (|s−r|)/2π. As we impose Dirichlet boundary conditions, the second term in the right-hand side of Eq. (18) is known. On the other hand, the normal derivative of θns at the boundary is obtained by solving the integral equation [38,39] ∂I ds[ν(s)·∇sθns(s)]G(s,r) =θns(r) 2+∂I dsθns(s)[ν(s)·∇sG(s,r)],(19) where now r∈∂I. In order to solve Eq. (19), we discretize the boundary in a set of straight segments (the boundary elements). 041706-4 GENERALIZED BERREMAN’s MODEL OF THE ELASTIC ... PHYSICAL REVIEW E 86, 041706 (2012) H L (a) (b) −1 0 1 1 ε2 ε1 ε R R Rc c c FIG. 2. (Left) Minimization cell for the evaluation of fI,ns. (Right) Mapping of the minimization cell under the Schwarz-Christoffel transformation in the limit H→∞. We use the constant boundary-element approach [39], so we assume that both θns and its normal derivative are constants along each boundary element. Introducing this approximation to Eq. (19), we obtain a set of linear algebraic equations for the normal derivatives of θns. Once we solve this set of equations, and introducing the same approximation in Eq. (18), we obtain the nonsingular orientational field θns inside region I. Once θns is evaluated, its contribution to the surface freeenergy density can be obtained from fns I=K(α−α∞) λC1 ν·∇θnsds,(20) =K 2λI dr|∇θns|2+K λC1 [α−θs(s)]ν·∇θsds,(21) where the first result is more appropriated for the boundaryelement technique, while the second is more appropriate for the finite-element method [note that the last term in the second result can be evaluated numerically with high accuracy by standard methods as we know analytically θs(r)]. The numerical minimization is performed in the cell shown in Fig. 2. In order to minimize finite-size effects, the cell height His taken to be at least 4 times the value of L(note that θs decays exponentially to α∞for y2Lcos α/π). We checked that the value of λf ns Iis independent of λ, as expected, so the nonsingular contribution to B(α,w), BI,ns(α), is related to fns I via BI,ns(α)=λf ns I. The numerical results are shown for both the N⊥and Ntextures in Fig. 3. The agreement between the results obtained from the finite-element method and boundary-element method is excellent, as for the evaluation of θns(r), which takes non-negligible values only above the contour C1and vanishes close to the wedges and apexes (see Fig. 4). Our results show that the nonsingular contribution corresponding to the N⊥(respectively, N) texture is smaller than the contribution associated to the N(respectively, N⊥) texture for α<π/4 (respectively, α>π/4). C. Exact evaluation of fI In the previous paragraphs we have outlined how to obtain θ(r)=θs(r)+θns(r) and, from that, to obtain fI=fs I+fns I, as well as BI(α)=BI,s (α)+BI,ns(α). However, it is possible to obtain fIdirectly without knowing the explicit form of θ(r) 0π/8 π/4 3π/8 π/2 α(rad) 0 0.5 1 λfI,ns/K N⊥ (FEM) N⊥ (BEM) N|| (FEM) N|| (BEM) FIG. 3. Plot of λfI,ns/K as a function of αfor the N⊥and N textures by use of the finite-element method (FEM) and the boundaryelement method (BEM). by using a Schwarz-Christoffel transformation [14,40], z=d¯ ζC (˜ ζ+1)1/2−α/π ˜ ζ2α/π (˜ ζ−1)1/2−α/π =Cζ1−2α π2F11 2−α π,1 2−α π,3 2−α π,ζ2+C,(22) where z=x+iy,2F1(a,b,c,z) is the Gauss hypergeometric function, and C,C, and C are complex constants. This conformal transformation maps the minimization cell for H→∞ into the upper half ζplane (see Fig. 2), transforming the origin into the origin and the edges z=±Lcos α+iLsin αinto ζ=±1, respectively. These conditions fix the values of Cand C, so the Schwarz-Christoffel transformation finally reads z=Leiα 3 2−α π1 2+α πζ1−2α π ×2F11 2−α π,1 2−α π,3 2−α π,ζ2,(23) where (x) is the complete gamma function. Equation (23) can be formally inverted, so ζ=ζ(z/L;α)=x(x/L,y/L;α)+ iy(x/L,y/L;α). We will consider the limit η/L →0 and Rc/L small but finite. In this approach, the boundary of zone II 0.00 0.05 0.10 0.15 0.20 x 0.0 0.2 0.4 y 0.00 0.01 0.02 0.03 θns x,y FIG. 4. (Color online) Plot of θns as a function of xand yobtained from the finite-element method for α=π/6andL=0.25. For the sake of clarity it is only represented in half a minimization cell 0 < x<Lcos α. 041706-5 O. A. ROJAS-G ´ OMEZ AND J. M. ROMERO-ENRIQUE PHYSICAL REVIEW E 86, 041706 (2012) becomes, under the Schwarz-Christoffel transformation, the real axis in the ζplane, rounded around ζ=0 and ±1. The expansion of Eq. (23) around these values show that the circles of radii Rcaround the origin and the edges z=±Lcos α+ iLsin αmap into circles [up to corrections of order of (Rc/L)2] of radii 1for ζ=±1 and 2for ζ=0 given by 1=1 21+2α π 1−2α π1 1 2+α πRc L3 2−α π1 2+α π 1 1 2+α π 2=Rc L3 2−α π1 2+α π 1 1−2α π.(24) As the Schwarz-Christoffel transformation is conformal, and θ(r) is harmonic, we are going to find the solution to the Laplace equation in the half ζplane, ˜ θ, subject to the boundary conditions ˜ θ=α∞for |x|>1+1and ˜ θ= α∞+(x/|x|)(α−α∞)for2<|x|<1−1. The solution ˜ θ(x,y) in the image of the region I on the ζplane, I,isgiven by Ref. [14] ˜ θ(x,y)=α∞+α−α∞ πarctan y x−1 −2(α−α∞) πarctan y x+α−α∞ πarctan y x+1. (25) From the solution ˜ θ, we can obtain θ(x,y)= ˜ θ[x(x/L,y/L;α),y(x/L,y/L;α)]. We note that we do not have an explicit expression for xand yas functions of x and y, so we cannot give an analytic expression for θ(x,y). However, we can evaluate exactly fIsince fI=K 2λI dxdy|∇θ|2=K 2λI dxdy|∇˜ θ|2 =K 2λB ˜ θ(ν·∇˜ θ)ds,(26) where Bis the image in the ζplane of the boundary of zone I. After a straightforward calculation, we get the expressions for fIand BI(α), fI=qK(α) 2π−ln qRccos α π−ln 3 2−α π1 2+α π −1 2−α πln π 2+α π 2−α (27) BI(α)=−K(α)1 2−α πln π 2+α π 2−α +ln 3 2−α π1 2+α π.(28) Figure 5shows an excellent agreement between the theoretical prediction for BI(α) from Eq. (28) and the results obtained in the previous subsections, except close to α=π/2. The latter may be due to numerical uncertainties in the evaluation of BI,ns, since either it diverges (N⊥texture) or vanishes (N texture) in that limit. As happened for BI,s, the dependence on the texture comes from K(α), leading to the same conclusions about the relative stability of the nematic textures with the angle α, in agreement with the results obtained in Ref. [14]. 0π/8 π/4 3π/8 π/2 α(rad) -0.4 -0.2 0 BI/K(α) Exact N⊥ (BEM) N|| (BEM) 0π/8 π/4 3π/8 π/2 α(rad) 0.001 0.01 0.1 1 10 -BI/K FIG. 5. Plot of BI/K(α) as a function of α. The continuous broad line corresponds to the exact expression, and the symbols correspond to the estimates from the numerical data for fI,ns (from the boundaryelement method): circles for the N⊥texture and diamonds for the N texture. (Inset) Plot of −BI/K as a function of α. The meaning of the symbols is the same as in the main panel. IV. EVALUATION OF fII In Ref. [27] it was observed that there is a dependence on w of the next-to-leading contribution to the surface free-energy density. However, in the previous section we have shown that the contribution from region I depends only on α.Sowe anticipate that this dependence comes from region II, where inhomogeneities of both the nematic order parameter Sand orientational θfields are observed. The free energy of region II can be evaluated as the sum of the contributions of the regions around each cusp (either wedge or apex). In each of these regions, and if RcL,we anticipate that the orientational field far from the cusp behaves asymptotically as that of the isolated disclination line which has been nucleated at the substrate wedge or apex: θ∼Iφ, where Iis the topological charge of the disclination line and (r,φ) are the polar coordinates taking as origin the cusp. The values of Iare fixed by the (strong) anchoring conditions on the substrate and depend on the nematic texture: I1= −α/(π/2−α) (respectively, I1=+1) at the wedge bottom and I2=α/(π/2+α) [respectively, −(π/2−α)/(π/2+α)] at the apexes for the N⊥(respectively, N)[14,27]. However, Sdecreases as r→0 from the bulk value far from the origin. Although we cannot solve analytically this problem, we can estimate, using an ansatz, the free energy associated to this combined distortion of θand Sassuming a melted core (i.e., S→0asr→0), as shown in the Appendix. However, we resort here to a full minimization of the LdG model [Eq. (5)] restricted to region II. At r=Rc, we impose Dirichlet boundary conditions to Q, where Stakes the value corresponding to the order-parameter profile corresponding to a planar wall at the distance between the boundary point and the closest substrate, the biaxiality parameter B=0 and θ=α∞+I(φ−π/2). Alternatively, we used free boundary conditions ν·∇Q=0, leading to similar results. Figure 6 shows typical textures obtained after minimization. These are similar to those obtained from full minimization close 041706-6 GENERALIZED BERREMAN’s MODEL OF THE ELASTIC ... PHYSICAL REVIEW E 86, 041706 (2012) FIG. 6. (Color online) Typical textures obtained from the full minimization of the LdG function in region II for τ=1, κ=2, w=1.0, L=90, and α=π/6 (upper panels) and α=π/3 (lower panels). The color map corresponds to the reduced nematic order parameter field ˜ S=S/Sband the segments correspond to the nematic director field n. to the cusps. At distances rlarger than a few correlation lengths, Sdecays to the bulk value except in the neighborhood of the substrates, where it takes approximately the value corresponding to the flat substrate profile. In fact, the nematic order parameter profile in general shows only slight variations close to the cusps. This situation would correspond to a surface disorientation wall, in contrast with true disclination lines which exhibit melted cores. An exception to this general situation occurs when α≈π/4, where a complex structure with a true defect line bound to the cusp is observed. On the other hand, the orientational field deforms continuously in order to satisfy anchoring conditions on each side of the wedge or the apex. So we anticipate that the contribution to the surface free-energy density fII will scale with Rc as λfII =4σnwRc+K(α)lnRc+BII.(29) The first contribution arises from the inhomogeneities of S close to the substrates and the second one from the asymptotic behavior of θfor r>1: By using the Frank-Oseen functional Eq. (11), and taking into account that ∇θ=Iuφ/r (where uφis the azimuthal unit vector), we find that the free-energy contribution is proportional to (K/2)I2φ ln Rc, where φ is the opening angle of the wedge or apex. From this result, we get the second contribution in Eq. (29) by noting that K(α)=KI2 1(π−2α)/2+KI2 2(π+2α)/2. The remaining contribution will give the core free energy per unit length associated to the disclination lines. Expression (29) can be used to extract BII from the numerical minimization. So fixing the value of α, and for each w(between 0 and 1.5), we considered a range of values of Rcbetween 10 and 90. The minimization was performed by 10 100 RC 4 5 6 7 8 9 λfII-4σNWRC 20 30 40 50 60 70 80 90 RC 0 1 2 3 BII FIG. 7. (Left) Representation of λfII −4σNWRcas a function of Rcfor w=1andα=π/6 (squares) and α=π/3(circles). The straight lines correspond to the logarithmic regressions of the numerical data. (Right) Representation of BII(α,w) as a function of Rc. The meaning of the symbols is the same as in the left panel. the mesh-adaptive finite-element method used previously [37], taking τ=1 and κ=2 (this choice is motivated to compare with results reported in the literature [27,29]). After substrating 4σnwRcto λfII, the numerical results clearly show a logarithmic dependence on Rc, with a slope approximately equal to K(α) (see left panel in Fig. 7). The next step was to subtract the subdominant contribution K(α)lnR∞. Now the numerical results are nearly independent of Rc(see the right panel in Fig. 7). The value of BII is estimated as the mean value of these results, with an error bar given by the dispersion of the numerical data around the average. V. DISCUSSION The results obtained in the previous sections can be combined as follows: f=fI+fII +fIII =σNW cos α+felastic,(30) where felastic =q 2πK(α)−ln qcos α π−1 2−α πln π 2+α π 2−α −ln 3 2−α π1 2+α π+q 2πBII(α,w) (31) up to corrections of order of q2. Consequently, B(α,w)=BI(α)+BII(α,w) =−K(α)1 2−α πln π 2+α π 2−α +ln 3 2−α π1 2+α π+BII(α,w). (32) This is the main result of our paper. We now can check this prediction by comparing these results with those obtained from the full minimization of the LdG model of a nematic in contact with a sawtoothed substrate [27]. In this reference, we considered NI coexistence conditions, κ=2, and a range of values of the fluid-substrate coupling 0.1<w<1.5. Consequently, we are considering substrates with extrapolation lengths between ξand 10ξ. These values are smaller than the experimentally relevant situations even for strong anchoring 041706-7 O. A. ROJAS-G ´ OMEZ AND J. M. ROMERO-ENRIQUE PHYSICAL REVIEW E 86, 041706 (2012) 00.5 11.5 w -3 -2 -1 0 1 B(α,w) L*=16 α=π/6 L*=16 α=π/3 L*=32 α=π/6 L*=32 α=π/3 L*=48 α=π/6 L*=48 α=π/3 L*=64 α=π/6 L*=64 α=π/3 L*=96 α=π/6 L*=96 α=π/3 This work α=π/6 This work α=π/3 FIG. 8. (Color online) Comparison between the values of B(α,w) obtained from the full minimization of the LdG model [27] and the results obtained in the present work. Symbols correspond to B(α,w) obtained from the full minimization of the LdG model for different values of Land α=π/6intheN⊥texture (black symbols) and α=π/3intheNtexture (orange or light gray symbols). The wide blue (black) and wide red (dark gray) lines with error bars correspond to the predictions from this work for α=π/6andπ/3, respectively. conditions (although under some circumstances they may be achieved [41]). On the other hand, the range of values of L, between 16 and 96, is also smaller than the experimentally feasible range. This range was dictated by computational limitations, as problems with the numerical procedure were observed for larger values of L. Our interest in the considered range of values of wis twofold. First, this range covered the typical values of wfor which wetting and related phenomena on rough substrates were observed in our model [29–31,35]. Second, the collapse of the data for B(α,w) obtained for the considered values of Lwas achieved only for w>0.5. However, we expect that this collapse will be observed also for smaller values of wif we increase the range of L, and the trends shown in the available data support this claim. Figure 8shows the comparison between the results reported in Ref. [27](τ=1, κ=2) and the calculated ones in this paper for α=π/6(N⊥texture) and α=π/3(Ntexture). For large w, the agreement is good, although our results slightly overestimate those from the full minimization. On the other hand, for small wthe curves obtained from the full minimization converge towards our prediction as Lincreases. So, our approximation is accurate even for moderate values of L, despite the assumptions involved in our approach. Consequently, the scheme considered in this paper is an alternative to the full minimization technique, which is quite expensive from a computational point of view, when wL > 1. Finally, it is interesting to note that the elastic contribution to the surface free-energy density, Eq. (31), can be expressed as follows: λfelastic =22πKI1 2 I2 2ln γ(α) L+BII(w,α),(33) where we note that K(α)=−4πK(I1/2)(I2/2), with I1and I2being the topological charges associated to the disclination lines at the wedges and apexes, respectively. The characteristic length γ(α), which absorbs the contribution to Bfrom region I, is defined as γ(α)=3 2−α π1 2+α ππ 2+α π 2−α1 2−α π .(34) This length scale decays continuously from π/2forα→0 to1forα→π/2. We can understand the first contribution in Eq. (33) as the interaction between the disclination line at x=0 with half its topological charge and the apexes at x= ±Lcos α, again with half their topological charges [42]. As the absolute values of the topological charges are always smaller or equal to 1, the nucleation of disclinations at the cusps of the surface is favorable with respect to bulk disclinations (which only may have half-integer values). With this intepretation, γ(α) can be used to define an effective core radius γ(α), and the sum of the core energies is given by BII(α,w). VI. CONCLUSIONS In this paper we have analyzed the elastic contribution to the surface free-energy density for a nematic in contact to sawtoothed substrates in the strong anchoring regime, i.e., wL > 1. We have extended the analysis done in Ref. [27]. So, in addition to the leading contribution proportional to −qln q, with qbeing the substrate periodicity wave number q, we have characterized the next-to-leading term. This term has two contributions: one associated to the deviation of the orientational field with respect to the contribution of the array of disclination lines nucleated at the cusps of the substrate and the core free-energy associated to them. We anticipate that our analysis can be generalized for other substrate shapes when the nematic texture presents topological defects induced by the structure, as, for example, in crenellated substrates [35]. Furthermore, our results may be used to predict accurately the location of first-order wetting transitions in nematic liquid 00.5 11.5 w -3 -2 -1 0 1 2 3 4 B II(α,w) N⊥, α=π/6 (Ansatz) N|| , α=π/3 (Ansatz) N⊥ , α=π/6 (full minimization) N|| , α=π/3 (full minimization) FIG. 9. Comparison between the values of BII(α,w) obtained from the full minimization of the LdG model and the ansatz (see text). Symbols correspond to BII(α,w) obtained from the full minimization of the LdG model for α=π/6intheN⊥texture (circles) and α=π/3intheNtexture (up triangles). The wide continuous and dashed lines correspond to the estimates from the ansatz for α=π/6 and α=π/3, respectively. 041706-8 GENERALIZED BERREMAN’s MODEL OF THE ELASTIC ... PHYSICAL REVIEW E 86, 041706 (2012) crystals in contact to general substrates [29,30,35]. However, we must note that our analysis is restricted to bulk nematics liquid crystals with in-plane deformations in the presence of grooved substrates. The effect of the substrate structure in the nematic texture in partially filled configurations (i.e., with an isotropic fluid in bulk), the influence of twist nematic deformations [19], and the effect of the structure in the z direction deserve further study, but this is beyond the scope of this work. ACKNOWLEDGMENTS The authors thank P. Patr´ ıcio and N. M. Silvestre for very stimulating conversations and technical advice on the finiteelement method and M. M. Telo da Gama, L. F. Rull, and A. O. Parry for enlightening discussions. We acknowledge the support from MICINN (Spain) through Grant No. FIS2009- 09326 and Junta de Andaluc´ ıa (Spain) through Grant No. P09- FQM-4938, both cofunded by EU FEDER. APPENDIX: ANSATZ FOR THE EVALUATION OF fII In order to estimate the free energy associated to region II, we can make use of an ansatz. As argued in the paper, the nematic order parameter Smust vanish at the cusps of the substrate to avoid the divergence of the free energy, but it must converge to the bulk value far from the wedges and apexes. So, in the region around each cusp, we can suppose that Sdepends only on the radial distance rfrom the wedge or apex. In particular, we use the following ansatz [43]: S(r)=Sb(1 −e−r β), where Sbis the bulk nematic order parameter and βis a length scale to be determined later. We will assume that the biaxiality parameter B=0 everywhere. Regarding the orientational field θ, we assume that has the expression corresponding to a disclination line placed at the origin, θ=α∞+I(φ−π/2), where Iis the winding number associated to the disclination line and α∞is the far-field value of the orientational field (0 for the N⊥texture and π/2for the Ntexture). Substituting this ansatz in the LdG functional [Eq. (5)], this expression reduces to a function of β, λf ≈13 144(π∓2α)β2−2w(Rc−β)+1 79I2(π∓2α)γ+ln Rc 2β+11 16(π∓2α) +1 8sin[(I−1)(π∓2α)] (I−1) ∓18Isin(2α)−3[sin(−2I(π∓α)+(I−1)π)−sin(−2I(±α)+(I−1)π)],(A1) where we have considered the case τ=1 and κ=2 (for other values, we can get straightforwardly analogous expressions). In this expression, the upper sign corresponds to the wedge situation (so I=I1), and the lower one to the apex (so I=I2), and γis the Euler constant. Minimizing Eq. (A1) with respect to β, we obtain the optimal value for this length scale: β=−72 13 w (π∓2α)+72 13 w (π∓2α)2 +648 91 I2. (A2) Substituting Eq. (A2) into Eq. (A1), we get the freeenergy estimate associated to each region around a cusp, which provides an upper limit to the real value of λf . In order to have an estimate of BII, we subtract to the obtained values of λf through this ansatz the surface and elastic contributions −2wRcand (K/2)I2(π∓2α)lnRc= (9/7)I2(π∓2α)lnRc, respectively. After this, the estimate of BII is obtained as the sum of the results obtained for the wedge and apex. Figure 9shows the comparison of this estimates with the results obtained from the full minimization for α=π/6 (N⊥texture) and α=π/3(Ntexture). They have the same order of magnitude, although our ansatz overestimates the value of BII, and, qualitatively, the dependence of BII with wis recovered. [1] B.-W. Lee and N. A. Clark, Science 291, 2576 (2001). [2] J.-H. Kim, M. Yoneya, and H. Yokoyama, Nature 420, 19 (2002). [3] S.Ferjani,Y.Choi,J.Pendery,R.G.Petschek,andC.Rosenblatt, Phys. Rev. Lett. 104, 257801 (2010). [4] C. V. Brown, M. J. Towler, V. C. Hui, and G. P. Bryan-Brown, Liq. Cryst. 27, 233 (2000). [5] C. Uche, S. J. Elston, and L. A. Parry-Jones, J. Phys. D: Appl. Phys. 38, 2283 (2005). [6] C. Uche, S. J. Elston, and L. A. Parry-Jones, Liq. Cryst. 33, 697 (2006). [7] A. J. Davidson, C. V. Brown, N. J. Mottram, S. Ladak, and C. R. Evans, Phys.Rev.E81, 051712 (2010). [8]C.R.Evans,A.J.Davidson,C.V.Brown,andN.J.Mottram, J. Phys. D: Appl. Phys. 43, 495105 (2010). [9] O. J. Dammone, I. Zacharoudiou, R. P. A. Dullens, J. M. Yeomans, M. P. Lettinga, and D. G. A. L. Aarts, Phys. Rev. Lett. 109, 108303 (2012). [10] N. M. Silvestre, P. Patr´ ıcio and M. M. Telo da Gama, Phys. Rev. E69, 061402 (2004). [11] T. Ohzono and J.-i. Fukuda, Nat. Commun. 3, 701 (2012). [12] D. W. Berreman, Phys.Rev.Lett.28, 1683 (1972). [13] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, 2nd ed. (Oxford University Press, Oxford, 1995). [14] G. Barbero, Lett. Nuovo Cimento Soc. Ital. Fis. 29, 553 (1980). 041706-9