Liquid crystal defect structures with Möbius strip topology

Topological solitons commonly appear as energy-minimizing field configurations, but examples of stable, spatially localized objects with coexisting solitonic structures and singular defects are rare. Here we use a nonpolar chiral liquid crystal system to show how twist domain walls can co-self-assemble with vortices to form spatially localized topological objects with spontaneous folding. These soliton–vortex assemblies, which we call ‘möbiusons’, have a topology of the molecular alignment field resembling that of the Möbius strip’s surface and package localized field excitations into folded structures within a confinement-frustrated uniform far-field background. Upon supplying energy in the form of electric pulses, möbiusons with different overall symmetries of structure exhibit folding-dependent rotational and translational motions, as well as topological cargo-carrying abilities that can be controlled by tuning the amplitude and frequency of the applied fields. We demonstrate on-demand transformations between various möbiusons and show examples of encoding information by manipulating folds in such structures. A model based on the energetics of solitons and vortices provides insights into the origins of the folding instability, whereas minimization of the Landau–de Gennes free energy closely reproduces details of their internal structure. Our findings may provide a route towards topology-enabled light-steering designs. Topological defect structures that swim have been realized in liquid crystals. Now, a range of structures with topology reminiscent of a Möbius strip swim and transform into one another.


Liquid crystal defect structures with Möbius strip topology
Hanqing Zhao 1 , Jung-Shen B. Tai  Topological solitons commonly appear as energy-minimizing field configurations, but examples of stable, spatially localized objects with coexisting solitonic structures and singular defects are rare. Here we use a nonpolar chiral liquid crystal system to show how twist domain walls can co-self-assemble with vortices to form spatially localized topological objects with spontaneous folding. These soliton-vortex assemblies, which we call 'möbiusons', have a topology of the molecular alignment field resembling that of the Möbius strip's surface and package localized field excitations into folded structures within a confinement-frustrated uniform far-field background. Upon supplying energy in the form of electric pulses, möbiusons with different overall symmetries of structure exhibit folding-dependent rotational and translational motions, as well as topological cargo-carrying abilities that can be controlled by tuning the amplitude and frequency of the applied fields. We demonstrate on-demand transformations between various möbiusons and show examples of encoding information by manipulating folds in such structures. A model based on the energetics of solitons and vortices provides insights into the origins of the folding instability, whereas minimization of the Landau-de Gennes free energy closely reproduces details of their internal structure. Our findings may provide a route towards topology-enabled light-steering designs.
Topological solitons and singular defects can enable diverse physical behaviour and technological utility, ranging from dislocation-mediated malleability of metals 1 to Skyrme-Witten descriptions of elementary particles as three-dimensional solitons [2][3][4] and models of cosmic string emergence after the Big Bang 5 . In liquid crystals (LCs), singular topological defects, which contain spatial regions where molecular ordering directions cannot be defined, typically appear as transient structures 1,2,[5][6][7][8][9][10][11] , although they can be stabilized by colloidal inclusions, patterned boundary conditions and confinement [12][13][14][15][16][17][18] . Solitonic structures in chiral LCs can also stabilize singular point and line defects (vortices) against shrinkage or annihilation 2,19,20 . However, the morphology of known field configurations containing such structures is limited to straight-line or looped configurations of twist walls and skyrmions coexisting with singular defects 2 . In this work, we show how solitonic twist walls in a chiral LC can co-self-assemble with vortices and then spontaneously fold into diverse spatially localized configurations. The nonpolar director field n(r), which describes the average local rod-like LC molecular alignment direction, exhibits highly elaborate, three-dimensional (3D) structures that comprise both solitonic and singular regions with the non-orientable n(r) tangent to a Möbius strip's surface. These spatially localized excitations, which we call 'möbiusons', spontaneously fold into diverse structures with different symmetries. We reconstruct and numerically model the 3D n(r) by using nonlinear optical microscopy Article https://doi.org/10.1038/s41567-022-01851-1 back to the mesophase ( Fig. 1c and Supplementary Video 1), a process that allows for embedding topologically nontrivial n(r) configurations. When a voltage of U ≈ 3 V is applied across the LC, this fibre-like twisted region spontaneously curves and folds ( Fig. 1d and Supplementary Video 1) until eventually relaxing into various stable möbiuson configurations (Fig. 1e-i) that emerge depending on the originally straight fibre's length. These emergent möbiusons are classified based on the number of folds N (Fig. 1e-i) and the symmetry of their polarizing optical micrographs (POMs), where such structures with mirror symmetric and right-or left-handed folding are denoted by N M , N R and N L , respectively. Micrographs of N M möbiusons contain a single mirror symmetry plane, characterized by the two-dimensional (2D) symmetry group D 1 (Fig. 1e). Each micrograph of N R or N L structure contains a two-fold rotation axis (2D symmetry group C 2 ), whereas N R and N L are mirror images of each other (Fig. 1f,g). Even-N möbiusons are either N R or N L , with their ends folding counter-clockwise or clockwise, respectively, whereas odd-N möbiusons are N M . Among special cases (Fig. 1), möbiusons with N = 2 exhibit all three types (2 M , 2 R and 2 L ), and the ones with N = 1 are torons with D ∞ symmetry ( Fig. 1h) 2 , denoted by 1 ∞ . Möbiusons can and by minimizing the Landau-de Gennes free energy 1,10-14 , respectively, showing how möbiusons emerge from competing free energy contributions. Möbiusons display inter-transformations, symmetry-dependent electrically powered rotational and translational motions and cargo-carrying abilities that exceed the known effects of out-of-equilibrium dynamics of defects and solitons [6][7][8][9] .

Emergence and inter-transformation of möbiusons
In a glass cell with perpendicular boundary conditions, confinement hinders chiral LC's tendency to twist, causing frustration (Fig. 1a). Möbiusons embed localized fibre-like regions of twisted n(r), which spontaneously fold and locally relieve the frustration (Fig. 1a). Across these fibre-like regions, n(r) twists by 180° while being tangent to a Möbius strip (as detailed below), the simplest non-orientable surface with only one side (Fig. 1b). To create a möbiuson, we first apply a voltage of U ≈ 3.5 V and generate an extended fragment of a localized fibre-like structure with 180°-twisted n(r) across its width by local melting of the chiral LC using laser tweezers and then quenching and its Möbius-like building blocks (top right), within which the director field is tangent to the surface of a Möbius strip. Here, n(r) is coloured according to orientations as shown in the colour-coded order parameter space (bottom right). b, A Möbius strip made from paper. c,d, POMs of a linear (c) and a spontaneously curved (d) fibre-like structure. e-h, POMs (left) and bright-field images (right) of möbiusons with mirror symmetry (e), right-(f) and left-handed folding symmetry (g), and D ∞ symmetry (h). i, Experimental encoding of 'C' and 'U' using the alphabetic ASCII codes, where the number of folds of the möbiuson is correlated to the ASCII table. Cell gap d = 10 μm, pitch p = 4.6 μm (d/p = 2.2) and U ≈ 3.0 V for all images. The directions of crossed polarizers in POMs are shown by double arrows. Scale bar 30 μm in a, 5 μm in c-h and 30 μm in i.

Article
https://doi.org/10.1038/s41567-022-01851-1 inter-transform through the creation and annihilation of folds, preserving or changing 2D symmetries. For example, upon increasing the voltage by ~0.02 V beyond the initial equilibrium value, a möbiuson can transform into its lower-N counterpart (Supplementary Video 2). This fold annihilation process reduces N by 2 with the following selection rules: mirror-symmetric möbiusons retain symmetry and transform as N M → (N − 2) M , whereas left-or right-handed ones change handedness, that is, N R → (N − 2) L and N L → (N − 2) R . However, low-N möbiusons with N ≤ 3 are exceptions, where 3 M →2 M and N = 2 of all symmetries can transform to 1 ∞ . These transformations are summarized in Extended Data Fig. 1. To add extra orders of folding to an existing localized excitation, one can grab one end of the initial möbiuson with laser tweezers and stretch it slowly. After turning off the laser, new folds emerge as the elongated möbiuson evolves towards equilibrium (Supplementary Video 3) by following the transformation rules in a reverse order. Stable möbiusons can be used to encode information, as shown in Fig. 1i for the example of 'CU'.

Three-dimensional structures of möbiusons
To understand the mechanism behind the formation of möbiusons, we investigate their internal structures. Before folding, the fibre-like translationally invariant feature (referred to as a cholesteric finger 21,22 ; Fig. 1c), comprises a twist wall and two vortex lines (twist disclinations) at the top and bottom substrates, which turn at the finger's ends and close into a single loop (Fig. 2a-c and Extended Data Fig. 2a). This nonsingular wall is a one-dimensional soliton embedded into the background n 0 of the unwound director field, where one-point compactification of its physical space 2,23 connects the equivalent director orientations after the 180° twist, much like in the case of making a Möbius strip from paper (Fig. 1b). A loop of the vortex line embeds the twist wall into the n 0 background, as seen within the cross-sections (Fig. 2a,b), where both the soliton and the vortex are labelled as elements of the first homotopy group, π 1 ( 2 /ℤ 2 ). The LC order is reduced within the vortex cores and the free energy density is much higher as compared with bulk regions 1,14 , whereas the twist wall is energetically favourable due to the LC's intrinsic chirality ( Fig. 2a and Extended Data Fig. 2b,c). The electric field applied along n 0 tends to constrain the extent of the twisted region by dielectrically coupling to n(r). The energy-minimizing shrinking tendency of the vortex line and expanding tendency of the twist wall, both tuneable by applied voltage, thus render the linear configuration of co-assembled soliton and vortex lines unstable, driving them to curve and fold (Figs. 1c-i and 2c-l). Optical micrographs obtained at different depths of focus reveal that the twist wall takes a longer spatial trajectory while the vortex line adopts a shorter one ( Fig. 2j and Supplementary Video 4). In high-N möbiusons (Figs. 1i and 2j-l), this leads to multiple folded regions, further lowering the total free energy by optimizing the geometric configurations of walls and singular vortex lines. The formation of möbiusons represents a symmetry-breaking process of the initial linear configuration with D 2 symmetry of the 2D optical image into structures with D 1 or C 2 symmetry. A special case emerges when the initial fibre-like finger fragment is short, so that the two ends interact while curving, eventually merging to yield the 1 ∞ toron (Fig. 2c-i and Supplementary Video 5). We reproduce this transformation in numerical simulations using Landau-de Gennes free energy minimization 1,10-14 and confirm that, when the two ends meet, the total free energy is reduced by merging and annihilating the two vertical segments of the high-energy π 1 ( 2 /ℤ 2 ) = ℤ 2 vortex line ( Fig. 2c-f). The process of merging the finger's ends transforms a single vortex loop into two loops upon reconnection, which then shrink into nanoscale loops identified as π 2 ( 2 /ℤ 2 ) point defects 2,24 . In the horizontal midplane, the ensuing solitonic structure is the π 2 ( 2 /ℤ 2 ) skyrmion, which emerges from looping and reconnecting the ends of a π 1 ( 2 /ℤ 2 ) twist wall (Extended Data Fig. 3). Thus, the overall topological transformation contains reconnection of a single π 1 ( 2 /ℤ 2 ) twist vortex loop of zero hedgehog charge into two loops after the finger's ends merge, which then become two π 2 ( 2 /ℤ 2 ) point defects of opposite relative hedgehog charge (±1), as well as the transformation of a π 1 ( 2 /ℤ 2 ) twist wall fragment into a π 2 ( 2 /ℤ 2 ) skyrmion tube terminating at the π 2 ( 2 /ℤ 2 ) point defects. Remarkably, this transformation between topological objects of different homotopy classes 23 is driven by free energy minimization (Fig. 2f) and involves embedding the solitonic structures into the background n 0 by singular defects of the same homotopy group, both before and after the transformation.
To further probe the 3D structures of möbiusons, we utilize three-photon excitation fluorescence polarizing microscopy (3PEF-PM) 25 synergistically combined with the structural analysis based on minimizing the free energy [10][11][12][13][14] . The 3PEF-PM imaging with different polarizations of the excitation laser light allows for the reconstruction of n(r), which then serves as initial conditions for the Landau-de Gennes free energy minimization, yielding the equilibrium 3D structures used to simulate 3PEF-PM images (Methods). The agreement between the computer-simulated and experimental images validates the n(r) reconstruction (  Fig. 6b) exhibits the two-fold rotation symmetry with respect to the vertical axis along its centre but is also asymmetric with respect to the cell's midplane. The handedness of folding, determined when viewed along the cell substrate's normal, the +z direction, correlates with the vortex loop symmetry, so that the 3D n(r) of 2 R and 2 L and all the right-left versions of möbiusons can be obtained from one another by flipping the LC sample upside down ( Fig. 3 and Extended Data Figs. 7 and 8). Although POM images of 2 R and 2 L distantly resemble twistions 19,20 , 3D imaging reveals different topology as these möbiusons comprise π 1 ( 2 /ℤ 2 ) = ℤ 2 vortex lines whereas twistions contain π 2 ( 2 /ℤ 2 ) point defects. The up-down asymmetry of high-N möbiusons manifests itself at the opposite ends of these structures, although they are symmetric with respect to the cell midplane within the central regions comprising many repeated folds (Extended Data Fig. 8b). For möbiusons with many folds, lateral displacement amplitudes of the twist wall D t and the vortex line D v from the average centre line and the periodicity of folding L vary with temperature, and the voltage U and the periodicity T U of the oscillating electric field (Fig. 4). By minimizing the free energy as a function of L, we find that the folded structure has lower free energy than the unwound state or the straight finger, with the equilibrium configuration featuring L/D t ≈ 1.76, close to the experimentally measured ratio of L/D t ≈ 1.84 (Fig. 4c).

Folding instability within the möbiuson
To gain insights into the möbiuson formation, we constructed an analytical model based on the minimization of the free energy including the vortex lines and the twist wall. The energy density of the vortex line F v contains contributions from its singular isotropic core (F c ) and elastic deformations of n(r) around it (F n v ). F c can be estimated based on the energy difference between the isotropic and nematic phase and the radius of the isotropic core (r c ). To evaluate F n v , we parameterize n(r) around the top vortex line along y in a möbiuson by using the ansatz Article https://doi.org/10.1038/s41567-022-01851-1 n(r) twists by 180° along a trajectory around the vortex line while smoothly connecting to the twist wall (Fig. 4d), consistent with the Möbius-strip-like structure. The free energy functional 1 contains the first three Frank-Oseen elastic energy terms with elastic constants K 11 , K 22 and K 33 describing the costs of splay, twist and bend deformations, respectively, as well as the dielectric contribution from the coupling to an external electric field |E| = U/d along the vertical (z) direction, where ε a is the LC's dielectric anisotropy and ε 0 is the vacuum permittivity. The ensuing free energy with respect to background per unit length is then derived via integrating equation (1) in the plane perpendicular to the vortex line ( Fig. 4d) thus where the radial integral from r c to p/4 is consistent with the observation that n(r) rotates by π/2 from the far-field alignment to its centre on either side of the vortex. The shrinking preference of the vortex line is manifested by F v > 0 (Fig. 4e). To calculate the energy for the twist wall, we assume an ansatz of n t = (0, sin( πx w ), cos( πx w )) for a wall with a uniform twist rate and width w (Fig. 4d). The associated linear energy density at an applied field reads with h being the height of the twist wall across the sample thickness (along z). Minimizing equation (3) gives the equilibrium width of an individual twist wall, w eq = √ stretching the vortex lines. The total free energy for the volume containing one fold of möbiuson within 0 < y < L is then calculated via integrating the linear energy densities along the respective paths of vortex lines and twist wall, F tot = ∫2F v dl v + ∫F t dl t , where The twist wall energy density in a möbiuson decreases with D t (Extended Data Fig. 9a), whereas that of the vortex lines increases with D v (Extended Data Fig. 9b). Minimization of equation (4) yields a structure with an equilibrium L = 11.2 μm and L/D t ≈ 1.75, consistent with experiments ( Fig. 4 and Extended Data Fig. 9c). Changes of twist wall displacements with varying |E| also qualitatively agree with experiments (Extended Data Fig. 9d). Increasing D t from 0 reduces the overall energy, revealing the folding instability of the linear fibre-like structure (D t = 0) and how möbiusons emerge from the competition between the voltage-dependent shrinking tendency of the vortex line and the expanding tendency of the twist wall (Extended Data Fig. 9).

Self-propulsion of möbiusons and control of topological cargo
When subjected to an oscillating field applied across the cell, along n 0 , möbiusons exhibit rotational and translational motions (Fig. 5a-c) dependent on their symmetry. In our LC with positive dielectric anisotropy, the oscillating electric field tends to rotate the möbiuson's n(r) towards n 0 , periodically squeezing it. The balances between the electric, elastic and viscous torques while periodically varying the applied voltage originate nonreciprocal dynamics of n(r) structures 2,6 and squirming motion of particle-like möbiusons dependent on their symmetry. N M and N R,L structures exhibit electrically powered translations and rotations, respectively. Meanwhile, 1 ∞ torons retain their axial symmetry in an oscillating electric field and display only Brownian motion, differently from LCs with negative dielectric anisotropy, where squirming motions emerge along with field-induced symmetry breaking of n(r) structures [6][7][8] . The linear and angular velocities of möbiusons with different symmetries (Fig. 5 and Supplementary Videos 7-9) depend on the frequency f of the applied field. For example, the translational speed of 2 M peaks at 0.4 μm s −1 for f ≈ 100 Hz (Fig. 5a,d and Supplementary Video 7). Above 100 Hz, the speed decays exponentially with f. Above 1,000 Hz, the translational motion of N M is indistinguishable from Brownian motion (Fig. 5a). At f ≲ 100 Hz, N M become unstable and shrink following the transformation rules ( Fig. 4b and Extended Data Fig. 1). The angular speeds of N R,L rotations depend on f, peaking at low frequencies and decreasing exponentially with f (Fig. 5b). The translational and angular speeds decrease as N increases. The direction Article https://doi.org/10.1038/s41567-022-01851-1 for 2 R and 2 L and for N R and N L (Fig. 5c), consistent with the fact that the left-right versions of möbiusons have similar 3D n(r) and differ only in terms of their 3D orientations relative to cell substrates (Supplementary Video 9).
Möbiusons can induce dynamics of high-symmetry topological solitons, such as torons, hopfions and skyrmions with axisymmetric n(r), which by themselves do not exhibit electrically powered motion in LCs with positive dielectric anisotropy. Möbiuson-soliton interactions (Fig. 6a) impart co-propulsions as shown by the solitons' trajectories for 2 M translating a toron (Fig. 6b), a skyrmion bag (Fig. 6c) and a hopfion (Fig. 6d) as topological cargos (Supplementary Video 10). Interesting dynamics emerge from the rotational motion of the 2 R interacting with a toron, which could be called 'topological dancing' (Fig. 6e, Supplementary Video 10). Compared with the self-propelling speed, 0.4 μm s −1 , of 2 M at 100 Hz, the speeds of 2 M -soliton complexes for torons, skyrmion bags and hopfions are 0.25, 0.07 and 0.15 μm s −1 , respectively, whereas the linear speed along the 'dancing' trajectory in Fig. 5e is 0.26 μm s −1 . The different co-propulsion speeds stem from the different effective viscous drag forces of these möbiuson-soliton assemblies, calling for their systematic analysis.

Conclusions and outlook
Möbiusons exhibit diverse folded configurations of co-assembled twist solitons and singular vortices in the order parameter field, both being topological mimics of the Möbius strip. Their non-orientable field configurations cannot emerge in polar systems, such as magnets or ferroelectrics, but potentially could have counterparts in other systems with nonpolar order parameters 1 , in polarization fields of light 26 and in other branches of science ranging from particle physics to cosmology [3][4][5]21 . Electrically powered emergent self-propulsions show möbiuson stability under out-of-equilibrium conditions and might help establish LCs as model systems to explore active matter behaviours and to transport cargo at the mesoscale 27 . While the confinement used in our work partially limits the diversity of folded configurations, even more possibilities could exist without confinement, including vortex-soliton assemblies with knotted configurations and Möbius-strip-like structures on larger scales. Such topology-based designs could enable novel mechanical responses in LC polymers and elastomers 28 and organize nanoparticles into mesoscale spatial patterns 29 . While chirality is known for stabilizing multi-dimensional solitons, we have shown that chirality can stabilize various co-assemblies of topological solitons and vortices, potentially enabling phases of matter analogous to twist grain boundary and blue phases in LCs and the A-phase in magnets [30][31][32][33][34][35][36][37] .

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41567-022-01851-1. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.

Materials and sample preparation
Chiral LCs were prepared by mixing 4-cyano-4′-pentylbiphenyl (5CB; EM Chemicals) or a commercially available nematic mixture E7 (Slichem) with a small amount of a left-handed chiral additive, cholesterol pelargonate (Sigma-Aldrich). To define the equilibrium pitch (p) of the ensuing chiral LCs, the weight fraction of the added chiral dopant was calculated using C dopant = 1/(h htp p), where the helical twisting power h htp = 6.25 m −1 for the cholesterol pelargonate additive. Additionally, for high-resolution, aberration-free 3D nonlinear optical imaging, a partially polymerizable chiral LC mixture was prepared by mixing E7-based chiral LC (69%), reactive mesogens RM-82 and RM-257 (at fractions of 12% and 18%, respectively; both from Merck) and ultraviolet-sensitive photoinitiator Irgacure 369 (1%; Sigma-Aldrich). The sample cells were assembled from indium tin oxide (ITO)-coated glass slides or coverslips treated with 25% solution of polyimide SE5661 (Nissan Chemicals) to obtain strong perpendicular (homeotropic) boundary conditions. The polyimide was applied to the substrates by spin-coating at 2,700 rpm for 30 s followed by baking (for 5 min at 90 °C and then for 1 h at 180 °C). The LC cell gap thickness was defined to be 10-15 μm by using silica spheres as spacers between two substrates, and the ratio between the cell gap and the pitch was d/p = 2.2. Metal wires were attached to the ITO and connected to an external voltage supply (GFG-8216A; GW Instek) for electric control.

Generating and imaging localized structures in chiral LCs
To generate a möbiuson, we utilized an ytterbium-doped fibre laser (YLR-10-1064, IPG Photonics), operating at 1,064 nm. First, we laser-generated an extended fibre-like fragment in a uniform unwound background under a sinusoidal electric field by using a laser power between 30 and 40 mW (Supplementary Video 1). Then, we switched off the laser and carefully adjusted the voltage to be around 3.0 V such that the fibre-like fragment did not shrink or extend. The two ends of the fibre-like fragment then spontaneously curved, forming möbiusons of various orders and symmetries, depending on the length of the fibre-like fragment and the spontaneous initial curving directions of its ends. High-N structures (N > 20) could be generated from a long fibre-like fragment or by adding extra orders to a low-N möbiuson (Supplementary Video 3). Torons, hopfions and skyrmion bags were also generated by laser tweezers by following procedures as reported previously 24,34,35 . Bright-field and polarizing optical microscopy images were obtained by using a multi-modal imaging setup built around an IX-81 Olympus inverted microscope (which is also integrated with the 3PEF-PM imaging setup described below) and charge coupled device cameras (Grasshopper and Flea FMVU-13S2C-CS, Point Grey Research or Spot 14.2 Color Mosaic, Diagnostic Instruments) with 100×, 40× and 20× Olympus objectives with high numerical aperture (NA) of 1.4, 0.75 and 0.4, respectively. The translational speed, angular velocity and trajectories of solitons were analysed by using ImageJ software (National Institutes of Health).

Three-dimensional nonlinear optical imaging
The 3D nonlinear optical imaging of möbiuson structures such as 2 M and 2 R was performed using the 3PEF-PM setup built around the IX-81 Olympus inverted optical microscope integrated with the ytterbium-doped fibre laser 25 . To reduce material birefringence that leads to 3PEF-PM imaging artefacts 35 , before imaging, the möbiusons stabilized in the partially polymerizable chiral LC mixture were photo-polymerized by weak ultraviolet illumination and the unpolymerized LC was washed away and replaced by infiltrating an indexmatching immersion oil. We used a Ti-sapphire oscillator (Chameleon Ultra II; Coherent) operating at 870 nm with 140 fs pulses at a repetition rate of 80 MHz as the source of the laser excitation light. An oil-immersion 100× objective (NA 1.4) was used to collect the fluorescence signal. The fluorescence signal was detected by a photomultiplier tube (H5784-20; Hamamatsu) after a 417/60 nm bandpass filter. The 3PEF-PM imaging involves a third-order nonlinear process 25 , during which LC molecules are excited via the three-photon absorption process and the signal intensity scales as cos 6 β, where β is the angle between the polarization of the excitation light and the long axis (transition dipole) of the LC molecule. Different polarization states of the excitation were controlled by a pair of half-wave and quarter-wave retardation plates. The patterns of strong 3PEF-PM intensity obtained for a linearly polarized excitation beam reveal the spatial regions where cosβ is large, reflecting where n(r) is close to parallel to the polarization of the laser beam (Extended Data Figs. 4 and 5). On the other hand, the laser excitation with circular polarization helps to visualize regions with n(r) parallel to glass substrates regardless of their azimuthal orientations (Fig. 3g,h). The experimental 3D fluorescence images were processed by applying uniform contrast enhancement and depth-dependent intensity correction and then compared with their numerically simulated counterparts.

Numerical modelling and visualization
We reconstructed the director field configuration of solitons by minimizing the Landau-de Gennes free energy of a continuum representation of a chiral LC, which reads F = ∫ ∫f surface d 2 r + ∫ ∫ ∫f bulk d 3 r, (5) where the bulk energy density ( f bulk ) was integrated over the 3D volume occupied by the LC whereas the surface free energy density ( f surface ) was integrated over the 2D surfaces at the top and bottom interfaces with substrates. The bulk energy density in the Q tensor formalism can be written as 10,11,38 where ε ikl is the Levi-Civita symbol, L i are the elasticity parameters, A, B and C are the nematic material parameters, ε is the average LC permittivity and ε mol a is the molecular dielectric anisotropy. Summation over all indices is assumed. The molecular tensor order parameter is defined by Q = S 2 (3nn − I ) where S is the scalar order parameter and I is the identity matrix. The first three terms in equation (6) correspond to the thermotropic energy that describes the nematic-isotropic transition of the LC, and minimization of the thermotropic energy yields the equilibrium scalar order parameter S eq = −B+√B 2 −24AC

6C
. The next five terms correspond to the elastic energy that comes from the distortion in n(r), including the terms related to molecular chirality. Together, the thermotropic and elastic terms reflect the energy penalties associated with deviating the molecular director away from its equilibrium configuration and order. The final two terms describe the dielectric coupling between Q and the electric field E. The surface free energy density f surface describing the surface anchoring at the substrates reads where W is the anchoring strength and Q (0) ij defines the preferred orientation and order (S eq ) of the LC at the surfaces, corresponding to perpendicular boundary conditions in our experiments.
An energy minimization routine based on the variational method was performed to minimize equation (5) 10,11,21 . To reconstruct molecular alignment configurations of möbiusons, n(r) orientations derived https://doi.org/10.1038/s41567-022-01851-1 from 3PEF-PM imaging were used as initial conditions. For all the simulations, the following parameters were used: 10-12 d/p = 2, A = −1.72 × 10 5 J m −3 , B = −2.12 × 10 6 J m −3 , C = 1.73 × 10 6 J m −3 , L 1 = 3.3 × 10 −11 N , L 2 = 5.3 × 10 −11 N , L 3 = 0 N , L 4 = 9.4 × 10 −11 N , L 6 = 3.5 × 10 −11 N, W = 5 × 10 −3 J m −2 , which yields S eq = 0.533 and the nematic correlation length ξ = 6.63 nm. The voltage applied across the cell gap was U = 7.0 V. Minimization of möbiusons' free energy was performed on a 3D grid with a spacing of 33.9 nm in all directions, and a typical simulation volume was a cuboid box with 64 × 112 × 50 grid points, corresponding to physical dimensions of 2.17 × 3.80 × 1.70 μm 3 . Doubling the grid spacing (to 67.8 nm) on a 128 × 224 × 100 grid (8.6 × 15.2 × 6.8 μm 3 ) did not yield qualitative differences in the minimized field configurations. We note that using grid spacings larger than ξ leads to an underestimation of the free energy of the defect cores. Nevertheless, variations in field configurations at micrometre length scales, which are three orders of magnitude larger than ξ and the size of the defect core, rather than the system's behaviour within the cores of defect lines, are most relevant to our study of möbiuson structure in this work. The excellent agreement between the simulated field configurations and their experimental counterparts justifies the used rescaling and grid spacing, which allowed for modelling of structures containing defects at the expense of nanometre-scale details within their cores. Minimizations of extended fragment structures were performed on a grid with 100 × 100 × 100 points and spacing of 10.17 nm. To gain insights into the origins of the folding instability, we simulated the snake-like fragment of the möbiuson using a computational box with 50 × L × 50 grid points, where L is the length of a repeating unit cell containing the periodically folded structure to be determined by energy minimization. The boundary conditions at the edges of this unit cell on planes orthogonal to the y axis in Fig. 2l are periodic. By minimizing the field configurations with different L values, we obtained the energy-minimizing period L, as well as the equilibrium director configuration and other geometric characteristics of the folded structure. The numerical results agree well with experimental observations. We simulated POMs by the Jones matrix method using the energy-minimizing n(r) of solitons 21,38,39 . We first spilt the cell into 50 thin sublayers along the z direction, then the Jones matrix for each pixel in each sublayer was calculated by finding the local optical axis and the ordinary and extraordinary phase retardation. The optical axis is along the local molecular direction, and phase retardation comes from the optical anisotropy of the LC, for example, n o = 1.58 and n e = 1.77 for 5CB. The Jones matrix for the whole LC cell was obtained by multiplying all the Jones matrices for each sublayer. The simulated singlewavelength POM was obtained as the analyser component of the product of the Jones matrix and the incident polarization. Coloured POMs were then obtained as the combination of single-wavelength POMs corresponding to red, green and blue primary colours at wavelengths of 650, 550 and 450 nm, respectively, with intensities matching the light source in the experiment.
The 3D field configurations of solitons were visualized by the isosurfaces of n z , the vertical component of n(r), at |n z | = 0.8, coloured according to the orientations of n(r), as well as by the isosurface of regions of reduced scalar order parameter (S = 0.45 for fibre-like fragments and S = 0.5 for möbiusons) to visualize singular vortex lines. Computer simulations of the 3PEF-PM images were based on the signal intensity I 3PEF−PM ∝ cos 6 β. In both our experiments and simulations, the intensity patterns obtained for different linear polarizations of excitation light were coloured as shown in Extended Data Figs. 4 and 5. Visualizations were performed in ParaView open-source freeware.

Data availability
Source data are provided with this paper. All other data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes used for the Landau-de Gennes modelling and the simulation of POMs are published in ref. 38 . Other codes used for the numerical calculations are available upon request.
Article https://doi.org/10.1038/s41567-022-01851-1 Extended Data Fig. 2 | Twist wall and near-surface twist disclinations of a fibre-like fragment. a, Schematic of a twist wall. The non-polar director field with head-tail symmetry is shown by cylinders that twist by 180°, matching the uniform unwound far-field background n 0 along the z direction. b, Free energy density of an extended fibre-like fragment stabilized at an applied electric field in the cross-sectional plane perpendicular to the direction of translational invariance. The average free-energy density of the twisted region (region I) is lower than that of the unwound far-field background (region II) and is plotted relative to the far-field background. The increased values of free energy density in the middle of the twist wall is higher due to the contribution of the dielectric coupling term. c, Free energy density along x coordinate in the z-midplane marked in (b). The blue dashed line corresponds to the free energy density of the far-field background.