Frank-Read Mechanism in Nematic Liquid Crystals

In a crystalline solid under mechanical stress, a Frank-Read source is a pinned dislocation segment that repeatedly bows and detaches, generating concentric dislocation loops. We demonstrate that in nematic liquid crystals, an analogous Frank-Read mechanism can generate concentric disclination loops. Using experiment, simulation, and theory, we study a disclination segment pinned between surface defects on one substrate in a nematic cell. Under applied twist of the nematic director, the pinned segment bows and emits a new disclination loop which expands, leaving the original segment intact; loop emission repeats for each additional 180$^\circ$ of applied twist. We present experimental micrographs showing loop expansion and snap-off, numerical simulations of loop emission under both quasistatic and dynamic loading, and theoretical analysis considering both free energy minimization and the balance of competing forces. We find that the critical stress for disclination loop emission scales as the inverse of segment length, and changes as a function of strain rate and temperature, in close analogy to the Frank-Read source mechanism in crystals. Lastly, we discuss how Frank-Read sources could be used to modify microstructural evolution in both passive and active nematics.

In crystalline solids, the most common mechanism for generating dislocations is a Frank-Read source [1,2,24].Here, a pre-existing dislocation line is pinned at two points along its length, e.g. by impurities or entanglement with other defects.Under applied stress, the dislocation segment bows out into a circular arc.As stress increases, the arc expands and eventually self-intersects; a new dislocation loop snaps off, while the segment between the pinning points remains intact.This process can repeat periodically, leading to formation of a series of expanding concentric dislocation loops, resulting in plastic deformation of the crystal under stress.This continuous process may arrest if outgoing loops are blocked, creating a back stress on the source.
Here we explore an analogous Frank-Read mechanism arising in nematic liquid crystals.
In crystals, the Frank-Read mechanism has been extensively investigated to understand material susceptibility to plastic deformation [2].To gain insight into the Frank-Read mechanism in nematics, we perform experiments on a pinned disclination segment under twist deformation, and perform simulations of disclination loop emission under both quasistatic and dynamic strain conditions.Next we develop an analytic theory using minimization of a two-dimensional (2D) free energy, and a more intuitive theory using the Peach-Koehler force, which is the force driving defect motion due to local stress.Lastly, we discuss the potential role of engineered Frank-Read sources to drive heterogeneous loop nucleation and control microstructural evolution in passive and active nematics.

II. EXPERIMENT
To study a Frank-Read source in a nematic, we perform an experiment using a dynamic cell [25,26] where the in-plane rotation angle between top and bottom glass substrates can be changed to impose twist.The upper substrate, which is a convex lens with radius of curvature of two meters, has uniform planar anchoring, with a polyimide layer aligned via uniaxial rubbing.The lower substrate has an anchoring layer with two point defects of topological charge +1/2 and −1/2 [3,15,26], as shown in Fig. 1 (a), patterned via a photoalignment technique [27].The region with defects on the lower substrate consists of three areas with distinct orientations, surrounded by uniform planar anchoring, producing the two point defects [3,15].The gap between the substrates is filled with a nematic liquid crystal and a vertical disclination arch spontaneously forms, pinned at its ends by the two surface defects, as shown schematically in Fig. 1 (a).Viewed from above, the arch appears as a straight segment, as seen in a micrograph Fig. 1 (b).At 0 • of twist, the length of the disclination line is approximately equal to the spacing between the two point defects which is 100 µm.The gap between the two substrates is approximately 30 µm.
An in-plane twist is then imposed between upper and lower substrates.After we increase the relative orientation angle, twisting the nematic director, the disclination bows and expands laterally into a horizontal arc as shown in Fig. 1 (c-g).If the twist angle exceeds a threshold value, the disclination arc expands enough to self-intersect and snap off a new loop, as shown in Fig. 1 (h-i).The newly emitted disclination loop expands, leaving the original pinned segment intact.Experimental details are provided in Methods and in Ref. [26].A real time video of the experiment (Movie S1) is provided in Supplemental Information.
The angle of 144 • was the smallest rotation for which we found loop emission within a reasonable waiting time, here about nine minutes.If we stop at a slightly larger rotation angle, loop emission occurs within a shorter waiting time.For a rotation angle of 160 • , the waiting time drops to zero.This sequence of defect bowing, self-intersection, and loop snap-off is nearly identical to that observed for dislocation loop emission from a Frank-Read source in a crystalline solid; but the properties of dislocations in crystals and disclinations in nematic liquid crystals are not identical.Therefore, we perform simulations and theoretical analysis to explore the properties and dynamics of a Frank-Read source in a nematic.

III. SIMULATIONS: QUASISTATIC TWIST DEFORMATION
To gain insight into the Frank-Read mechanism in nematics, we perform numerical simulations of a nematic cell using the geometry shown in Fig. 2, under conditions of quasi-static twist.The bottom substrate has strong anchoring with two point defects of topological charges +1/2 and −1/2, similar (although not identical) to the experiment shown in Fig. 1 [9, 28,29].The top surface has strong anchoring in a uniform planar direction.The four side surfaces are free boundaries.The spacing between surface disclinations is w = 4 (arbitrary units), with cell size 120 × 120 × 10.Inside the cell, orientational order is characterized by the nematic order tensor Q(r, t), with free energy given by a Landau-de Gennes series expansion.We implement relaxational dynamics using an overdamped equation of motion, solved numerically via time-dependent finite-element modeling, using the software package COMSOL.Details are given in the Methods section.
Inside the cell shown in Fig. 3a, the top surface planar anchoring is initially aligned with the y axis, which is also the far-field orientation of the bottom surface.We initialize the order tensor as a uniform nematic state aligned with the y axis, and it relaxes into a stable configuration with a single disclination line connecting the two surface defects in an arch.Because the surface anchoring imposes parallel alignment with the two substrates, the resulting director field is also parallel to both substrates to minimize the elastic energy.As a result, the rotation axis of the director around the disclination ( Ω) is always perpendicular to the substrates, causing the disclination line to transform continuously from a +1/2 wedge disclination, to a twist disclination, and eventually to a −1/2 wedge disclination along its length.
Next we rotate the anchoring direction of the top surface through an angle δϕ.Loading is quasistatic; at each value of δϕ, we begin with the previous configuration of nematic order, and run the simulation until the nematic order relaxes to a new steady state.As δϕ increases from 0 • to 150 • , the disclination bows laterally into a curve, as seen in Fig. 3a-d.When δϕ reaches 150 • , the curve is approximately a semicircle.The disclination separates regions with different twist angles between bottom and top substrates.As discussed below, it is stabilized by the balance between the line tension (which favors a short disclination) and the twist energy (which favors moving the disclination to enlarge the low-twist region and reduce the high-twist region).We note that since the rotation of the top substrate does not introduce any z component into the director field, the Ω vector of the bowed disclination line [9] remains perpendicular to the substrates, implying that the topological feature of the bowed disclination line is still the same as the one for δϕ = 0.
As we increase applied twist, the system eventually reaches a threshold where the behavior changes radically.At 160 • , the disclination moves through the series of shapes shown in Fig. 3e-j via overdamped relaxation.Here the disclination does not find a stable or metastable position.Instead, it expands; arms adjacent to the two fixed ends rotate toward each other, intersect, and merge (Fig. 3g), snapping off to form two new disclination lines (Fig. 3h).One of these is a newly emitted pure twist disclination loop which expands (Fig. 3i) and eventually moves out of the simulation box by annihilating at the free boundaries (Fig. 3j).Through emission, expansion, and escape of this loop, the net twist deformation in the cell is reduced by π radians, which demonstrates a mechanism of stress release in nematic liquid crystals by generation and motion of topological defects.The remaining disclination segment connects the two surface defects, and this source is ready to emit more loops under further imposed twist.A new loop is emitted for each additional 180 • of applied twist.An animation of the simulation (Movie S2) is provided in Supplementary Information.
Because our model has mirror symmetry about the vertical plane crossing the two defect points on the bottom substrate at δϕ = 0, and the nematic liquid crystal is achiral, we expect the same behavior for the pinned disclination line when δϕ is negative, except that the pinned disclination line bows towards the opposite direction.

IV. SIMULATIONS: FINITE TWIST RATE AND THERMAL EFFECTS
To examine effects of finite twist rate, we simulate a nematic liquid crystal in a cell with the same geometry as Fig. 2, using a Lebwohl-Lasher rotor model [30][31][32].The liquid crystal is represented as a uniaxial nematic with a unit director field, ni , defined on each site, i, of a cubic lattice, interacting through the nearest-neighbor potential The director field evolves via overdamped rotor dynamics, relaxing along the potential energy gradient.Temperature is controlled via a Langevin thermostat.We implement the model for fast execution using graphics processing unit (GPU) acceleration.
Model details are provided in Methods.
Figure 4(a-c) shows disclination loop emission at finite temperature and under continuously increase twist angle; color indicates local potential energy density.As twist increases, the disclination bows and expands into an arc.At a twist angle δϕ * , the arc self-intersects (Fig. 4b), snapping off to form a new loop which expands and moves out of the simulated cell via annihilation at the free surfaces (Fig. 4c).We note that the Lebwohl-Lasher model introduces anisotropic line tension which causes curved disclinations to facet, an artifact of the underlying cubic lattice.An animation is shown in Supplementary Information (Movie S3).
We vary the applied twist rate of the top surface anchoring orientation.As shown in Fig. 4d, the threshold for loop emission δϕ * (identified here as the angle corresponding to peak potential energy) rises with increasing twist rate.This finding is in qualitative agreement with the behavior of Frank-Read sources in crystalline solids, where the strain threshold for loop emission typically rises with increasing strain rate [33].
Next we model thermal effects.We performed simulations at different temperatures and find that the threshold angle for loop emission rises with temperature (Fig. 4h).This increase in strain at loop emission as a function of temperature is qualitatively similar to behavior of Frank-Read sources in crystalline solids [34].
As another measure of thermal effects, we study the behavior of a bowed disclination segment under gradually increasing temperature.We hold the twist angle at a fixed, subcritical value where the bowed disclination segment is metastable.As temperature rises from k B T = 0.1ϵ to 0.5ϵ, the disclination arc expands outwards (Fig. 4e-g).This spontaneous expansion is similar to the disclination "blowout" at the nematic-isotropic transition, reported in previous Lebwohl-Lasher simulations [30], but occurs within the nematic phase.It is also consistent with the change in disclination shape on heating observed experimentally by Modin et al [35].
We note that the Lebwohl-Lasher model only takes account of director relaxation without hydrodynamics.In future work on finite strain rate studies, we will add hydrodynamics [11,22,36,37] to investigate backflow effects in loop emission and resulting rheological response.

V. ANALYTIC THEORY: 2D FREE ENERGY
For an analytic theory of the Frank-Read source, we consider the geometry in Fig. 5a.
We assume the nematic director lies in the (x, y) plane so that n = (cos ϕ, sin ϕ, 0).On the lower substrate, it is anchored at orientation ϕ B (x, y), given in Methods, which has surface defects at (±w/2, 0).On the top substrate, it is anchored at ϕ T = (π/2) + δϕ, with δϕ ≥ 0, independent of x and y.This geometry induces formation of a disclination shown by the thick red line in Fig. 5a.
We approximate the disclination shape as a circular arc between the surface defects, passing through (0, y top ).For y top = 0, this shape is a straight line between the two surface defects.For y top = w/2, it is a semicircle connecting the surface defects.For y top → +∞, it is an almost-complete circle in the positive y region, just touching the surface defects.
The free energy of this configuration has two components.One component is the line energy, proportional to the length of the disclination.The second component is the elastic energy associated with director twist between the two substrates.As y top increases and the disclination grows, the line energy increases but the elastic energy decreases, because more area is transformed from high twist to reduced twist.The balance between these two competing terms is worked out in detail in Methods.The minimum free energy occurs at which implicitly shows how disclination shape evolves with δϕ.Here, ȳtop = y top /w is a dimensionless parameter characterizing y top normalized by defect separation w.Likewise, α = (Id)/(πK 22 w) is a dimensionless parameter representing line tension I times cell thickness d, compared with twist elastic constant K 22 and defect separation w.
Hence, the disclination grows from a straight segment at δϕ = 0 into a semicircle, then into an almost-complete circle connecting the surface defects at δϕ = π/2.All of these shapes are stable minima of the free energy.If δϕ rotates beyond π/2, then there is no stable shape in the positive y region.Rather, the disclination must grow to infinity.
If α > 1/4, the behavior is different.As δϕ rotates from 0 to π/2, y top increases from 0 to a finite value, which depends on α.If δϕ rotates beyond π/2, then y top continues to increase, and the disclination grows larger in the positive y region as a metastable structure (while the stable structure is a disclination in the negative y region).The metastable disclination can grow up to If δϕ rotates beyond δϕ max , the disclination becomes unstable; it finds no steady state in the positive y region.Instead, it must expand to infinity.
Predictions of this analytic theory can be compared with the quasistatic nematic-ordertensor simulation presented above.Figure 5c shows simulation results for defect separation w = 4, 8, and 18.For each value of δϕ, we analyze the simulated shape to extract y top .These results can be fit to Eq. ( 1), with the dimensionless parameter α = 0.908, 0.454, and 0.202.
The prediction and simulations agree well over the full range of δϕ, and they give consistent values for the limit of metastability δϕ max .As expected, α is inversely proportional to w.
From the values of α and w, along with the cell thickness d = 10 and twist elastic constant K 22 = 3.0, we extract the line tension I = 3.4 (all in the arbitrary units defined in Methods).
Line tension is assumed to remain constant throughout the dynamic expansion process.This value is reasonable in terms of Landau-de Gennes theory.
The predictions can be used to interpret the experiment presented in Fig. 1 and Supplementary Video S1.Because the disclination line expands rapidly past δϕ = 144 • , we estimate the dimensionless parameter α = (Id)/(πK 22 w) = 0.9.Because the surface defect spacing is w = 100 µm, and the cell gap is d = 30 µm, we have the ratio I/K 22 = 10.
For example, the twist elastic constant might be K 22 = 4 pN, and disclination line tension I = 40 pN.These numerical estimates can be compared with a recent experimental and theoretical study by Modin et al [35].Our dimensionless parameter α is equivalent to their parameter 2γ, and our estimate I/K 22 = 10 is similar to their reported ratio of 12 to 18, depending on temperature.

VI. ANALYTIC THEORY: PEACH-KOEHLER FORCE
In crystalline solids, the standard theory of a Frank-Read source is based on the balance between the Peach-Koehler force (which pushes the dislocation forward) and the dislocation line tension (which pulls it back).The Peach-Koehler force is a force per length acting on dislocations in crystalline solids under stress [1,2,38].It can be written as where b is the Burgers vector, σ is the local stress tensor, and t is the local tangent vector of the dislocation.An estimate of its magnitude is f PK = bσ.By contrast, the line tension provides a force per length of I/r, where I is the dislocation energy per length and r is the local curvature radius.These forces balance when bσ = I/r.As an estimate, the dislocation energy per length is I ≈ Gb 2 , where G is the shear modulus of a solid.Furthermore, the minimum curvature radius is r ≈ w/2, where w is the spacing between the pinned defects.
Hence, the curved dislocation is only stable up to a maximum stress of σ ≈ 2Gb/w.Beyond that stress, the Frank-Read source becomes active and produces dislocation loops.
We can apply the same argument to the Frank-Read source in a nematic liquid crystal.
Early work by Kléman [3] showed that nematic disclinations experience a Peach-Koehler force, but this effect has not been studied much in the liquid-crystal literature.In a recent paper [9], our group developed a theory for this force as Here, Ω is the rotation vector of the director field around the disclination, with π Ω analogous to a Burgers vector.The effective stress tensor σ eff is the difference of torques between the Hence, the curved disclination is only stable up to the maximum effective stress of σ eff ≈ 2Kπ/w.Beyond that stress, it produces disclination loops.Critical stress for dislocation loop emission in crystalline solids also scales inversely with the distance between pinning points.
To make our theory more precise, we return to Eq. (3).The tangent vector t runs along the disclination.It is in the (x, y) plane, except in small regions near the two surface defects.The rotation vector Ω is vertical everywhere.We must make an arbitrary choice for the signs of t and Ω, and these two signs must be compatible with each other.Let us choose that t goes from the +1/2 to the −1/2 surface defect, and Ω = + ẑ.From Ref. [9], assuming equal Frank elastic constants, the effective stress tensor is As a simple approximation, we consider only the variation of n in the z direction.The effective stress of rotating the top alignment from ϕ T = π/2 to ϕ T = (π/2) + δϕ then becomes , with all other components zero.Hence, the Peach-Koehler force per length acting on the disclination is In our geometry, t is mostly in the + x direction, and hence f PK is mostly in the + ŷ direction.
It pushes the disclination toward the positive y region, as seen in both experiment and simulations.By comparison, line tension provides a force per length of with direction inward toward the center of curvature.By balancing these forces, we predict the curvature radius Using expressions for r and alignment angle ϕ B from Methods, we see that this equation is exactly equivalent to Eq. (1) derived from free energy minimization.It predicts the maximum stable stress of σ eff zz = (2I)/(πw), consistent with the estimate above, in the limit of α = (Id)/(πKw) ≫ 1.
One further benefit of the Peach-Koehler force concept is that it provides a way to estimate the velocity of an expanding disclination loop.In general, a disclination moves with overdamped dynamics, following the equation of motion f total = ζv, where f total = f PK + f tension is the total force (per length), ζ is the drag coefficient (per length), and v is the velocity.In Methods, we use this dynamic equation to determine the velocity of a large expanding disclination loop.We find v = 0.81(δϕ − π/2), in the arbitrary units.For example, δϕ = 160 • gives v = 1.00, and δϕ = 108 • gives v = 0.26.
To check this estimate, we monitor loop expansion after δϕ exceeds the critical threshold, as in Fig. 3e-j, tracking y top as a function of time.Figure 6a shows the motion in the simulation with defect spacing w = 4 at angle δϕ = 160 • .After an initial phase of slow growth, presumably limited by line tension, y top increases at the steady velocity of 1.10 arbitrary units.Figure 6b shows corresponding results in the simulation with w = 8 at δϕ = 108 • .In this case, y top increases at the steady velocity of 0.256 arbitrary units.These simulation results are quite consistent with the prediction.
We note that the present analysis using the Peach-Kohler force was carried out under the approximation of equal splay and bend Frank elastic constants.In this limit we find that the analogy to the Peach-Kohler mechanism in crystalline solids works well.The case of unequal Frank constants is more challenging and will be the subject of future investigation.

VII. DISCUSSION AND OUTLOOK
In closely related experimental and theoretical work, Jiang et al [39] studied a lattice of surface-pinned disclination segments in a nematic cell.Under director twist imposed by optically rewriting surface anchoring orientation on the opposite substrate, the defects underwent programmable transformations via expansion, merging, and reconfiguration.We interpret this geometry as an array of Frank-Read sources that are closely spaced, such that expanding disclination segments merge and reconnect, pre-empting loop emission.We predict that spacing the pinned segments further apart will enable emission of expanding loops, forming an array of Frank-Read sources.Expanding loops will then collide and reconnect.Independent of defect spacing, this mechanism creates a localized slip plane for twist deformation, in analogy to slip planes in metals.We speculate that this technique could be implemented to control microstructural evolution and rheological response (e.g.lubrication) for flow of nematics under shear between patterned substrates.The Frank-Read mechanism could arise more generally in nematics via other stimuli that deform the director field, such as fluid flow, electric or magnetic fields, or by internal stresses in active nematics.In nematic flow through a microfluidic channel [11,[40][41][42] or in a Couette cell [10], disclination nucleation has been observed at confining walls due to the presence of surface irregularities which pin disclinations.
In nematics, we have so far only considered Frank-Read sources with pinning points located on a confining surface.In crystals, by contrast, dislocations can pin at random along their length via jog formation, interaction with impurities, or entanglement [2, 33], generating new Frank-Read sources.Disclinations in pure nematics do not spontaneously pin by any analogous mechanism.However addition of colloid particles can induce disclinations to form knots and complex entangled networks [11], suggesting that a nematic-colloid composite might demonstrate spontaneous Frank-Read source formation under imposed flow.We speculate that resulting microstructural evolution might show complex behaviors in analogy to mechanisms governing plasticity in crystalline solids.On the other hand, thanks to the interaction between colloids and disclinations, the Frank-Read source also provides a controllable approach for redistribution of colloids immersed in nematic liquid crystals responding to external stimuli.From recent studies on manipulation of light using nematic disclinations [43][44][45], we speculate that disclinations controlled and created via the Frank-Read mechanism might find use in novel optical applications.
In 3D active nematics, active stress drives both fluid flow and director rotation.Studies by Duclos et al. [22] revealed 3D details of microstructural evolution.They showed that colliding disclinations merge and reconnect without forming long-lived pinning points.While they observed that disclination self-intersection can generate new loops, in the absence of defect pinning, no Frank-Read source formation was observed.In future work, we will investigate whether a network of Frank-Read sources can be created in an active nematic by inscribing defect pinning points on confining substrates.
It is also interesting to note that the homogeneous defect loop nucleation mechanism in active nematics, as shown in Duclos et al. [22], is not typically observed in plastic deformation of ductile crystals.Dislocation loop generation from Frank-Read sources occurs at far lower critical stress, and thus heterogeneous defect nucleation preempts homogeneous nucleation.
Reasoning by analogy, we speculate that a high density network of Frank-Read sources in an active nematic could produce a more ordered pattern of defect generation and flow, with stresses never reaching the threshold for homogeneous loop nucleation.Any mechanism that produces disclination pinning-e.g.surface defects, colloids, or micropillars -could thus modify pattern formation and rheology of both active and passive nematics.
In conclusion, we have presented experimental, theoretical, and simulation studies of the Frank-Read mechanism in nematic liquid crystals.The experimental micrograph in Figure 1(b-i) shows the signature shape evolution of a Frank-Read source.Quasistatic simulations using a nematic order tensor and Landau-de Gennes free energy demonstrate the mechanism.Finite strain rate simulations via the Lebwohl-Lasher rotor model show that threshold twist for loop emission rises with increasing strain rate, analogous to behavior in metals.We also observe spontaneous loop expansion under increase of temperature.We develop analytical theory of loop emission by minimizing a 2D free energy, and also via balance of Peach-Koehler and line tension forces, showing close agreement with quasistatic simulations.
As a historical aside, we note that Sir Frederick Charles Frank made crucial contributions to the study of both crystal plasticity and nematic liquid crystals.He co-discovered the Frank-Read source mechanism together with Thornton Read in 1950 [24]; and also formulated the Frank Free energy for liquid crystals with contributions from splay, twist, and bend, in 1958 [46].The work described here brings together and builds upon Frank's contributions to both fields.

A. Dynamic Cell Experiment
The dynamic cell [25,26] is comprised of two glass substrates, top and bottom, with a nematic liquid crystal (in this case 5CB) sandwiched in between (Fig. 1a).The bottom substrate is fixed on a rotating stage whose rotation axis is along the substrate normal (zaxis).The top substrate is fixed on a stage which translates along the z-axis.This setup allows the bottom substrate to twist with respect to the top substrate, and for the top substrate to translate to change the cell gap between substrates.The cell gap is measured using a capacitance sensor.
The bottom substrate is a 1 inch × 2 inch glass substrate.We used a plano-convex lens (1 inch diameter; with a radius of curvature of two meters) for the top substrate.The top substrate is coated with PI-2555 polyimide and is rubbed 10 times with a velvet cloth and marked to indicate the rubbing direction.The bottom substrate is coated with an azo-based photoalignment compound (SD-1) that can be photoaligned by polarized UV light [47].
Photoalignment of the bottom substrate is done in two steps.First, the entire substrate is exposed to polarized UV light (8 mW/cm 2 for 1 minute) to induce uniform planar anchoring.Surface defects are then generated in the surface anchoring pattern using a maskless photoalignment system [25,27,48].
The patterned anchoring is comprised of four regions with different planar alignment angles (Fig. 1a).The central rectangular region (100 µm×50 µm) is split into three segments as shown in Fig. 1a.The maskless system uses a rotating polarizer to expose the substrate with polarized UV light.Each region is aligned in turn: first the polarizer is rotated to generate the desired alignment angle.Then the maskless system creates an image of the region focused on the substrate surface.
Once the substrate preparation is complete, the dynamic cell is assembled.First a small amount of 5CB is deposited on the bottom substrate surface.The disclination line (Fig. 1b) spontaneously forms between the two surface defects.Then the top substrate is translated downwards, put in contact with the liquid crystal, and is further translated downward until a cell gap of approximately 30 µm is attained.

B. Simulations of nematic order tensor
Nematic order is represented by the tensor Q(r, t), which is a function of position and time.In uniaxial regions, outside of disclination cores, this tensor is related to the nematic director n and scalar order parameter S by The free energy is given by the standard Landau-de Gennes series expansion where the first three terms are the bulk free energy, and the last term is the distortion free energy.We choose the coefficients A = −1, B = −12.3,C = 10, and L = 2.32 in arbitrary units.The ratios of these coefficients correspond to the ratios in the liquid crystal 5CB, provided that the unit length in our simulations corresponds to 4.45 nm.With these coefficients, the bulk scalar order parameter is the minimum of the bulk free energy S = 0.53.The single elastic constant L for the nematic order tensor implies a single Frank elastic constant for the director K = 9 2 S 2 L = 3.0.To model pure relaxational dynamics, we use the overdamped equation of motion for where Γ is the rotational viscosity coefficient in the tensor representation.We choose units of time such that Γ = 1.
In the dynamic evolution, we impose strong planar anchoring conditions on the bottom and top surfaces, with n = (cos ϕ, sin ϕ, 0).On the bottom, the director configuration is This pattern has a +1/2 defect at (−w/2, 0) and a −1/2 defect at (+w/2, 0), which enables a pinned disclination line connecting the two defects through the interior of the liquid-crystal cell.On the top, the director has the uniform alignment The alignment of the top surface is initially along the y direction, and δϕ is the angle by which it is rotated from the initial state.
Disclinations are visualized by the biaxiality parameter In particular, Fig. 2 shows the isosurface of β = 0.7.The thickness of the disclination corresponds to the core radius, which can be estimated from Landau-de Gennes theory as in the arbitrary units.We recognize that this value is unrealistically large, in comparison with the defect spacing and system size.By exaggerating this parameter, we reduce the ratio of length scales, which allows the simulations to run more quickly.That is a drawback of simulations of the nematic order tensor, which is addressed in simulations of the lattice model.

C. Simulations of lattice model
The Lebwohl-Lasher model is a simplified version of the order-tensor representation of a nematic liquid crystal, where we assume equal elastic constants, no biaxiality, and scalar order parameter of uniform magnitude.Spins on a N × N × N cubic lattice represent the director field of a nematic liquid crystal.These spins can rotate freely, but they do not have any translational motion [49,50].The Lebwohl-Lasher model is effective at simulating large systems on the order of 10 6 spins, where each spin could represent a region of similarly aligned molecules [30,50].Spins interact through the potential where i and j are neighboring lattice sites, and ϵ is a positive interaction constant [32,50].
The Lebwohl-Lasher potential effectively brings neighboring spins into parallel alignment and accurately reproduces disclination behavior of liquid crystal systems [30,31,50].Spins can be fixed at the boundary and patterned to represent a strongly anchored layer.In this study, we use the anchoring patterns of Eqs.(9-10) on the bottom and top surfaces, respectively, just as in the simulations of the nematic order tensor.
Many Lebwohl-Lasher models in the literature use a stochastic Monte-Carlo algorithm to study equilibrium behavior at finite temperature [30,32,50].However, we use an overdamped torque relationship between the neighboring spins to minimize our energy, creating a Lebwohl-Lasher rotor model [32,50].This dynamical system allows us to study microstructural evolution where the system continuously moves towards the lowest energy state.The torque between neighboring spins is calculated from the derivative of the energy with respect to orientation, To allow the system to overcome energy barriers at constant temperature and reach the global energy minimum, we also add another term to the torque associated with an angular-momentum-conserving Langevin thermostat scaled with a temperature parameter . This torque is then applied equally and oppositely to the neighboring spins ni and nj .
We set each spin's angular velocity proportional to its total torque, resulting in overdamped rotation.The angular velocity then becomes ω = Cτ , where C is a mobility factor inversely proportional to the rotational viscosity.Hence, the spins in the Lebwohl-Lasher in each time step dt.
We note that the underlying cubic lattice in the 3D Lebwohl-Lasher model gives rise to anisotropy in disclination line tension.The lattice favors disclinations running in orthogonal directions, leading to rectangular shapes.Additional interactions with the edges of the square simulation box enhance the rectangular shape of the disclination when it is allowed to completely relax.This rectangular shape is not as pronounced when under a constant twist rate.Furthermore, the motion of a disclination requires overcoming a small energy barrier between lattice sites, in analogy to the Peierls-Nabarro barrier for dislocations in crystals [52].

D. Minimization of 2D free energy
We consider a nematic liquid crystal in the geometry shown in Fig. 5a.The bottom and top substrates have strong anchoring with the configurations given in Eqs.(9-10), respectively.Hence, the bottom substrate has two surface defects at positions (±w/2, 0), and these defects are connected by a disclination line.We model this disclination line by a circular arc passing through (0, y top ).The total free energy is the sum of two terms: the line tension associated with the disclination, and the twist free energy of the director field between the bottom and top substrates.We wish to minimize the total free energy over the parameter y top .
The tension free energy is the line tension I multiplied by the arc length of the disclination, which is the circular radius times the angle subtended.The radius is r = (y top /2) + (w 2 /(8y top )), and the angle is θ = 4 arctan(2y top /w).Hence, the tension free energy is The twist free energy is associated with twist of the director between the bottom and top We combine Eqs. ( 15) and ( 17) to obtain the derivative of total free energy with respect to y top , and then set it equal to zero to find the minimum, giving the result in Eq. (1).
Formally, that equation gives δϕ as a function of y top .However, it can be inverted either numerically or graphically to find y top as a function of δϕ.
top and bottom surfaces, normalized by the area over which the torques are applied.Hence, a simple estimate of the Peach-Koehler force per length is f PK = πσ eff .It competes with the disclination line tension, which provides a force per length of I/r, and the forces balance when πσ eff = I/r.The disclination energy per length is approximately I ≈ Kπ 2 , assuming all Frank elastic constants are equal to K. The minimum curvature radius is r ≈ w/2, is half the spacing between the surface defects.

FIG. 1 .FIG. 3 .FIG. 4 .
FIG. 1. Experimental micrographs showing a Frank-Read source in a dynamic liquid crystal cell, with a disclination segment pinned at two point defects on the substrate, before and after an imposed twist strain.(a) Diagram of the dynamic cell.The surface of a convex lens (upper substrate) imposes uniform planar alignment.The anchoring pattern with two point defects is shown on the lower substrate.(b) Viewed from above, the arch appears as a straight line at t = 0 s.(c)-(i) Disclination shape evolution after the angle between the two substrates is increased via mechanical rotation by 144 • , showing the Frank-Read mechanism.The defect segment initially bows outward, shown at times (c) t = 210 s (d) t = 300 s (e) t = 400 s (f) t = 540 s (g) t = 740 s.(h) Self-intersection occurs at t = 770 s. (i) At t = 788 s, a new loop has snapped off, leaving the initial segment intact.

FIG. 5 .
FIG. 5. (a) Geometry for the theoretical free energy of the Frank-Read source.The thick red line represents the curved disclination, viewed from above.(b) Stability diagram in terms angle δϕ and dimensionless parameter α, indicating whether the curved disclination is stable, metastable, or unstable (so that it must grow to infinity).Along the gray dashed line, the shape is a semicircle.(c) Comparison of theoretical predictions with nematic-order-tensor simulation results, for different defect spacings w.

FIG. 6 .
FIG. 6. Motion of y top during the unstable expansion process, from nematic-order-tensor simulations.(a) Defect spacing w = 4 and angle = 160 • .(b) Defect spacing w = 8 and angle δϕ = 108 • .The black dashed lines show linear fits in the limit of y top ≫ w.
and its derivative with respect to y top is substrates, across the cell thickness d.It is difficult to calculate the total twist free energy, integrated over the entire liquid crystal, but easier to calculate its derivative with respect to y top .If the center of the circle moves upward from y top to y top + dy top , then a narrow arc is transformed from Region II to Region I. Remarkably, the anchoring orientation ϕ Combining these factors, the derivative of twist free energy with respect to y top is B is uniform along this arc, with the value ϕ B = arctan(2y top /w).In Region II, the twist angle is (ϕ T − ϕ B ), and hence the twist free energy density per area is (K 22 /(2d))(ϕ T − ϕ B ) 2 .In Region I, the twist angle is reduced to (ϕ T − ϕ B − π), and hence the twist free energy density per area is reduced to (K 22 /(2d))(ϕ T − ϕ B − π) 2 .Under this change in y top , the area transformed from Region II to Region I is