Computing with non-orientable defects: Nematics, Smectics and Natural patterns

Defects, a ubiquitous feature of ordered media, have certain universal features, independent of the underlying physical system, reflecting their topological, as opposed to energetic properties. We exploit this universality, in conjunction with smoothing defects by"spreading them out,"to develop a modeling framework and associated numerical methods that are applicable to computing energy driven behaviors of defects across the amorphous-soft-crystalline materials spectrum. Motivated by ideas for dealing with elastic-plastic solids with line defects, our methods can handle order parameters that have a head-tail symmetry, i.e. director fields, in systems with a continuous translation symmetry, as in nematic liquid crystals, and in systems where the translation symmetry is broken, as in smectics and convection patterns. We illustrate our methods with explicit computations.


Introduction
Macroscopic physical systems consist of large numbers of interacting (microscopic) parts and are described by thermodynamic principles. A central tenet of thermodynamics is that the equilibrium state, and the relaxation to equilibrium, are described by an appropriate free energy [Gre95]. While the details differ, free energies describing systems that spontaneously generate ordered/patterned states have certain universal features independent of the underlying physics. These include 1. nonconvexity of the free energy and the existence of multiple ground states for the system.
2. regularization by a singular perturbation (an "ultraviolet cutoff") to preclude the formation of structures on arbitrarily fine scales.
These features are present in free energies that describe many systems including liquid crystals [Vir95], micro-magnetic devices [DKMO00] and solid-solid phase transitions [KM94]. This collection of ground states, which we denote by V, is given by the cosets G/H where G is the symmetry group of the system (e.g. the Euclidean group for systems without a preferred origin and orientation) and H is the subgroup of symmetries that are unbroken in an equilibrium state [TK76,Kle95].
If H is a proper subgroup of G, i.e. if some symmetries are indeed broken in the equilibrium states, V is nontrivial. If V is nontrivial, the system naturally forms defects, by which, we mean the occurrence of localized regions where the microscopic degrees of freedom have not (and perhaps cannot) relax to an equilibrium state in V.
Why should defects occur in such systems? For example, in the case of a large Prandtl number convection, fluid is heated uniformly from below, whereas the emerging flow pattern self-organizes into rolls/stripes with a preferred wavelength reflecting the breaking of the continuous translation symmetry of the ambient space/forcing through heating by a discrete symmetry of translation by one wavelength perpendicular to the stripes. The rotational symmetry of the system is however unbroken so that there is no preferred orientation of the stripes. In large aspect ratio systems where the box size is very large compared to the preferred wavelength, the local orientation is chosen by local biases (such as boundary effects). Thus, the emerging pattern is a mosaic of patches of stripes with preferred wavelength with different orientations, patches which meet and meld along (in 2D) grain boundaries which themselves meet at points. These lines and points in 2D (planes, loops and points in 3D) constitute the defects in the convection pattern.
This argument is not peculiar to convection patterns. Defects arise as well in crystalline elastic solids, and complex fluids in the nematic and smectic phases, including dislocations, disclinations, grain boundaries, and twin (phase) boundaries. Fundamentally, defects arise in these phases due to the presence of microscopic structural symmetries. In crystalline solids the possibility of nontrivial deformations that preserve lattice periodicity locally gives rise to defects, and in nematic liquid crystal phases it is due to the head-tail symmetry of the director field.
The macroscopic state of an extended system is therefore best understood as a patchwork of domains that meet at various types of defects [Mer79]: disclinations, dislocations, monopoles, walls, etc. It is thus of interest to develop tools that allow us to understand, predict, control and manipulate energy driven pattern formation. A natural approach to analyzing these systems is thus 'two-pronged', namely describing the system at the level of domains corresponding to ground state configurations from V, and patching together multiple such domains across defects. Within each domain we therefore seek to obtain effective pdes for the macroscopic order parameters, i.e., the element of V, that describe the local structure. For energy driven systems, these equations arise from averaging the free energy over all the microstructures consistent with the macroscopic order. Since these order parameter equations depend (largely) only on the relevant broken symmetries, they are universal, i.e., the same equations arise in a variety of physical contexts [PN94]. The topologies of allowed defects are also universal [KMT77,Kle95], and they are captured by discontinuous and/or singular 'solutions' of the averaged pdes that are initially derived to describe the domains [EINP00]. This universality is one motivation for the work presented in this paper, namely, the idea that there is a common modeling and computational methodology which can be applied to study defects in systems with vastly different physics at widely separated scales. Because of the nature of the singularities involved, often involving serious non-uniqueness of 'solutions,' practicality forces slightly more refined considerations resulting in the enhancement of the basic model with defect fields and equations, nevertheless of a universal nature; in this paper we develop such models.
In elastic solids, there is a historically systematic way of interpreting the geometric and some elastic aspects of dislocation and disclination defects beginning from Volterra [Vol07], that has been recently generalized and incorporated into modern continuum thermomechanics, generating models for practical application. In this viewpoint, a disclination is interpreted as a terminating curve of a surface on which the elastic rotation is discontinuous -this idea has been generalized to consider terminations of surfaces of elastic distortion discontinuity (including strain), with the resulting defect called a generalized disclination [ZA18,ZAP18]; a dislocation is understood as the terminating curve of a surface on which the (inverse) deformation is discontinuous. The core idea lends itself to useful generalization, as has been done to define higher-order branch point defects [AV19]. By these physical definitions it becomes clear that a dislocation induces a lower energy state in a body than a single disclination, as shown in Fig. 1. It can also be shown, both mathematically and physically [ZA18], that a disclination dipole pair in a solid may be interpreted, for small inter-disclination spacing, as a dislocation. It is no surprise then that dislocations (in the bulk) and disclination dipoles (at interfaces) are the defects that are most often observed in solids but, as as seen in experiment, individual disclinations in isolation are found in nature at the junction of twin and grain boundaries. The relevant field variable for elasticity, in the broader setting including defects, is the inverse deformation w t : Ω t → R 3 , a vector field on the deforming body Ω t at time t, and the natural order parameter is the inverse deformation gradient Dw t (a 3 × 3 matrix). As we discussed above, defects in elastic materials are encoded by discontinuities of w t and Dw t . The theory of elasticity and of plasticity and other defect mediated behaviors in solids results from defining an appropriate dynamics for the elastic and plastic parts of Dw t , which in essence correspond to the absolutely continuous part and physically mollified, singular part of Dw t , respectively, along with the motion of the body in ambient Euclidean space 1 . We will use this framework as a template to develop models that adequately describe the gross continuum mechanical features of nematic and smectic liquid crystals. These models can describe the slow evolution of defects 2 and are thus also useful in studying the long term coarsening dynamics of natural stripe patterns.
A nematic liquid crystal consists of rod-like molecules and the relevant microscopic field is the (distribution of) orientation of these molecules. The natural order parameter is a director field n(x) that encodes the long-range ordering of the nematic. A smectic is defined by the nematic director n(x) and an additional (nearly periodic) field ρ(x) that describes the maxima and minima of the density of the liquid crystal molecules. While n and ρ are independent in principle, for smectic-A liquid crystals, n is oriented along ∇ρ motivating de Gennes' definition of a complex order parameter ψ = Ae iθ with ρ ∼ (ψ). For stripe patterns, the microscopic field is a locally periodic function u = u(θ(x, y)) = Ae iθ + A * e −iθ . Near defects, A and θ are independent so the relevant order parameter is Ae iθ . Away from defects, however, A is slaved to θ, A = A(|Dθ|) so the director field k = ±Dθ is the appropriate order parameter k.
In all these cases, the order parameter is a director field n or k that should be compared with the elastic (or absolutely continuous) part, W , of the inverse deformation gradient Dw t for solids. We can make this analogy if we consider a scalar field θ : Ω t → R (instead of a vector field), so that Dθ, the analog of Dw t , will give a gradient/vector field whose regular/absolutely continuous part is a representative of the director field n or k. At the same time, the dominant elasticity in solids is related to Dw t whereas that of nematics is related to Dk. It is this thread of conceptual unification, while being fully cognizant of the essential physical differences between different systems, that we pursue in this paper, with the formulation of mathematical models that exploit this analogy and the determination of approximate solutions of this model.
What is the motivation for developing a new modeling framework for nematics, smectics and natural patterns? There is, of course, the state-of-the-art model for the mechanics of liquid crystals -the Landau-DeGennes (LDG) model -that can predict a variety of defect related phenomena in liquid crystalline phases [dGP95,SV12]. However, it is fair to say that the LDG model sheds no light on any possible connection that might exist between defects in elasticity of solids and those in liquid crystals, even though such connections have been known to exist from almost the inception of the modern theory of liquid crystals starting with Frank's seminal paper [Fra58] and the fundamental one-to-one correspondence between the elastic fields of the screw dislocation in an elastic solid and the wedge disclination in a nematic liquid crystal. Our paper aims to establish this connection and, in this sense, may be thought to be in the spirit of the work of Kleman [Klé73,KF08]. Abstracting some of the mathematical ideas of defects (and incompatibility) in elastic solids [ZAWB15], a first successful attempt at representing the energetics and dynamics of a planar director field in nematics along these lines has been demonstrated in [ZZA + 16], which is capable of representing higher-strength defects than ± 1 2 as fundamental entities, including their annihilation and dissociation. However, this model relies strongly on the angle parametrization of a planar director field and does not generalize cleanly to 3D, i.e. without any reliance on a specific paramaterization of the director field, as shown in Appendix A. The full 3-d model we present in this paper overcomes this shortcoming, while being physically more in line with elastic defect theory.
Stripe patterns form in extended systems when a homogeneous state loses stability to periodic state(s) with a preferred wavelength, but no preferred orientation [SH77]. As we argued above, in large systems whose size is much bigger than the preferred length scale, defects are both inevitable and ubiquitous. The Cross-Newell equation [CN84], and its weak/singular solutions [NPB + 96], describe the long term evolution of stripe patterns and the defects in them. The parallels between defects in stripe patterns and in elastic sheets have been explored in recent work [NV17]. In this paper we establish a further parallel with the elastic theory of defects in solids, and treat smectic liquid crystals and stripe patterns within a unified framework. Although the defects we consider have universal features, it is important to recognize that defects have important differences, as well, in different physical systems. Therefore, our goal is to develop a framework that can be adapted to the particular physics of various systems, and thus has wide applicability.
We conclude this introduction with a guide for the reader. In earlier work by the first two authors, we have discussed the relationship between defects in solids and those in liquid crystals [ZAWB15, ZZA + 16], and used these ideas to develop numerical approaches to computing the be-havior of defects in these systems. We expand this circle of ideas in Sec. 2 where we review, and provide a fresh perspective on the mathematical theory of stripe patterns. We identify the natural topological defects in this theory and discuss the adequacy of this simplified framework for studying the energetics of patterns, but also of nematic and smectic-A liquid crystals. The solutions of this theory naturally possess non-integrable singularities in the energy density which motivates our approach, outlined in Sec. 3, to model defects in scalar fields from the perspective of defects in elastic solids leading to plasticity and phase transitions. A primary feature of such models is the replacement of the singular parts of solutions to 'classical' defect models by new independent fields (bearing much similarity to gauge fields) that are nevertheless smooth (but localized), and the utilization of integrable energy densities. The model is constrained by thermodynamics, the balance laws of mechanics, and statements reflecting conservation of topological charge of defects, and it has a 'good' numerical formulation in terms of smoothed defect fields. In Sec. 4 we demonstrate the scope of our model through numerical computations of equilibrium features of various defect configurations in nematics and in smectics/patterns. Section 5 presents an overview of our work along with a concluding discussion of its implications and possible future extensions.
2 Defects in a scalar field and analysis of natural stripe patterns A microscopic model for the evolution of convection patterns in large Prandtl number fluids is the Swift-Hohenberg equation [SH77], where R > 0 is the forcing parameter (analog of the Rayleigh number for convection), u is a proxy for the vertically averaged temperature and k 0 is the preferred wavenumber of the roll patterns. Although convection patterns are two dimensional, the Swift-Hohenberg equation is defined in an arbitrary number of dimensions. For our purposes, we consider domains Ω ⊂ R d , d = 2, 3, and let x denote an arbitrary point in Ω. The L 2 gradient flow u t = − δ δu E for the SH energy is the Swift-Hohenberg equation (1) with natural boundary condtions. E[u(x, t)] is monotonically non-increasing along solutions. The asymptotic states are critical points for the non-convex energy (2). Ignoring the boundaries, the ground states for this energy form a manifold V that consists of periodic stripe patterns with a range of allowed wave numbers in the vicinity of k 0 and all possible orientations [CE90]. Cross and Newell [CN84] demonstrate that, away from defects, stripe patterns are modulations of these stable critical points, u(x, t) = w 0 (θ(x, t), |k| 2 ), where k = Dθ varies slowly in space and time. w 0 (q · x + θ 0 , |q| 2 ) is 2π periodic in its first argument, and is "normalized" so that it's maxima (resp. minima) occur at q · x + θ 0 = 2nπ (resp. (2n + 1)π). w 0 describes the profiles of stable stripe patterns in V with a constant wave-vector q. Away from defects, Dk O( ) 1 where is the ratio of the pattern wavelength to the size of the system. The wave vector k is the appropriate macroscopic (i.e. slowly varying) order parameter, and a given k is consistent with distinct "microstates" corresponding to different choices of θ 0 , the constant of integration needed to recover the phase from the equation Dθ = k. This translation invariance θ(x, y) → θ(x, y) + θ 0 implies that the linearization of the dynamics in (1) about a modulated stripe pattern has a non-trivial kernel, and the corresponding solvability condition yield, at lowest order, the (unregularized) Cross-Newell equations [CN84], a gradient flow, that describes the macroscopic dynamics for the wave-vector field k(x, t). These equations lead to the formation of shocks, so they need to be regularized by higher order effects in the small parameter [NPB + 96, EINP00]. An alternative to employing the Fredholm alternative/solvability is to directly compute an effective energy E[k(x, t)] by averaging the energy (2) over all the microstates that are consistent with a given macroscopic field k(x, t) [NV17]. This is equivalent to averaging over the phase shift θ 0 ∈ [0, 2π], and yields the Regularized Cross-Newell (RCN) energy where W is a nonconvex "well potential" in k. For our purposes, W (|k| 2 ) = (|k| 2 −1) 2 is an adequate approximation. In the context of patterns, we also need that k ⊥ should yield a measured foliation [Poé81], corresponding to a stripe pattern with constant width. Consequently, a necessary condition for k to describe a smooth stripe pattern is that curl k = 0 [FLP12]. With the substitution k = Dθ (equivalent to the constraint curl k = 0), we recognize that the RCN energy is closely related to the Aviles-Giga functional [AG87], which was introduced as a model for smectic liquid crystals.
On the other hand, we can also consider the energy (3) without the constraint curl k = 0. In this setting |k| → 1 in the long time dynamics without k necessarily being a gradient, and the RCN energy is related to the Ginzburg-Landau functional. The Ginzburg-Landau functional is defined for a complex order parameter ψ = Ae iθ , and away from defects, A ≈ 1, so the connection with (3) is through setting k = e iθ , a unit vector. (1) starting from small random fluctuations of the homogeneous state u = 0 ( Fig. 2(d)) and the (numerically suggested but mathematically unproven) global minimum of the Swift-Hohenberg energy (Fig. 2(c)).
Note the significant differences in the patterns along the major axis of the elliptical domain for the AG and SH minimizers in Fig. 2. These differences are entirely due to fact that the local phase of a pattern is "multi-valued" [Kle08]. The phase θ is not directly observable, unlike the pattern field u(x, t) in (1). Thus, we have to identify different phase functions θ that give the same field u = w 0 (θ, |k| 2 ) where k = Dθ. Since w 0 is an even, 2π periodic function of the first argument, we have the identifications θ → θ + 2nπ, k → k where n is any integer (periodicity or symmetry under translations by multiples of 2π) and θ → −θ, k → −k (evenness or head-tail symmetry). In particular, the values θ = mπ with integer values for m are distinguished, since we can apply a combination of the two symmetries to achieve θ → −θ → −θ + 2mπ = θ, k → −k → −k, so the contours θ = mπ are the locations which can support disclinations, i.e. flips θ → θ, k → −k as illustrated in Figs. 3(a) and 3(b). Disclinations are thus point defects with nontrivial monodromy for the map x → k = Dθ, so they are necessarily codimension 2, i.e. points in 2D and lines in 3D. These arguments were used to develop a variational theory for disclinations in stripe patterns [EV09]. These conclusions were also obtained independently in [CAK09, PSS14, AMK17] for smectic liquid crystals using different arguments.
It is an interesting question as to how one computes gradient flows for functionals that depend on such "multi-valued" fields θ and (the director) k. One approach is to introduce branch cuts to obtain a single valued phase θ and a vector field k = Dθ on a branched double cover [LZ04, Chap. 1] of Ω where the ramification points are the locations of the convex and concave disclinations that carry the nontrivial monodromy of the map x → k(x) [EINP03].  Let γ : [0, 1] → Ω be a simple, closed, smooth, oriented curve in Ω, i.e. γ(0) = γ(1). Ifγ is a continuous lifting of γ to the double cover, it is no longer necessarily closed, because γ(0) and γ(1) might end up in different sheets of the branched covering. We can always assume that γ intersects the branch cut set S transversally at finitely many points, x 0 , x 1 , x 2 , . . . , x n−1 , x n = x 0 , indexed consecutively with the same orientation as γ. We denote the segment ofγ between x i and x i+1 byγ i and the oriented jump of a quantity q at denotes the limiting values of q as we approach x i from either side along the curveγ. The integral of the distributional gradient of the single valued lifting of θ around a closed curveγ is zero. Consequently, ifγ(0) =γ(1) (they are on the same sheet), the curve γ in Ω must intersects the branch cut set S an even number of times. We can now decompose the distributional gradient of Dθ into its regular part k and its singular (jump) part to obtain [θ] γ is the negative of the sum of the (oriented) jumps in θ at the points where γ intersects S. Since the branch cut set S consists of curves where θ a multiple of π when approached from Ω \ S, it follows that [θ] γ = mπ, where m is an integer that depends on γ. θ is thus the analog, in layered/foliated media, of the deformation in a crystalline solid and [θ] γ is the "Burgers scalar" [KM16,AMK17]. Note that, the argument makes it clear that [θ] γ is only defined if the total number of flips within γ is even.
We also note that k can jump at the intersections of γ with S, and the net jump is given through A, the absolutely continuous part of the distributional gradient Dk, by γ A · dx = [k] γ . In analogy with the simplest, and most practical, definitions of elastic defects (cf. [ZA18,§1]), the quantity [k] γ is non-zero if the curve γ encloses a net disclination density. In what follows, we adopt the viewpoint that arbitrary (composite) defects and defect distributions can be decomposed into sums of elementary defects, corresponding to convex and concave disclinations. This approach has  The dislocation and disclination densities given by [θ] γ and [k] γ are "linear" defect measures, in that they depend linearly on (jumps in the fields) θ and k = Dθ outside S. For patterns, and in 2-d, there is a further defect measure J = Det(D 2 θ) = θ xx θ yy − θ 2 xy , the Gaussian curvature of the phase surface (x, y) → θ(x, y) [NPB + 96, NV17]. If k = 1 a.e., J is related to the field curl λ in the angle parametrization (cf. Appendix A) which was introduced in [ZZA + 16] as a representation of defects in planar nematic director fields. For a domain D whose boundary γ = ∂D intersect S transversally, relates the total twist of the angle parameterization along γ to the mass of T of the measure J, which equals the area enclosed by the image of γ, normalized by π, under the mapping Dθ : Ω → R 2 . J is one of an entire family of such measures that can encode this twist. For any (nonlinear) scalar valued function f (k), that has a branch cut intersecting the unit circle |k| = 1 transversally, integrating the absolutely continuous (regular) part of Df •k along a closed curve in the real domain picks up a jump of f • k when Df • k has a singularity in the domain inside the curve. The function f (k) = tan −1 ( k·e 2 k·e 1 ), where (e 1 , e 2 ) is an arbitrary orthonormal frame, gives the defect measure J. This specific function f , the angle parametrization of the planar director field, was used as the fundamental field in the model in [ZZA + 16]. In the current work, we instead use the entire vector k, because the resulting theory generalizes easily to 3 dimensions, and has closer connections to defects in solids, as we discuss in Sec. 3. For completeness, we present an extension of the angle parameterization to 3 dimensions in Appendix A.
As a final remark, we observe that the regularized Cross-Newell energy (3) is defined in terms of k, the regular part of Dθ, and away from defects, we have the constraint curl k = 0, like for the Aviles-Giga energy functional. However, if we allow for defects, then, we have k · dx = mπ is no longer necessarily zero, so that the distributional curl of k can be non-zero, albeit quantized. Natural patterns and smectic liquid crystals are therefore in a seperate universality class that lies in between the Aviles-Giga and the Ginzburg-Landau energy functionals. in Ω such that Ω\S is simply connected. Given a smooth, symmetric second-order tensor fieldÃ with vanishing curl on Ω, i.e., satisfyingÃ ij =Ã ji ,Ã ij,k =Ã ik,j , the question is to characterize the jump of the fieldθ across S if the second derivative ofθ equalsÃ on Ω\S. Thus, we seek solutions to the system Dθ =k and we are interested in evaluating where x ± are any sequences of points that approach x ∈ S from the ± sides of S at x, respectively. For any closed loop l in Ω that cannot be continuously shrunk to a point while staying within Ω, is a constant sinceÃ is continuous and curl-free in Ω. We define the constant t to be the disclinationstrength of the fieldÃ in the domain Ω.
SinceÃ is curl-free in Ω\S, the latter being simply connected, it is possible to construct a field k on the same domain satisfying (5) 2 . By (6),k, in general, has a non-vanishing jump across S.
For a fixed S and its correspondingk, it is further possible to construct a scalar fieldθ satisfying (5) 1 on Ω\S due to the symmetryÃ ij =Ã ji .
Let c be a curve belonging to S joining points x, y ∈ S, and c ± two curves near c on the ± sides of S. Noting that t = Dθ = k is a constant on S, computing c ± Dθ dx along the curves c ± and considering the limit of the difference of the result as c ± → c, we obtain Unlike the disclination-strength t, (7) shows that the dislocation-strength, lk dx, which equals the jump in the phase field across S, is not a well-defined topological constant (independent of the loop l) when the disclination-strength t = 0. This is analogous to the statement for layered structures that [θ] γ is not well-defined for a closed loop γ unless the enclosed number of flips is even. When the disclination-strength t = 0, θ is a constant on S and the dislocation-strength is a well-defined topological constant. Another situation when this is so is when S is a plane normal to t. Of course, it is possible in these situations forθ to be continuous as well.
The assumption thatÃ is a continuous field on Ω merits further discussion based on what is physically observed in elastic solids in contrast with nematics, smectics and patterns. In the presence of a dislocation in an elastic solid, i.e. considering a terminating surface S on which a constant displacement discontinuity occurs, it can be seen that there is no jump in the limiting values of the displacement gradient field as the surface S is approached. Consider now Frank's celebrated solutions for the entire family of straight, nematic wedge disclinations of integer multiples of 1 2 strength (mathematically identical to the solution for the screw dislocation in a solid). In this case, writingk(x) = cosφ e 1 + sinφ e 2 , where (e 1 , e 2 ) is a fixed orthonormal frame and φ is a half integer multiple of the anglek(x) makes with e 1 (plus a constant), it can be checked that while across any surface S whose trace on the x 1 − x 2 -plane is a straight radial ray from the origin there is no jump in the limiting values of Dφ, there is a jump in the limiting values of Dk = ∂ φk ⊗ Dφ as S is approached. For an S that is not necessarily planar, e.g., in smectics and stripes where it can be chosen as a contourθ = mπ and k = (1 − e 2πiJ ) ν (cf. Eq. (4)) is along ν, the direction of the local normal to the contour,Ã has a nonzero distributional curl on S (cf. Appendix A).
Consequently, we also consider the case whenÃ is symmetric, smooth, and curl-free in Ω\S. Repeating the above arguments, we now find that and In case Ã is of the form a ⊗ ν where a is a vector field on S and ν the unit normal field S, we note that the jump k is still constant on S, even thoughÃ is possibly discontinuous across S and the jump inθ still satisfies (7). All the computational examples in this paper satisfy this condition which is, however, not a necessary feature of our theoretical or computational formalism.
Motivated by the simple arguments above, it is clear that in the class of defects that are integer multiples of 1 2 , our formalism is adapted to representing the ± 1 2 -strength defects as fundamental (based on kinematic grounds) and all others are necessarily represented as composites of these fundamental defects. We compute examples of the field of such composite defects in Secs. 4.2.2 and 4.2.3.
3.1 Continuum mechanics of defects in scalar fields: natural patterns, smectics, and nematics (NPSN) Let us denote the region of the interior hole (with boundary) in Fig. 4 as the 'core' C. The considerations of Sec. 3 show that the phase field is in general discontinuous in non simply connected domains or, alternatively, if the fieldÃ was prescribed in the simply connected domain Ω ∪ C, but now with non-vanishing curl supported in the core C. It can also be seen, by considering Ω to be a punctured domain (i.e., C consists of a single point in 2-d or a curve in 3-d), that in such situationsk = Dθ is not in general a square-integrable field on Ω ∪ C -for simplicity, consider the case when t = 0 and the constant jump θ = 0. Even when C is a set of full measure, Dθ is not an integrable field on Ω and therefore its presence in any governing pde would be problematic in the presence of defects (characterized by non-curl-freeÃ and/ork fields). Consequently, we think of allowing fields with at most (smoothed) bounded discontinuities and removing all (smoothed) concentrations, referred to as 'singular' parts, from their gradients, these rehabilitated 'gradients' being called 'regular' parts. Roughly speaking, in classical governing equations for these phenomena, developed for situations not containing defects, we admit the appearance of only the regular parts of fields. We are also interested in modeling possibly large collections of moving defects that interact and possibly intersect, and tracking the topology of the defected body with cores modeled as excluded 'cylinders' is clearly impractical. Thus we seek a model that can be posed in simply connected domains, but nevertheless is descriptive of the topological properties of the line defects we are interested in. With this understanding, we consider a phase field θ with the regular part of its gradient denoted by k. The regular part of Dk in turn is denoted by A, with singular part by B so that A = Dk − B. The terminology of 'singular' is in the sense described above; when viewed at a microscopic scale these are smoothed concentrations on sets whose far-field identities are those of lower (than 3) dimensional objects.
The energetic physics of the nematic director or the pattern phase gradient (far from onset of roll instabilities) is based on energetic cost of director gradients. The director has head-tail symmetry, and this is modeled by assigning null energetic cost to values of the field B (the singular part of Dk) that reflect local changes in director orientation by 180 • over a small coherence length, typically of the order of a linear dimension of a disclination core. The terminating 'curve' of a 'surface' on which B is non-vanishing, say a constant, is a region where −curl B = curl A =: π is supported, and such a region corresponds to a disclination and we refer to the field π as the disclination density.
With reference to Fig. 4, if the disclination density field π were to be supported in the region C, then a curl-free field A and discontinuous fields k and θ satisfying (5) without the˜can certainly be defined. Moreover, if region C were to contain two separate concentrations of the disclination density of opposite sign, i.e. a disclination dipole, such that s πν da = 0 where ν is the unit normal field to any surface s that transversally intersects C, then the field θ would have a constant jump on any admissible surface S. As well, if the field A were to vanish for the moment and curl k were to be supported in C, then again θ would have a constant jump on any S. And, of course, if a θ field had to be defined at least locally in some region, curl k would have to vanish therein. When modeling smectics and natural patterns, we will energetically penalize curl k strongly and refer to regions that contain concentrations of curl k as a dislocation and the field −curl k =: γ as the dislocation density (this is slightly different from the definition for a solid). The considerations of Sec. 3 suggest that a region containing a disclination dipole may also be considered as an effective dislocation in that the far-field topological identity of both cases, measured by integrating the director field along any closed loop encircling the region, has to be a constant.
The definition of the disclination density field as a curl, i.e., π = curl A, associates a 'charge' with any closed curve l in the body, given by s πν da where s is any surface whose bounding curve is l, ν being the unit normal field on the surface. The flux of this charge across the bounding curve can also be kinematically characterized and is given by −π × V , where V is admitted as a velocity field, relative to the material, of the disclination density, resulting in the conservation laẇ These considerations lead to the following kinematics of our model: where (10) 5 follows from (9) up to a gradient. If 'grain boundaries,' or thin regions with a 2d skeleton on which B has a concentration, are allowed to move transversely to themselves with velocity V ⊥ , then (10) 5 would be modified to readḂ = curl B×V +grad BV ⊥ (without disturbing (9)). For simplicity, we do not consider this extra mechanism in this paper. A sufficient condition for the construction of a scalar field θ, corresponding to the fields k, B, is that it be possible to remove the support of the disclination density and dislocation density fields from the body and the resulting body (say Ω) be amenable to being rendered simply connected by the removal of a connected surface in it. In that case, a field θ may constructed satisfying Dθ = k in Ω which is generally discontinuous. In general, it is unclear if there is a single connected surface whose removal will permit the construction of a single valued phase θ on the complement. An example is the random stripe pattern in Fig. 2(d). Even if Ω can be rendered simply connected by removing a surface, it is natural to expect that there would be more than one surface with the same property and each such surface would correspond to a different θ field on Ω. Thus, when θ can be constructed, on an appropriately 'reduced' domain, in the presence of dislocations and disclinations, it can be expected to be 'massively' non-unique. This non-uniqueness of θ is a price one has to pay in going from a microscopic model for the system, which resolves behaviors on the scale of the underlying periodicity, to a macroscopic phase description.
Let ξ be a length-scale corresponding to the linear dimension of a disclination core. A typical example of a free-energy density function for the model is The vector field k is physically non-dimensional. The material constant K characterizes the elasticity of director gradients, P 1 ξ 2 K 1, P 2 K 1 are penalizing constants, f is a nondimensional, nonconvex function of |B| of the type defined below, α > 0 is a nondimensional number that tunes the strength of the nonconvexity of f , K * ξ 2 K ≈ 1, and ε Kξ 2 ≈ 1. Since the function f reflects the head-tail symmetry of the director, there should be approximately vanishing elastic cost for pointwise values of the director gradient of the type grad k ≈ n−(−n) aξ ⊗ l, where 0 < a ≤ 1 and n, l are unit vectors, the latter representing the direction along which the jump of n occurs, and therefore the two wells of f should be at |B| = 0, 2 aξ (based on rough energy-minimization arguments disregarding constraints of compatibility on Dk). Finally, to model pure nematics one sets P 2 = 0.
In what follows, we refer to a simply-connected domain or body within which the mechanics of interest takes place as Ω.
Beyond the specification of the energy density of the system, the disclination velocity V has to be specified. Guidelines for that specification arises from demanding that, up to contributions from the boundary of the body, the evolution of B results in a decrease in the total energy Ω ψ dv. We consider a free energy density with the following dependencies: ψ(k, Dk, B, π), noting that −curl k = grad k : X. Theṅ Consequently, requiring with M B a positive, scalar, mobility constant, is sufficient for the contribution to the rate of change of the total free-energy of the body to be non-positive due to the evolution of the field B, up to contributions from the boundary. For nematics and smectics, the first term on the right-and-side corresponds to increase or decrease of stored work in tune with externally supplied power, where k satisfies the balance of angular momentum (here, we are considering no material motion, in which case balance of linear momentum would be involved, along with viscous dissipation and material inertia [Les92,Ste04]). Convection patterns, in contrast are driven by organized, collective motion of materials, and k is not directly associated with a conserved quantity, i.e. mass, momentum or energy. Nonetheless, the late stages in the evolution of a convection pattern can be written as the gradient flow for the reduced Cross-Newell energy (3) that can be expressed in terms of k. When the primary concern is to understand dynamics close to local minima of the system free-energy, it therefore suffices to consider, in addition to (13), the following 'gradient flow' governing k:

Continuum mechanics of defects in vector and tensor fields: elastic solids
For defects in scalar fields, it sufficed to consider k as a fundamental field and consider only the regular and defect parts of Dk. This was primarily dictated by the nature of the energy density function for such systems, in particular elasticity arising due to director gradients, with a nonconvex contribution accommodating energetically preferred states of Dk.
To understand the similarities and differences between defects in elasticity of solids and NPSN, it is useful to first consider 'anti-plane' deformations of elastic solids, i.e., a body undergoing displacement in the out-of-plane direction as a function of in-plane coordinates. Then the displacement vector field has one non-trivial component, analogous to the phase field θ. A primary difference, θ, p, B fundamental kinematic fields p defect part of displacement gradient =: plastic distortion k = D θ − p regular part of displacement gradient =: elastic distortion B defect part of Dk =: eigenwall field A = Dk − B = D 2 θ − Dp − B regular part of elastic distortion gradient π = curl A = −curl B g.disclination density γ = A : X + B : X = curl p dislocation density (customarily defined as A : X) . p = −curl p × V γ evolution of plastic distortioṅ B = curl B × V π evolution of eigenwall field however, arises, from the energetics. In elastic solids, the primary elasticity arises from displacement gradients and it is important to consider regular and defect parts of the displacement gradient. Moreover, crystal periodicity dictates energetically preferred displacement gradient states. Hence, the kinematics of defected 'anti-plane' elastic solids is given in Table 3.2. A typical energy density function for defects in 'anti-plane' elastic solids (screw dislocations with Burgers vector and line direction in the out-of-plane direction, and twist disclinations with axis and rotation vector in the in-plane directions) looks like where µ is the elastic shear modulus, K is a modulus related to couple-stress elasticity, g 1 is a non-convex function reflecting preferred strain states due to lattice periodicity, g 2 represents a non-convex grain boundary energy reflecting preferred lattice misorientations, ε 1 is a material parameter characterizing dislocation core energy, and ε 2 is a constant characterizing (g.)disclination [ZA18, ZAP18] core energy. The balance laws of linear and angular momentum (involving secondderivatives in time) provide the governing equations for the evolution of θ, and Table 3.2 8,9 represent the evolution of the fields p, B. Constitutive guidance for V γ , V π for closing the model are deduced from thermodynamic arguments following similar argument as in deriving (13) [AF15]. A primary difference between elastic solids and NPSN is reflected in the scaling µ Kξ −2 , where ξ is assumed to be a typical linear dimension of a core for NPSN defects. The elastic modulus µ depending on p allows the modeling of earthquake dynamics [ZAWB15].
In 3D elasticity, all fields in Table 3.2 are tensors of one higher order than for the anti-plane case (and the elastic modulus is a 4 th -order tensor). Nonlinear elasticity requires the energy density to depend on k T k (where k now is the elastic distortion field) [ZAP18, AA19], and it can be a smooth function which is at least rank-one convex. And, of course, elasticity with defects in solids is fundamentally about material deformation and motion and singularities (at a macroscopic scale) in such deformation.
Finally, we note that if P (t) is the power supplied to, and K(t) the kinetic energy of, the body at time t, then our thermodynamic formalism ensures that P ≥K +Ė for all t, so thatĖ is not necessarily ≤ 0, allowing for externally driven, strongly out-of-equilibrium phenomena involving rapid material motion in our model.
That the type of model discussed above is realistic for elastic solids and earthquake rupture dynamics, even at the level of being robustly computable in dealing with objects that are macroscopically viewed as nasty singularities is demonstrated in [ZZA + 16, ZA18, ZAP18, ZAWB15, AF15,GAM15]. Connections of such models to NPSN are shown in [ZZA + 16] and alluded to in [NV17].
Similar models applicable to NPSN, with intriguing analogies to cosmology and the Standard Model of particle physics are discussed in [New12,NV17].

Illustration of theory
In this Section we demonstrate salient equilibrium features of the theory developed above through particular examples.

Nondimensional gradient flow dynamics
To non-dimensionalize Equations (14), (13), and (11), we introduce the following dimensionless variables, and assume M 1 = M 2 ξ 2 without loss of generality here, since the gradient flow equation for k can be treated as simply a device to get equilibrium of k with B fixed. The non-dimentionalized gradient flow equations of the energy (11) read as: +˜ (e tmn e tsjBin,ms ).
For convenience, we remove all tildes in remaining work and use the following nondimensional evolution equations in the rest of the paper: Note that the evolution equation (16) for B is different from (13). It is shown in [ZZA + 16] that while the L 2 -gradient flow dynamics (16) for the energy density (11) describes defect equilibria well, it is not able to adequately describe defect interaction and evolution in important situations, e.g. the elastic interaction and annihilation of a pair of positive and negative half-strength disclination. On the other hand, (13), a dynamics based on kinematics of topological charge conservation and thermodynamics, succeeds in this task, as demonstrated in [ZZA + 16]. While a theoretical explanation for this inadequacy of the gradient flow dynamics for these co-dimension 2 defects remains to be addressed (speculation is provided in [ZZA + 16]), with possible relation to similar phenomena for the equal well-depth case for the co-dimension 1 case analyzed in [RSK89] (recognizing that the mutual elastic interaction of co-dimension 2 defects is much stronger than for co-dimension 1), in this paper we simply rely on the gradient flow dynamics to predict approximate equilibria, and reach physical conclusions based simply on comparisons of total energy content of various defect configurations.

Computational Examples
In this Section, we assume a = 1, ξ = 0.1 and the size of the domain to be 20ξ × 20ξ. Unless specified otherwise, for all results pertaining to modeling nematics, we use the following default values for the (non-dimensional) material constants: P 1 = 100, P 2 = 0, α = 10, K * = 5, and = 1.
In the following computed examples, k is specified at one point (that eliminates rigid translation in pure statics), along with the natural boundary condition corresponding to (16) 1 , ((Dk − B) · n + P 2 curl k × n) = 0, where n is the outward normal to boundary of the domain. Also, the natural boundary condition for (16) 2 , curl B × n = 0 is applied. The one-point specification of k allows the prediction of distinct director patterns for defects with identical magnitude of strength, e.g., the 'target' and the 'source' for the strength +1 defect, utilizing identical B fields, as well as the ± 1 2 defects which involve initial conditions on B with differing sign. In the following calculations, results from gradient flow (16) with both k and B evolving are referred as equilibrium, and results of evolving k with specified B are referred as constrained equilibrium. In the constrained equilibrium calculations, B is not evolved from its specified initial condition. The acceptance criterion for a (constrained) local equilibrium state for all calculations is |Es−E s−1 | E s−1 ∆s < 10 −5 , where E s is the total energy at discrete time s, and ∆s is the time step at time s. 4.2.1 Strength ± 1 2 defects As mentioned in Sec. 3, f (|B|) has two wells at 0, 2 aξ . In this section, we prescribe initial conditions for the gradient flow calculations for the B field as non-zero within a layer. For a positive half strength disclination, In addition to different initial conditions for B, the director field k is specified, for all times, at one point on the top of the boundary of the layer. Fig 6 shows numerically computed equilibria of k obtained from the gradient flow equations (16) and Fig 7 shows energy density plots, for both the positive half and negative half disclinations respectively. Noteworthy is the fact that although k dramatically changes both direction and length within the layer, the energy density localizes only around cores. Figure 8 shows the Frank energy contribution (K|Dk − B| 2 ) for the positive half disclination and the comparison between the Frank energy density along x 2 = 0 with the function 1/r 2 , where r is the distance of a point from the origin. The Frank energy density is also localized around the core, and it follows the 1/r 2 decay rate outside the core, while yielding finite energy density inside the core.
To illustrate the effect of the prescription of B, we model − 1 2 defect with a different layer field B as follows, 2aξ e 1 ⊗ e 1 + √ 2 2aξ e 2 ⊗ e 1 , if |x| < aξ 2 and y < 0 0, otherwise.          Fig. 9(b) and Fig. 9(c) are the static results for director k and energy density field respectively. Although the prescription of B is very different compared to Fig. 5(b), the static equilibria of k and energy density outside cores and the energy density core shapes are similar.

Strength +1 defect
A strength (+1) defect can be represented as a composite defect in our model by putting two + 1 2 defects close together, as shown in Fig 10(a). The initialization of the k field is shown in Fig  10(b). Since the strength one defect is energetically unstable, we increase α to 50 in this example to constrain the diffusion of B. Both B and k are evolved following the gradient flow dynamics (16). The equilibrium of the director field k and the energy density are shown in Figure 11. In this section the color legends of all energy density plots are normalized by their maximum values. Since the strength one defect is energetically unstable, it tends to split into two half-strength defects with opposite signs which, subsequently, repel each other. As mentioned in Section 4.1, the gradient flow dynamics (16) is not capable of capturing this behavior. Thus, we calculate the equilibrium configurations and total energies corresponding to two strength-half defects at different separation distances to approximate the strength-one defect splitting process. In Fig  12(a) and 12(b), we show equilibria of the k field corresponding to two opposite half-strength defects at specified distances. Fig 12(c) and 12(d) are their corresponding energy density fields. The total non-dimensionalized energies for the positive one strength defect is 1.812 × 10 5 . After being normalized by the positive one strength defect's total energy, the total energies for the two opposite half-strength defect configurations at small and large separation distances are 0.508 and 0.498, respectively. Thus, a pair of + half-strength defects are energetically preferable in comparison to a single strength +1 defect, and the elements of the pair repel each other.
In addition to the 'source' pattern of the +1 disclination, we also calculate the 'target' pattern of the same strength +1 disclination, whose director field k and energy density are shown in Figure  13. For the 'star' pattern, k flips horizontally across the layer, and B 12 is nonzero within the layer. For the 'target' pattern, k flips vertically across the layer, and B 22 is nonzero within layer.
We note that that in all the examples solved in this paper, the width of the layer(s) for the      specification of B through initial conditions can be made arbitrarily small without affecting the qualitative properties of the solutions.

Strength −3/2 defect
Here we demonstrate a −3/2 strength defect as another interesting example of the capability of our theory in modeling composite defects of higher strength. Figure 14(a) shows the |B| field for the initial condition B(x, y). The initial condition is a piecewise-constant field with three different constant values of B in the layers, all with |B| = 2. In this calculation, B is not allowed to evolve from its initial conditions (this is an energetically unstable defect, and we are simply interested in demonstrating a negative-strength composite here), and k evolves following the gradient flow equations (16) 1 . Figure 14(b) shows the constrained equilibrium of the director field. Figure 15(a) shows the non-dimensionalized energy density, where the layers are completely invisible and the energy density shows a strong, non-singular (by design) concentration at the core.

Defect loop in 3D
A square half strength defect loop in 3D case is demonstrated in this part. Fig. 16(a) shows the prescription of B given as follows, where d represents the half length of defect square side. And Fig. 16(b) shows the corresponding prescription of π field. The one-point specification of k is applied at (−1, −1, −1) in addition to      zero moment boundary condition. The initial prescription of k is given as otherwise.
(21) Fig. 17(a) shows the constrained equilibrium of director k nearby defect layer. Director colors represent the norm of projection on x axis, where red means k 1 > 0 and blue means k 1 < 0. In this case, directors on the upper surface are e 1 (red) while the ones on the bottom are −e 1 (blue). Fig. 17(b) shows an illustration of director constrained equilibrium on multiple x-z cross sections. In Fig. 17(a), we draw two circuits corresponding to ones in Fig. 17(c) and in Fig. 17(d) respectively. Fig. 17(c) shows the director projection on x-z section (y = 0) and Fig. 17(d) shows the director projection on y-z section (x = 0). To make visualization cleaner, we draw black solid arrows along the circuit in each figure, zooming in their director arrows at each point. If we follow a circuit from upper surface to bottom surface in x-z plane, the directors transit from e 1 to −e 1 by varying director and size in plane. On the other hand, the director transition happens out of plane if following a circuit in y-z plane. For example, in Fig. 17(d), black arrows represent director projections on y-z cross section with meaning director pointing out and ⊗ meaning director pointing in. Fig. 18 shows the energy densities at y-z section (x = 0) and at x-z section (y = 0), indicating that the energy is localized around core where π = 0.
In the disclination loop being considered, the segments parallel to the x-axis in Fig. 16(b) are of twist character, (with axis of rotation along the y-axis) and the segments parallel to the y-axis are of wedge character. The solution demonstrates that such a loop produces minimal far-field director distortion away from the loop (this is similar to what happens for the dispalcement field corresponding to dislocation loops in solids). Closer examination also reveals that the director field is planar -in the x-z plane -almost everywhere except very close to the loop, with spatial variation   in all three directions. It is instructive to compare this almost planar director field with the field of the loop computed in [PAD15] which is of twist character everywhere and comprises a (strictly) planar director field with 3-d spatial variation. As shown there, such a loop with the same character everywhere induces significant far-field director variations, implying significantly more total energy content than the wedge-twist loop computed in this paper.

Smectic boundary
A grain boundary in a smectic is a special 'canonical' pattern of layered material systems that arises as a first bifurcation from a homogeneous state under external forcing, giving rise to piecewise homogeneously oriented domains separated by the boundary. The boundary can further reduce its energy by inducing defects within it. The existence of layers in a smectic implies that deviations of the curl k field from 0 are strongly energetically penalized. Thus, the parameters used in this calculation are P 1 = 1, P 2 = 1, α = 50, K * = 5, and = 0.1. Smectic boundaries can be modeled as a series of point defect pairs. To compare energies for different defect configurations, we model the smectic boundary as two defect dipoles at various distance, as illustrated in Fig 19. For each configuration, B is initialized following the procedure in Sec 4.2.1, given as (22): 2 aξ e 1 ⊗ e 2 , if |y| < aξ 2 and x d2 < x < x d3 2 aξ e 1 ⊗ e 2 , if |y| < aξ 2 and x > x d4 0, otherwise, where x d1 , x d2 , x d3 , and x d4 represent x coordinates of four defect cores from left to right. Both B and k evolve following (16). Fig 20 shows equilibria of k corresponding to different defect dipole configurations.
Since layers are of more interest in understanding smectic boundaries, we calculate the phase field θ from k from the following equations, div Dθ = div k on B with Dirichlet boundary condition Dθ · t = k · t on ∂B, where t is the unit tangent field on ∂B (i.e., a Helmholtz decomposition of k). The contour plots of the phase field θ of grain boundaries without and with defects are provided in Fig. 21. In Fig. 21, black lines represent phase field layers, red dots represent + 1 2 disclination cores, and blue dots represent − 1 2 disclination cores. In this example, the total energy for the defect-free smectic boundary is 8.596 × 10 4 . The non-dimensionalized total energies for the defected boundary configurations with small, medium, and large inter-dipole separations are 0.892, 0.891, and 0.881 respectively, normalized by the total energy of the defect-free boundary configuration. In addition, we consider another smectic boundary representation where defect dipoles repel and finally move out of body leaving a through layer field across the body. Figure 22(a) shows the corresponding prescription of B. Instead of defect dipoles, B is nonzero within the entire layer. Figure 22(b) shows its static equilibrium of smectic layers. The normalized non-dimensionalized total energy for this static equilibrium is 2.2 × 10 −3 .
Thus, larger inter-dipole separations are energetically favorable and it indicates that defected configurations are more energetically stable in the considered example.

Discussion
Patterns (ordered microstructures) and defects (breaking of an ordered pattern) are ubiquitous in extended systems. They arise from the interplay between two "universal" mechanisms, the tendency of systems towards order as their energy/temperature is lowered, and the tendency towards disorder from entropic considerations and the likelihood of "getting stuck" in "local" metastable states, that precludes perfect ordering. These defects play a big role in the properties of extended systems and understanding the birth, disappearance and dynamics of defects is critical in explaining a range of phenomena from plasticity, solid-solid phase transitions, fracture, convective transport, and complex fluids. It is thus of considerable interest to develop modeling and numerical methods to explain, analyze and predict the behaviors of defects and the extended systems they live in.
Fortunately, defects have many universal features, independent of the underlying physical system, reflecting their topological origins. A primary thrust of this work is to exploit this universality to develop a modeling framework and associated numerical methods that are applicable to computing defect driven behaviors in a wide range of systems of interest in materials science and continuum mechanics. We describe a common language for defects in natural patterns, smectics, and nematics, that draws on classical ideas for defects in solids [Vol07,Wei01] which have been incorporated into practically computable modern continuum mechanical theory recently [AF15, ZA18, ZAP18, AA19, AV19]. In this language, we develop a modeling framework that captures the dynamics of defects in terms of integrable energy densities, an important consideration for having a good numerical formulation. Our models can handle order parameters that have a head-tail symmetry, i.e. director fields, in systems with a continuous translation symmetry (e.g. nematic liquid crystals) and in systems where this symmetry is broken and replaced by a discrete translation symmetry (e.g. smectics and convection patterns). The framework we develop gives entire classes of models, and it allows for a natural incorporation of thermodynamic principles, the balance laws of mechanics and/or other physical principles like conservation laws for topological charges of defects.
We illustrate our methods with explicit computations for equilibrium configurations of nematic and smectic/pattern systems in 2-d (cf. Sec. 4.2) and in 3-d, illustrating the scope of our method. Now that we can adequately capture the equilibrium behavior of defects, an outstanding challenge is to capture the dynamics of defects. Some success has been shown in [ZZA + 16], but much remains to be done in this regard, for example, in modeling the process by which the random stripe pattern    Along the grain boundary, the red dots represent + 1 2 disclination cores and the blue dots represent − 1 2 disclination cores. Each ± 1 2 dipole may be considered a dislocation.  in Fig, 2(d) relaxes to the global ground state in Fig. 2(c). Another area of exploration is the coupling of our model to fluid flow to model liquid crystalline polymer flows with a high density of defects.
(a) Parameterization of k in terms of η and φ.  angle between the projection of k on e 1 − e 2 plane and e 1 axis. Similarly, we define φ as the angle between k and e 3 , ranging from −π to π.Thus, given a director k, η and φ can be calculated as η = arctan(k · e 2 , k · e 1 ) φ = sign(k · e 2 ) arccos(k · e 3 ), where arccos is the inverse cosine function whose range is from 0 to π, arctan(y, x) is the inverse tangent function returning angle ranging between −π and π whose tangent value is y x , and sign is a function returning the sign of φ. On the other hand, given a pair of (η, φ), the director k can be written as k = cos η sin |φ|e 1 + sin η sin |φ|e 2 + cos φe 3 .
Based on above parameterization, the jump of director k in a defect can be interpreted in terms of φ. For example, for a half strength defect, the director k changes its direction shown in Fig. 23(c), and the difference between φ 1 and φ 2 is π. To demonstrate the connection with planar cases discussed in [ZZA + 16], we adopt same notation λ to represent director discontinuity in φ. Then the regular part of elastic distortion gradient A can be written as It is easy to verify that A = Dk in defect-free cases where λ = 0. In addition, ∂ φ k can be calculated as ∂ φ k = tanh η cos η cos φe 1 + tanh η sin η cos φe 2 − sin φe 3 = k 1 k 3 tanh k 2 e 1 + k 2 k 3 tanh k 2 e 2 − |k − k 3 e 3 |e 3 .
The model above, while confirming to conventional intuition on thinking about disclinations in nematics and smectics, however is not, at least manifestly, invariant to the choice of the arbitrary orthonormal frame used in the definition of A.
When λ is of the form λ = a ν with support on a terminating layer around a surface with unit normal field ν and a is a constant, then it can be shown that curl λ = 0, except at the termination of the surface. In addition, if the k field does not have any longitudinal variations along the layer, then both the defect densities π and β = curl λ are localized at the same location. If k has longitudinal variations, then π is distributed all along the layer.