PHYSICAL REVIEW RESEARCH 6, 023284 (2024)

Approach and rotation of reconnecting topological defect lines in liquid crystal
Yohei Zushi ,1 Cody D. Schimming ,2 and Kazumasa A. Takeuchi

1,3,*

1

2

Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
3
Institute for Physics of Intelligence, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan(Received 1 April 2024; accepted 30 May 2024; published 17 June 2024)
Topological defects are a universal concept across many disciplines, such as crystallography, liquid-crystalline
physics, low-temperature physics, cosmology, and even biology. In nematic liquid crystals, topological defects
called disclinations have been widely studied. For their three-dimensional (3D) dynamics, however, only recently
have theoretical approaches dealing with fully 3D configurations been reported. Further, recent experiments haveobserved 3D disclination line reconnections, a phenomenon characteristic of defect line dynamics, but detailed
discussions were limited to the case of approximately parallel defects. In this paper, we focus on the case of two
disclination lines that approach at finite angles and lie in separate planes, a more fundamentally 3D reconnection
configuration. Observing and analyzing such reconnection events, we find the square-root law of the distance
between disclinations and the decrease of the interdisclination angle over time. We compare the experimental
results with theory and find qualitative agreement on the scaling of distance and angle with time, but quantitative
disagreement on distance and angle relative mobilities. To probe this disagreement, we derive the equations of
motion for systems with reduced twist constant and also carry out simulations for this case. These, together withthe experimental results, suggest that deformations of disclinations may be responsible for the disagreement.
DOI: 10.1103/PhysRevResearch.6.023284
I. INTRODUCTION

Where there is order, its mismatch exists. Topologicaldefects are prime examples of mismatching order and areomnipresent across many materials [1,2], including crystals[1], liquid crystals [1,3,4], superfluids [5], spacetime in cosmology [6], and living things [7,8]. Despite various orders
depending on systems, corresponding topological defects are
expected to share common behaviors. For example, lineshaped defects often experience reconnections [6,9–17] inwhich two defect lines approach, collide, exchange endpoints,
and separate.
Liquid crystals, especially nematic liquid crystals, are amaterial where topological defects have been heavily studied
because of their optical properties and controllability. Nematic
liquid crystals are typically composed of building blocks withelongated shapes which tend to align with each other. They arecharacterized by a unidirectional order represented by a unit
vector n̂ with an apolar nature (n̂ ≡ −n̂) called the director.
Its singularities are topological defects, with the most basic
kind being disclinations. Many studies on nematic disclinations include, for instance, the generation of defects followedby the ordering process [16], defects’ interaction with microparticles [18–21] or light [19,22], operation of molecules
by defects [23,24], and control of defects by alignment of*Contact author: kat@kaztake.org

Published by the American Physical Society under the terms of theCreative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)and the published article’s title, journal citation, and DOI.
2643-1564/2024/6(2)/023284(10)the surface of the containers [20,21,25–28]. In recent years,
nematic ordering and patterns have also been found in living
systems [7,8], and topological defects are suggested to berelated to some biological functions, including cell extrusion[29], promotion of bacterial colonies’ vertical growth [30,31],
and organization centers of Hydra morphogenesis [32]. As fortheoretical approaches to disclination dynamics, although thegoverning equations of liquid crystal are well-known [4], it is
not straightforward to obtain equations of motion of defects,
especially for fully three-dimensional (3D) configurations.
Recently, there has been some theoretical work to better understand the interaction and velocity of 3D topological defects[28,33–35].
Despite the intensive study of nematic disclinations, previous observations of disclination dynamics are mainly bytransmitted light, and only two-dimensional (2D) information
can be obtained. In our previous work [17], making use of
the dye localization at the disclination core [24] to visualize
disclination lines, we successfully observed its 3D dynamics,
in particular, reconnections and loop shrinkage [17]. However,
the analysis was limited to the dynamics of reconnections
occurring to essentially parallel disclinations in an approximately single plane.
In this paper, we study the dynamics of intersecting reconnections, in which two disclinations lying in separate planes
approach not in parallel but with an angle. We used confocal
microscopy and observed disclinations relaxing from an electrically driven turbulent state. We extracted the positions of
disclinations from acquired images, investigated the time evolution of the distance and the angle between two reconnecting
disclinations, and discussed the result by comparing it withtheory and numerics.

023284-1

Published by the American Physical Society

ZUSHI, SCHIMMING, AND TAKEUCHIPHYSICAL REVIEW RESEARCH 6, 023284 (2024)

II. EXPERIMENTAL METHOD(a) an intersecting type

We used a nematogen MLC-2037 (Merck, a discontinued product) because of small optical anisotropy (n =0.0649) for precise optical observation and negative dielectric anisotropy ( = −3.1) for generating disclinations, as
described below. We added 0.01 %wt of electrolyte, tetran-butylammonium bromide, to increase its conductivity and0.005 %wt of a fluorescent dye, coumarin 545T, to label
disclinations [17,24]. The sample was put into a 130 − µmthick cell consisting of a coverslip and a glass plate withpolyimide tape spacers. Both glass substrates were coatedwith indium tin oxide and imposed unidirectional planaralignment by polyvinyl alcohol coating and rubbing with velvet cloth. Similarly to our previous work [17], by applying
an ac voltage of root-mean-square amplitude 150 V and frequency 35 Hz, we induced a turbulent state called dynamic
scattering mode 2 [4,37] to generate a high density of disclinations. Then, by switching off the applied voltage, the system
exhibited a relaxation process in which disclinations interacted, shrinked, and disappeared.
To obtain 3D dynamics of disclinations, we used a laserscanning confocal microscope, Leica SP8 (objective 40x, NA
1.30, oil immersion) with resonant scanner working at 8 kHz
and a piezo objective scanner. After starting relaxation, we
first searched disclination pairs likely to reconnect by moving the field of view, which was fixed thereafter. Due tothe imposed homogeneous alignment at the cell surfaces,
disclinations moved away from the surfaces and were often
found near the midplane of the cell thickness. In such a circumstance, reconnection events are conveniently classified to
two types: an in-plane reconnection, in which two reconnecting disclinations are almost in a single plane and approach
approximately in parallel [Fig. 1(b)], and an intersecting reconnection, in which disclinations are not in a single planeand form a finite angle [Fig. 1(a)]. Here, we observed 12
intersecting reconnection events.
A series of experimental snapshots is shown in Fig. 1(c)(see also Supplemental Material videos 1–4 [36]). The fluorescent dye was excited at 488 nm by laser light polarized
in the direction perpendicular to the nematic easy axis. The
directions of the laser polarization and the nematic easy axis
are denoted by the x and y axes, respectively. The fluorescencesignal in the range between 500 and 600 nm was confocally detected by a photomultiplier tube detector (pinhole size
23 µm). The voxel size in the xy plane was 0.455 µm and thespacing between z slices was 1 µm. The number of voxels was
512, 96, 30 in the x, y, and z directions, respectively. The timeinterval between consecutive confocal images was 0.330 s. We
extracted the coordinates from the 3D xyz images by applying
the snake method [38] (see Appendix A).

III. EXPERIMENTAL RESULTS

We start with the minimum distance δ(t ) between two
reconnecting disclinations at time t [Fig. 2(a)]. The result in
Fig. 3(a) shows the square-root law,
δ(t )  C|t − t0 |1/2 ,

(1)

(b)

an in-plane type−2.0 s(c)

20 μm−1.0 s0s0.3 sFIG. 1. Reconnectionsofnematic
disclinations.
(a),(b) Schematics of an intersecting reconnection (a) and an in-plane
reconnection (b). (c) Fluorescence observation of an intersecting
reconnection. Top views (left column) and side views (right column)
are displayed. The reconnection occurred at time t = 0. See also
Supplemental Material videos 1 and 2 [36]. Another reconnection
event is also shown in Supplemental Material videos 3 and 4 [36].

where t0 is the reconnection time. This scaling is the same asthat observed for in-plane reconnections of essentially parallel
disclination pairs [17]. Theoretically, the scaling was derived
from the effective force inversely proportional to the distance

FIG. 2. Schematics of two reconnecting disclinations and theirdirector structures. (a) Two disclinations with the respective tanˆ 1, 
ˆ 2 are shown. The
gent vectors t̂1 , t̂2 and the rotation axes 
two disclinations separate a minimum distance δ and form an angle
ψ = arccos(t̂1 · t̂2 )  π /2 with the tangent vectors at the two closest
points. (b) Various director configurations of 3D disclination lines,
seen at a cross-section perpendicular to the line. All these are topologically equivalent, i.e., homeomorphic. The thick arrows indicate
ˆ a unit vector normal to the plane in which thethe rotation axis ,
director rotates by π along a closed loop around the disclination.

023284-2

APPROACH AND ROTATION OF RECONNECTING …

PHYSICAL REVIEW RESEARCH 6, 023284 (2024)

FIG. 3. Time dependence of the distance δ (a) and the angle
ψ (b) between two reconnecting disclinations. Results for all 12
reconnection events are shown with different colors.

δ [33] under the assumption that the force is balanced to a
drag force proportional to the disclination velocity. It was also
predicted by the disclination density tensor method [35] calculated from the nematic tensor Q, which contains the directorn̂ and the scalar order parameter, and observed in simulations[34]. From the viewpoint of universality, we note that a similar
scaling was also observed for quantum vortices in superfluid
helium 4 [9–13] and in the Gross-Pitaevskii equation [39].
For intersecting reconnections, the angle between the twodisclinations is another quantity characterizing them. We consider the two closest points of the reconnecting disclinations
1 and 2, with tangent vectors t̂1 and t̂2 for the respective
disclinations [Fig. 2(a)]. We define the interdisclination angle
ψ by ψ = arccos(t̂1 · t̂2 ), as shown in the inset of Fig. 2(a).
Since the sign of t̂1 and t̂2 is arbitrary, we choose it in such a
way that t̂1 · t̂2  0, i.e., 0  ψ  π /2. Our analysis showsthat the angle ψ is not constant during each reconnection
event but decreases over time [Fig. 3(b)]; in other words,
each disclination pair tends to be closer to parallel as timegoes on. This behavior is consistent with earlier numericalobservations [34,35] and theoretical prediction [35].
Let us quantitatively compare the time evolution of distance δ and angle ψ, which we obtained experimentally, withthe theory reported in Ref. [35]. According to the theory, thetime evolution of the distance δ and the angle ψ are describedby the following nondimensionalized simultaneous differential equations:dδˆ1·ˆ 2 ) cos ψ ,
= 4(
(2)
dt
δ
dψ
ˆ1·
ˆ 2 ) sin ψ .
= 4(
(3)
dt
δ2
ˆ 1 and ˆ 2 are the rotation axes, each
Here, the unit vectors 
normal to the plane in which the director rotates by π along
a closed loop around the disclination [Fig. 2(b)]. Dividing
Eq. (2) by Eq. (3) and integrating the result, we obtainsin ψ ∝ δ.

(4)This relation is tested with the experimental data in Fig. 4(a).
The result shows that, while δ and sin ψ are indeed related bya power law, the exponent value was estimated at 0.26 ± 0.08
by averaging over all events (the error indicates the standard
deviation), clearly smaller than 1, the prediction of Eq. (4).
This indicates a quantitative limitation of Eqs. (2) and (3).
To probe the dynamics of the disclinations, suppose
Eqs. (2) and (3) are replaced by the following more general

FIG. 4. Comparison with the theory. (a) Log-log plot of sin ψand δ. (b) Histogram of the coefficient ratio μ2 /μ1 . The averagevalue over events is μ2 /μ1 = 0.26 ± 0.08 (the error indicates thestandard deviation). (c), (d) Evaluation of each of the two equations,Eq. (8) (c) and Eq. (9) (d). Images from eight to two time frames
before the reconnection were used for all events. For (a), (c), and (d),results for all 12 reconnection events are shown with different colors.
(e) Histograms of the respective coefficients μ1 and μ2 . The values
are obtained by the fitting of the data in (c) and (d). (f) Scatter plotof the coefficient ratio obtained from the individual values of μ1 andμ2 in (e) against that directly obtained from the power law in (a).
Error bars indicate the 95% confidence intervals evaluated from therespective fittings.

form, with two different coefficients −μ1 , −μ2 < 0 insteadˆ1·ˆ 2 ) in Eqs. (2) and (3):of the common one, 4(
cos ψ
dδ
= −μ1
,
(5)
dt
δ
dψ
sin ψ
= −μ2 2 .
(6)dt
δThe coefficients may be thought of as effective mobilities forthe time evolution of δ and ψ. We remind that we chose the
ˆ1·
ˆ2 <
signs of t̂1 and t̂2 that satisfy cos ψ  0; as a result, 
0 is expected for an attracting disclination pair [33,34]. This
is why we set the negative coefficients −μ1 , −μ2 < 0. Then,similarly to the derivation of Eq. (4), by dividing Eq. (5) by
Eq. (6) and integrating the result, we obtainsin ψ ∝ δ μ2 /μ1 .

(7)This power law is to compare with the experimental data
in Fig. 4(a). The values of the power-law exponent, μ2 /μ1 ,023284-3

ZUSHI, SCHIMMING, AND TAKEUCHIPHYSICAL REVIEW RESEARCH 6, 023284 (2024)

varied a little among reconnection events [Fig. 4(b)], with theaverage being μ2 /μ1 = 0.26 ± 0.08 as already noted.
Equations (5) and (6) were also tested, respectively. Since
derivatives of experimental data are noisy, not Eqs. (5) and (6)themselves but the following integrated forms are considered: tcos ψ (t  )dt δ(ti ) − δ(t ) = μ1,
(8)δ(t  )
ti t
sin ψ (t  )
dt 
.
(9)
ψ (ti ) − ψ (t ) = μ2
δ(t  )2tiHere, ti is set to be the first frame used in the analysis.
By using images from eight to two timeframes before thereconnection, we indeed confirm that the left-hand side is proportional to the integral in the right-hand side for both Eqs. (8)and (9) [Figs. 4(c) and 4(d)]. This indicates that Eqs. (5) and(6) describe intersecting reconnections. Moreover, from the
proportionality coefficients, the individual values of μ1 andμ2 are obtained [Fig. 4(e)]. For each event, the coefficient
ratio is consistent with the value obtained directly from Eq. (7)[Fig. 4(f)].
To summarize the results so far, Eqs. (5) and (6) successfully describe the evolution of the distance and the angle ofintersecting reconnections, and the mobility ratio μ2 /μ1 =0.26 ± 0.08 was obtained. Physically, this mobility ratio represents the degree to which the angle between disclinations
changes with respect to the distance. The inequality μ2 /μ1 <1 indicates that the angle between disclinations changes relatively slower than the distance. However, the theory ofRef. [35] predicts Eqs. (2) and (3) with the common coefficient, hence μ2 /μ1 = 1, disagreeing with the experimental
result. What then determines the mobility ratio μ2 /μ1 in theexperiment?
IV. THEORY

To investigate the discrepancy in the mobility ratio μ2 /μ1
between the experimental data and the theory of Ref. [35], we
consider here the case of unequal nematic elastic constants.
Nematic distortions may be described by the elastic energy
densityK2
K3
K1|n̂ × (∇ × n̂)|2 , (10)(∇ · n̂)2 +(n̂ · (∇ × n̂))2 +2
22
where K1 , K2 , and K3 are the elastic constants for splay, twist,
and bend, respectively. In the theory derived in Ref. [35],it was assumed K1 , K2 , and K3 were equal; however, in theexperimental system this is not the case. MLC-2037, themesogen used in this experiment, has a smaller twist elastic
constant, K2 = 6.1 ± 0.5pN, than splay, K1 = 11.6pN, andbend, K3 = 13.2pN [17]. Therefore, to more closely align thetheory with experiments, we derive equations akin to Eqs. (2)and (3) but for a generalized case of K1 = K3 = K2 .
We first remark the director structure around a straightline disclination in its normal plane for the case K1 = K3 =K = K2 , discussed already in the literature, e.g., in Ref. [40].
For concreteness, we assume the tangent vector t̂ = ŷ, theˆ = ẑ, and the director around the disclination
rotation vector 
is given by n̂ = cos θ x̂ + sin θ ŷ, where θ = θ (x, z) is a yet
undetermined function. Substituting this into Eq. (10) andtaking a variational derivative yields the following equation for θ :
K

∂ 2θ
∂ 2θ
+ K2 2 = 0.
2∂x∂z

A solution that accommodates a defect is given byxθ (x, z) = m arctan √
,
εz

(11)

(12)

where m is a half-integer multiple winding number and ε =K/K2 is the ratio of elastic constants. We note that sinceEq. (11) is linear, the full director solution with multiple
defects is a superposition of Eq. (12) for each defect.
To derive the equations for δ and ψ, we use the methodsof Ref. [35] to approximate the structure of the nematic tensor
order parameter Q at the location of the disclination line. Thisapproximation is then used to approximate the velocity of thedisclination line. To do this, we assume that the dynamics of
the Q tensor is solely due to free-energy relaxation: ∂Q
1 δF T S
(13)
=−
∂t
γ δQ
where γ is a rotational viscosity, F is the free energy of thesystem, and [·]T S denotes the traceless, symmetric part of amatrix. The relevant elastic free energy density in terms of Q
isfe = L1 |∇Q|2 + L2 |∇ · Q|2 ,(14)

where L1 and L2 are elastic constants that may be mapped to
Ki :


K1
L2
=2
−1 .
(15)
L1
K2
We note that including the L2 term in Eq. (14) does not breakthe degeneracy K1 = K3 . Since we only consider relaxational
dynamics, the equations we derive do not take into accounthydrodynamic effects such as backflow. We assume the system is composed of two, infinitely long disclinations that
ˆ1·
ˆ2 =
remain straight with antiparallel rotation vectors −1 aligned along the z axis [Fig. 2(a)]. The resulting equations of motion for δ and ψ are (see Appendix C for details of
the calculation)(4L1 + 2L2 ) cos ψ
dδ
=−
,
(16)
γ
dt
εδ
dψ
4L1 sin ψ
γ
=−
.
(17)
dt
δ2
The predicted ratio of effective mobilities μ2 /μ1 is
μ2
2L1 ε
=
= 1,
μ1
2L1 + L2(18)

where Eq. (15) is used in the second equality. We thus predict the same mobility coefficient ratio as Eqs. (2) and (3),independent of the ratio of the elastic constants. We notethat this prediction does not take into account deformationsof the disclination line and the resulting director field. For a
deformed line disclination, Eq. (12) no longer holds, and amore complicated director configuration will be assumed by
the system.

023284-4

APPROACH AND ROTATION OF RECONNECTING …

PHYSICAL REVIEW RESEARCH 6, 023284 (2024)

of μ2 /μ1 . Since μ1 and μ2 correspond to mobilities for the
distance δ and the angle ψ, respectively, we can say that
disclinations with larger angles are harder to rotate to beparallel. This may suggest the effect of the deformation of
disclination lines, induced locally or globally by the rotationof disclinations, which is not considered in the theoretical
approach.
B. Simulations

FIG. 5. Scatter plots of the coefficient ratio μ2 /μ1 [as determined directly by Eq. (7)] against several properties of thedisclination pairs. (a) Against the time-averaged angle of the disclination pair θpair . The angle θpair is obtained by projecting the vectort̂1 + t̂2 to the xy plane and measuring the angle between the projected
vector and the easy axis (y), the direction of the surface alignment of
the cell. (b), (c) Against the speed vM (b) and the angle θdrift (c) of
the drift. The drift is a constant flow that the two disclinations are
exposed to, which is supposedly induced extrinsically (e.g., by otherdisclinations present outside the field of view). We evaluate the drift
velocity by using the midpoint of the two closest points. Then theangle θdrift is obtained by projecting the drift velocity in the xy planeand measuring the angle it makes with the easy axis. The drift isalso used to define the comoving frame to discuss the symmetryof dynamics in Appendix B. (d) Against the angle θdrift,pair formedby the drift velocity and the disclination pair. (e) Against the initial
angle ψi formed by the two disclinations, which is determined here
eight timeframes before the reconnection. (f) Against the curvatureof the disclinations. The values κ1 and κ2 are for the two respective
disclinations. The curvatures at the closest points right before thereconnection are used.
V. INFLUENCING FACTORS OF THE MOBILITYRATIO μ2 /μ1
A. Experiment

Since the reduced twist elastic constant alone turned out
to be unable to explain the small mobility ratio observed in
the experiment, we test if the experimental values of μ2 /μ1
for individual events [Fig. 4(b)] may be correlated with any
property of the disclination pairs (Fig. 5). Among the inspected properties, the only significant dependence we found
was a slightly negative correlation with the initial angle ψiformed by the two disclination lines [Fig. 5(e)]. This indicatesthat disclination pairs with larger angles have smaller values

We also probe the mobility ratio μ2 /μ1 by performing 3Dcomputations of the evolution of the Q tensor governed by
Eq. (13) for systems with a reduced twist elastic constant as
considered in the previous section. For the free energy, we
use the Maier-Saupe Ball-Majumdar potential for the bulk
free-energy density [41] and Eq. (14) for the elastic freeenergy density. Instead of applying Dirichlet or Neumannconditions, we directly solve the Euler-Lagrange equations for
Q at the boundary. The computations are carried out using
the MATLAB/C++ finite element package FELICITY [42],while matrix inversions are performed using AGMG [43–46].
We nondimensionalize the system by scaling lengths in termsof the nematic correlation length, and times in terms of thenematic relaxation time. We fix the computational time stepto t = 0.2. Further details involving the computational algorithm can be found in Ref. [47].
For the computations, we fix the ratio L2 /L1 = 2, corresponding to an elastic constant ratio of K2 /K1 = 0.5, similar
to that of the experiment. To consider a pair of straight line
disclinations that are initially apart by a distance δ0 in the z
direction and form an initial angle ψ0 in the xy plane aroundthe y axis [see also Fig. 2(a)], we initialize the system in a 3D
domain of dimensions [Lx , Ly , Lz ] with director fieldn̂ = cos(θ1 − θ2 )ŷ − sin(θ1 − θ2 )x̂,x cos ψ0 /2 − y sin ψ0 /2
1
θ1 (x, y, z) = arctan,√2
ε(δ0 /2 + z)1
−x cos ψ0 /2 − y sin ψ0 /2
θ2 (x, y, z) = arctan, (19)√2ε(δ0 /2 − z)which satisfies the Euler-Lagrange equation for the director,
Eq. (11). We fix the initial distance at δ0 = 4 and vary ψ0such that cos ψ0 ∈ {0.3, 0.4, 0.5, 0.6, 0.7, 0.8}. We computeon two different system sizes [Lx , Ly , Lz ] = [10, 10, 10] and[20,20,10] so, in the latter size, disclinations are initially twice
the length of those in the former. In Fig. 6, we show several
time snapshots of a simulated disclination reconnection where
the contours represent the locations of the defects.
We measure δ(t ) and ψ (t ) [Figs. 7(a) and 7(b)] and confirm the square-root law (1) for small |t − t0 |, as expected.
The agreement with the square-root law was better for larger
system sizes. Then we plot in Fig. 7(c) the log-log scalingof sin ψ versus δ for δ  2 and find a power-law relationshipbetween sin ψ and δ, just as observed in the experiments andpredicted by the theory [Eq. (7)]. For δ < 2, the disclination
cores begin to overlap and a clear power-law relationship was
no longer observed, which is reasonable because the estimation of ψ becomes less accurate and the theory is not expected
to hold for strongly overlapping disclination cores. Although
the square-root law (1) for δ(t ) is not yet established for δ  2,023284-5

ZUSHI, SCHIMMING, AND TAKEUCHIPHYSICAL REVIEW RESEARCH 6, 023284 (2024)

FIG. 6. Time snapshots of a simulated disclination reconnection
for initial angle cos ψ0 = 0.6, system size [Lx , Ly , Lz ] = [20, 20, 10],and time steps t/t = 5 (a), 50 (b), 88 (c), 97 (d). The contours
indicate points where the nematic scalar order parameter S = 0.3SN ,where SN is the equilibrium value of S in the nematic phase.

the smaller system size, μ2 /μ1 ∼ 0.9 and increases slightly
as ψ0 increases. The values are close to the theoretically predicted mobility ratio of μ2 /μ1 = 1, but are smaller, possibly
due to the deformations of disclinations as they reconnect. For
the larger system size, in which disclination lengths are doubled, we find that μ2 /μ1 is systematically smaller than that of
the smaller system size for each ψ0 . Further, μ2 /μ1 decreases
with increasing ψ0 for the larger system size, similarly to theexperiment [Fig. 5(e)].
The numerical results indicate that, for the parameters
explored here, the effective mobilities of disclination reconnections are sensitive to large-scale deformations of thedisclinations. The theoretically predicted μ2 /μ1 = 1 is mostclosely assumed by the numerical results at small systemsizes, when the disclinations cannot appreciably deform andremain better approximated by straight lines. When doubling
the length of disclinations, the mobility ratio decreases withincreasing angle, indicating that disclinations slow their rotation rate due to deformations along the disclination. These
deformations can be seen in Fig. 6 as the disclinations reconnect.
VI. CONCLUDING REMARKS

we note that it is not required for the power law (7) to hold,since theoretically Eq. (7) holds for all t as long as Eqs. (5)and (6) are valid, while to obtain Eq. (1) one needs to take thelimit |t − t0 | → 0.
Now, using the power law (7), we evaluate the effective
mobility ratio μ2 /μ1 by determining the exponent from eachcomputation result shown in Fig. 7(c). In Fig. 7(d), we plot
μ2 /μ1 as a function of initial angle ψ0 for the set of computations with system size Lx = Ly = 10 and Lx = Ly = 20. For

FIG. 7. Numerical results for reconnecting disclinations. (a),(b) Time dependence of the distance δ (a) and the angle ψ(b) between reconnecting disclinations for various simulations with
Lx , Ly = 10. (c) Log-log plot of sin ψ versus δ for various initial
angles ψ0 and two system sizes, Lx , Ly = 10 (circles) and Lx , Ly =20 (squares). Note that the data for Lx , Ly = 20 is vertically offset
to distinguish the data sets. (d) Mobility ratio μ2 /μ1 versus initial
disclination angle ψ0 for systems sizes Lx , Ly = 10 and 20. Lengths
and times are given in dimensionless simulation units.

We investigated intersecting reconnections of disclinations.
Experimentally, the square-root scaling was determined forthe time dependence of the distance δ between reconnecting
disclinations. The angle ψ between disclinations, which decreases over time, is also important to describe the dynamics.
It was found that Eqs. (2) and (3) predicted in Ref. [35] candescribe the experimentally observed time evolution of the
distance and the angle, except that the mobility coefficients,
μ1 and μ2 in Eqs. (5) and (6), were found to be different,
resulting in a value of the ratio μ2 /μ1 significantly smaller
than the theoretical prediction μ2 /μ1 = 1. We extended thetheory to the case of a reduced twist elastic constant andconsidered the change in the equilibrium director field, but thepredicted mobility ratio μ2 /μ1 turned out to remain the same.
We therefore searched for influencing factors of μ2 /μ1 in bothexperiments and simulations, and found that μ2 /μ1 tends todecrease with an increasing initial angle between the disclination pair. Numerically, the mobility ratio was also found to besmaller for longer disclination lines. Since disclinations may
not rotate while maintaining their straight shapes as assumed
in the theory, the rotation results in local deformation of thedisclinations, and this effect is stronger when two disclinations are longer or form a larger angle. We therefore consider
that the deformations of the disclinations may be relevant to
the reduced value of μ2 /μ1 .
Our results suggest a few interesting directions for futurestudies. First, it is important to develop a theoretical framework to deal with the dynamics of deformable disclinations.
It may also help to extend the analysis of experimental andnumerical data, to analyze not only the vicinity of the closest points but longer parts of the disclinations. The surfacealignment of the cell may also influence the dynamics of
disclinations. Second, even though the theory predicted thatthe mobility ratio μ2 /μ1 does not change by a reduced twist
elastic constant, previous work [35] has shown that numerical
data were closer to theoretical predictions, i.e., Eqs. (2) and

023284-6

APPROACH AND ROTATION OF RECONNECTING …

PHYSICAL REVIEW RESEARCH 6, 023284 (2024)(3), indicating μ2 /μ1 = 1, if equal elastic constants were assumed. This suggests that the reduced twist elastic constant
may still have a nontrivial effect, presumably affecting theway the deformation is involved. Therefore, it may be interesting to conduct experiments using liquid crystals with different
ratios of the elastic constants, K2 /K1 , by using, e.g., large K2
expected near the nematic-smectic transition [4] or for nematic discotic liquid crystals [48]. We also note that deviationof the director field from the equilibrium one [Eq. (12)] canchange the mobility ratio μ2 /μ1 . For example, if we replace
Eq. (12) by that for the one-constant case (ε = 1), we obtain
μ2 /μ1 = K2 /K1 < 1 in the theory. This suggests the potential
importance of observing the director field around disclinations, by methods such as the fluorescence confocal polarizing
microscopy [49,50], two- or three-photon excitation fluorescence polarizing microscopy [51–53], and the tomographic
measurement of the dielectric tensor [54].
As described in the Introduction, since topological defects
provide useful means to control microparticles and light in liquid crystal medium, better understanding of defect dynamics
can contribute to developments in this direction. Moreover,as topological defects also appear in various scientific fields
other than liquid crystals, it is also important to unravel thegeneral behavior of topological defects beyond liquid crystals.
We hope our paper will contribute to these and trigger further
investigations to elucidate fully 3D dynamics of topologicaldefects.

ACKNOWLEDGMENTS

We thank S. Shankar for his suggestion to use Eq. (7)to compare the experimental data with the theory. We acknowledge the material data of MLC-2037 provided by Merckand their permission to present them. This paper is supported in part by Japan Science and Technology Agency(JST) Precursory Research for Embryonic Science and Technology (PRESTO) (Grant No. JPMJPR18L6), by KAKENHI
from Japan Society for the Promotion of Science (Grants
No. JP22J12144, No. JP22KJ0843, No. JP19H05800, No.
JP19H05144, No. JP20H01826, No. JP23K17664, and No.
JP24K00593) by JSR Fellowship (The University of Tokyo),and by FoPM, WINGS Program (The University of Tokyo).
C.D.S. acknowledges support from the U.S. Department ofEnergy through the Los Alamos National Laboratory.
K.A.T. designed research. Y.Z. performed experiment. Y.Z.
and K.A.T. analyzed data. C.D.S did theoretical calculationand numerical computations. Y.Z., C.D.S., and K.A.T. wrote
the paper.

APPENDIX A: IMAGE ANALYSISThe 3D coordinates of disclinations were extracted fromthe obtained data of fluorescence intensity at 3D positions(x, y, z) via a method called snakes [38]. It is a way to findsmooth contours. Since we know disclinations are smoothlines, the method is suitable for extracting disclinations’ positions. In the snake method, the position of a disclination line
is represented by v (s) = (x(s), y(s), z(s)) with a parameter s,and its shape is determined to minimize the total cost function,or energy, given byE ∗ = E(v (s))ds
= (wint Eint (v (s)) + wim Eim (v (s)))ds.

(A1)

Here, Eint is the internal energy, Eim is the image energy, andwint and wim are the weights for the internal energy and theimage energy, respectively. The internal energy is given by
 


 2
 ∂v (s) 2
 ∂ v(s) 2
1

 ,


Eint =+ βint (A2)αint 
2∂s ∂s where αint and βint are coefficients. The first and the second
terms of Eq. (A2) correspond to the energy cost due to thelength and the roughness of the line, respectively. For theimage energy, we adopted the intensity I (x, y, z) itself here,
since we assume that the fluorescent intensity is higher nearer
to the disclination core:
Eim = −I (x, y, z).

(A3)

Practically, the contour v was expressed as a series of points,
with intervals of approximately 2 µm, and the minimizationof Eq. (A1) was implemented by the gradient descent method.
Since this implementation can, in principle, only reach local
minima, it is important to start from an initial condition close
to the desired result. For that purpose, disclination shapes
roughly estimated from the obtained images were used as
initial conditions for the first timeframe. From the second
timeframe, the result of the previous timeframe was used as
the initial condition for the subsequent frame. To reduce noise
and improve convergence, the images were filtered with a3 × 3 × 3 median filter and a Gaussian filter (a kernel witha standard deviation of 0.5). Concerning the time stamp, although the fluorescence intensity at different positions wasrecorded at slightly different times in the laser scanning confocal microscopy, we used a unique time value for each 3D
image, represented by the time at the z coordinate of the
midpoint of the two closest points. The reconnection time t0
was determined not from confocal images but from the 2D
image from the transmitted excitation laser to benefit from thefiner time resolution.
APPENDIX B: SYMMETRY

We discussed the symmetry of defect dynamics in ourprevious work [17]. In 2D, a pair of pointlike defects withwinding numbers of ±1/2 is well-known to approach asymmetrically and annihilate, but 3D line-shaped defects were
found to approach in a symmetric manner and reconnect, atleast for the in-plane reconnections analyzed in Ref. [17]. Thissymmetry restoring is explained on the basis of the topological equivalence of various structures of 3D disclinations andthe lower energy of twist defects, which result in symmetricdynamics.
Here we inspect the symmetry of the dynamics for the caseof intersecting reconnections. As discussed in Ref. [17], thedisclination dynamics in the laboratory frame may include
extrinsic effects due to other disclinations that exist outside
the field of view. Such effects appear as a drift that is uniform

023284-7

ZUSHI, SCHIMMING, AND TAKEUCHIPHYSICAL REVIEW RESEARCH 6, 023284 (2024)

comprises two infinitely long disclinations at distance δ andangle ψ with respect to one another. We fix a coordinate
system so one disclination is along the y axis, while the otherdisclination lies in the xy plane, and that the shortest points
between disclinations are along the z axis. We also assume theˆ 2 = ẑ, as sketched in Fig. 2(a).
ˆ 1 = −rotation vectors are 
Finally, we assume the director structure around a single isolated disclination is given byˆn̂ = cos θ x̂ + sin θ y,
FIG. 8. Histogram of the asymmetry parameter A.

in the field of view and has a constant velocity for each reconnection event. Therefore, these extrinsic effects were removedby evaluating the dynamics in the frame comoving with theconstant drift. The comoving frame was defined for eachreconnection event in a manner similar to that adopted forin-plane reconnections [17]. The midpoint of the two closest
points of the respective disclination lines was located at eachtime, and the time evolution of the midpoint was fitted by alinear function of time. The slope of this function determines
the constant velocity of the comoving frame.
In the comoving frames, two reconnecting disclinations
turn out to approach each other symmetrically for all intersecting reconnections (Fig. 8). Therefore, the symmetry-restoringmechanism discussed in Ref. [17] also applies to intersecting
reconnections, and disclinations are considered to take a twist
configuration during the events. Asymmetry parameter A wasdetermined for each reconnection event as follows. First, we
obtained the distance D1 (t ) and D2 (t ) between respective
disclinations and the reconnection point X 0 , and fit Di2 (t ) by
the equation Di2 (t )  Ci2 (t0 − t ) with constants Ci . Then, A isdefined by
A :=max{C1 , C2 }.
min{C1 , C2 }(B1)Here, the estimate of the reconnection point X 0 was refinedby using the time-dependent coordinates of the disclinations
before the reconnection. Specifically, X 0 = (X0 , Y0 , Z0 ) was
determined in such a way that the scaling Di (t )  Ci |t −t0 |1/2 is satisfied more precisely in a time period before thereconnection. This was done by evaluating Di (t ) with thereconnection point varied over six neighboring positions in
3D, fitting it to Di (t )2 = ai |t − t0 | + bi , choosing the directionthat minimizes b21 + b22 , and iterating this to reach the (local)
minimum. The point was moved at an interval of 1/5 ofthe voxel size. The final results were at most 2.5 µm awayfrom the first rough estimate, which was located from theseries of transmitted and confocal images (X0 and Y0 from thetransmitted images, Z0 from the confocal images).
APPENDIX C: DERIVATION OF THEORETICAL
EQUATIONS OF MOTION

Here we give more details leading to the equations of
motion of Eqs. (16) and (17). As mentioned above, we assume the dynamics are purely relaxational and that the system(C1)

where θ is given by Eq. (12).
Following the methods of Ref. [35], we use the kinematic
equation for the velocity of a disclination line,ˆ · g t̂ × 
v=2
(C2)
 ,

|D|
r=R

where t̂ is the tangent vector of the disclination, gγ k =εγ μν ∂t Qμα ∂k Qνα is related to the topological current, Dγ i =εγ μν εik ∂k Qμα ∂ Qνα , and we have assumed summation on
repeated indices and simplified notation so that ∂k ≡ ∂/∂xk .
Note that these quantities need only be computed at the location of the disclination core, r = R.
The Q-tensor dynamics areγ ∂t Qμα = 2L1 ∂k ∂k Qμα2+ L2 ∂α ∂k Qμk + ∂μ ∂k Qαk − δμα ∂ ∂k Qk .
3(C3)
Upon substituting Eq. (C3) into Eq. (C2), we find that thelast term of Eq. (C3) does not contribute to disclination
motion, while the first two L2 terms give the same contribution, so it suffices to simplify and consider only γ ∂t Qμα =2L1 ∂k ∂k Qμα + 2L2 ∂α ∂k Qμk . Given the new director structureof Eq. (12), the linear core approximation for Q is
√

1
1
εz
(x̂ ⊗ x̂ − ŷ ⊗ ŷ)Q ≈ SN I − ẑ ⊗ ẑ +
6
2
2a

x
+ (x̂ ⊗ ŷ + ŷ ⊗ x̂) ,(C4)2awhere SN is the magnitude of Q in the nematic phase and ais the disclination core radius. Substituting this into Eqs. (C3)and (C2) and taking into account the rotation of the director
caused by the other disclination line, the velocity of the first
disclination line along the y axis isv 1 (y) =2L1 yε sin ψ(2L1 + L2 )δ cos ψẑ − 2x̂.
222εδ + y sin ψεδ + y2 sin2 ψ(C5)

A similar equation may be derived for the velocity of theother disclination, v2 . To obtain equations of motion for δ andψ, we use the identitiesdδ(C6)= ẑ · (v2 (0) − v1 (0)),
dtdv1 
dψ
dv2 
· t̂1 +· t̂2 ,
(C7)
− sin ψ=dt
dy y=0dy y=0
leading to Eqs. (16) and (17).

023284-8

APPROACH AND ROTATION OF RECONNECTING …

PHYSICAL REVIEW RESEARCH 6, 023284 (2024)[1] P. M. Chaikin and T. C. Lubensky, Principles of Condensed
Matter Physics (Cambridge University Press, Cambridge,2000).
[2] M. Nakahara, Geometry, Topology and Physics, 2nd ed. (Taylor& Francis, Boca Raton, 2003).
[3] S. Chandrasekhar, Liquid Crystals, 2nd ed. (Cambridge University Press, Cambridge, 1992).
[4] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, 2nd ed., International Series of Monographs on Physics,
Vol. 83 (Oxford University Press, New York, 1995).
[5] W. H. Zurek, Cosmological experiments in condensed mattersystems, Phys. Rep. 276, 177 (1996).
[6] A. Vilenkin and E. P. S. Shellard, Cosmic Strings and Other
Topological Defects, revised ed. (Cambridge University Press,
Cambridge, 2000).
[7] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F.
Sagués, Active nematics, Nat. Commun. 9, 3246 (2018).
[8] A. Doostmohammadi and B. Ladoux, Physics of liquid crystals
in cell biology, Trends Cell Biol. 32, 140 (2022).
[9] G. P. Bewley, M. S. Paoletti, K. R. Sreenivasan, and D. P.
Lathrop, Characterization of reconnecting vortices in superfluid
helium, Proc. Natl. Acad. Sci. USA 105, 13707 (2008).
[10] M. S. Paoletti, M. E. Fisher, K. R. Sreenivasan, and D. P.
Lathrop, Velocity statistics distinguish quantum turbulencefrom classical turbulence, Phys. Rev. Lett. 101, 154501 (2008).
[11] M. S. Paoletti, M. E. Fisher, and D. P. Lathrop, Reconnection
dynamics for quantized vortices, Physica D 239, 1367 (2010).
[12] E. Fonda, D. P. Meichle, N. T. Ouellette, S. Hormoz, and
D. P. Lathrop, Direct observation of Kelvin waves excited byquantized vortex reconnection, Proc. Natl. Acad. Sci. USA 111,
4707 (2014).
[13] Y. Minowa, S. Aoyagi, S. Inui, T. Nakagawa, G. Asaka, M.
Tsubota, and M. Ashida, Visualization of quantized vortex reconnection enabled by laser ablation, Sci. Adv. 8, eabn1143(2022).
[14] S. Serafini, M. Barbiero, M. Debortoli, S. Donadello, F.
Larcher, F. Dalfovo, G. Lamporesi, and G. Ferrari, Dynamics
and interaction of vortex lines in an elongated Bose-Einsteincondensate, Phys. Rev. Lett. 115, 170402 (2015).
[15] S. Serafini, L. Galantucci, E. Iseni, T. Bienaimé, R. N. Bisset,
C. F. Barenghi, F. Dalfovo, G. Lamporesi, and G. Ferrari, Vortex
reconnections and rebounds in trapped atomic Bose-Einsteincondensates, Phys. Rev. X 7, 021031 (2017).
[16] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Cosmology inthe laboratory: Defect dynamics in liquid crystals, Science 251,
1336 (1991).
[17] Y. Zushi and K. A. Takeuchi, Scaling and spontaneous symmetry restoring of topological defect dynamics in liquid crystal,
Proc. Natl. Acad. Sci. USA 119, e2207349119 (2022).
[18] I. Muševič, Nematic liquid-crystal colloids, Materials 11, 24(2018).
[19] I. I. Smalyukh, B. I. Senyuk, S. V. Shiyanovskii, O. D.
Lavrentovich, A. N. Kuzmin, A. V. Kachynski, and P. N. Prasad,
Optical trapping, manipulation, and 3D imaging of disclinations
in liquid crystals and measurement of their line tension, Mol.
Cryst. Liq. Cryst. 450, 79 (2006).
[20] J. Jiang, X. Wang, O. I. Akomolafe, W. Tang, Z. Asilehan,
K. Ranabhat, R. Zhang, and C. Peng, Collective transport andreconfigurable assembly of nematic colloids by light-driven

cooperative molecular reorientations, Proc. Natl. Acad. Sci.
USA 120, e2221718120 (2023).
[21] X. Wang, J. Jiang, J. Chen, Z. Asilehan, W. Tang, C. Peng, andR. Zhang, Moiré effect enables versatile design of topologicaldefects in nematic liquid crystals, Nat. Commun. 15, 1655
(2024).
[22] C. Meng, J.-S. Wu, and I. I. Smalyukh, Topological steering oflight by nematic vortices and analogy to cosmic strings, Nat.
Mater. 22, 64 (2023).
[23] X. Wang, D. S. Miller, E. Bukusoglu, J. J. de Pablo, and N. L.
Abbott, Topological defects in liquid crystals as templates formolecular self-assembly, Nat. Mater. 15, 106 (2016).
[24] T. Ohzono, K. Katoh, and J.-I. Fukuda, Fluorescence microscopy reveals molecular localisation at line defects innematic liquid crystals, Sci. Rep. 6, 36477 (2016).
[25] B. S. Murray, R. A. Pelcovits, and C. Rosenblatt, Creating
arbitrary arrays of two-dimensional topological defects, Phys.
Rev. E 90, 052501 (2014).
[26] H. Yoshida, K. Asakura, J. Fukuda, and M. Ozaki, Threedimensional positioning and control of colloidal objects
utilizing engineered liquid crystalline defect networks, Nat.
Commun. 6, 7180 (2015).
[27] J. Jiang, K. Ranabhat, X. Wang, H. Rich, R. Zhang, and C. Peng,
Active transformations of topological structures in light-driven
nematic disclination networks, Proc. Natl. Acad. Sci. USA 119,e2122226119 (2022).
[28] C. Long, M. J. Deutsch, J. Angelo, C. Culbreath, H. Yokoyama,J. V. Selinger, and R. L. B. Selinger, Frank-read mechanism in nematic liquid crystals, Phys. Rev. X 14, 011044(2024).
[29] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S.
Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, andB. Ladoux, Topological defects in epithelia govern cell deathand extrusion, Nature (London) 544, 212 (2017).
[30] K. Copenhagen, R. Alert, N. S. Wingreen, and J. W. Shaevitz,
Topological defects promote layer formation in Myxococcus
Xanthus colonies, Nat. Phys. 17, 211 (2021).
[31] T. Shimaya and K. A. Takeuchi, Tilt-induced polar order andtopological defects in growing bacterial populations, PNAS
Nexus 1, pgac269 (2022).
[32] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E.
Braun, and K. Keren, Topological defects in the nematic orderof actin fibres as organization centres of Hydra morphogenesis,Nat. Phys. 17, 251 (2021).
[33] C. Long, X. Tang, R. L. Selinger, and J. V. Selinger, Geometry
and mechanics of disclination lines in 3D nematic liquid crystals, Soft Matter 17, 2265 (2021).
[34] C. D. Schimming and J. Viñals, Singularity identification for thecharacterization of topology, geometry, and motion of nematic
disclination lines, Soft Matter 18, 2234 (2022).
[35] C. D. Schimming and J. Viñals, Kinematics and dynamics of
disclination lines in three-dimensional nematics, Proc. R. Soc.
A 479, 20230042 (2023).
[36] See Supplemental Material at http://link.aps.org/supplemental/
10.1103/PhysRevResearch.6.023284 for videos 1–4, whichshow fluorescence observation of intersecting reconnections.
Videos 1 and 2 show the event displayed in Fig. 1(c), fromthe top and side, respectively. Videos 3 and 4 show anotherreconnection event, again from the top and side, respectively.

023284-9

ZUSHI, SCHIMMING, AND TAKEUCHIPHYSICAL REVIEW RESEARCH 6, 023284 (2024)[37] S. Kai and W. Zimmermann, Pattern dynamics in the electrohydrodynamics of nematic liquid crystals, Prog. Theor. Phys.
Suppl. 99, 458 (1989).
[38] M. Kass, A. Witkin, and D. Terzopoulos, Snakes: Active contour models, Int. J. Comput. Vis. 1, 321 (1988).
[39] A. Villois, D. Proment, and G. Krstulovic, Universal andnonuniversal aspects of vortex reconnections in superfluids,
Phys. Rev. Fluids 2, 044701 (2017).
[40] J. Geurst, A. Spruijt, and C. Gerritsma, Dynamics of s = 1/2
disclinations in twisted nematics, J. Phys. France 36, 653(1975).
[41] J. M. Ball and A. Majumdar, Nematic liquid crystals: From
Maier-Saupe to a continuum theory, Mol. Cryst. Liq. Cryst. 525,
1 (2010).
[42] S. W. Walker, Felicity: A Matlab/C + + toolbox for developing
finite element methods and simulation modeling, SIAM J. Sci.
Comput. 40, C234 (2018).
[43] Y. Notay, An aggregation-based algebraic multigrid method,
Electron. Trans. Numer. Anal. 37, 123 (2010).
[44] A. Napov and Y. Notay, Algebraic analysis of aggregationbased multigrid, Numer. Linear Algebra Appl. 18, 539(2011).
[45] A. Napov and Y. Notay, An algebraic multigrid method withguaranteed convergence rate, SIAM J. Sci. Comput. 34, A1079(2012).
[46] Y. Notay, Aggregation-based algebraic multigrid forconvection-diffusion equations, SIAM J. Sci. Comput. 34,
A2288 (2012).

[47] C. D. Schimming, J. Viñals, and S. W. Walker, Numericalmethod for the equilibrium configurations of a Maier-Saupebulk potential in a Q-tensor model of an anisotropic nematic
liquid crystal, J. Comput. Phys. 441, 110441 (2021).
[48] M. Osipov and S. Hess, The elastic constants of nematic andnematic discotic liquid crystals with perfect local orientational
order, Mol. Phys. 78, 1191 (1993).
[49] I. Smalyukh, S. Shiyanovskii, and O. Lavrentovich, Threedimensional imaging of orientational order by fluorescence
confocal polarizing microscopy, Chem. Phys. Lett. 336, 88(2001).
[50] O. D. Lavrentovich, Fluorescence confocal polarizing microscopy: Three-dimensional imaging of the director, Pramana
J. Phys. 61, 373 (2003).
[51] T. Lee, R. P. Trivedi, and I. I. Smalyukh, Multimodal nonlinear
optical polarizing microscopy of long-range molecular order inliquid crystals, Opt. Lett. 35, 3447 (2010).
[52] R. P. Trivedi, I. I. Klevets, B. Senyuk, T. Lee, and I. I. Smalyukh,
Reconfigurable interactions and three-dimensional patterningof colloidal particles and defects in lamellar soft media, Proc.
Natl. Acad. Sci. USA 109, 4744 (2012).
[53] P. J. Ackerman and I. I. Smalyukh, Reversal of helicoidal twist
handedness near point defects of confined chiral liquid crystals,
Phys. Rev. E 93, 052702 (2016).
[54] S. Shin, J. Eun, S. S. Lee, C. Lee, H. Hugonnet, D. K. Yoon,
S.-H. Kim, J. Jeong, and Y. Park, Tomographic measurementof dielectric tensors at optical frequency, Nat. Mater. 21, 317(2022).

023284-10