Sculpting liquid crystal skyrmions with external flows

We investigate, using experiments and numerical simulations, the distortions and the alignment of skyrmions in the liquid crystal under external flows for a range of average flow velocities. The simulations are based on the Landau-de Gennes $Q$ tensor theory both for isolated as well as for systems with many skyrmions. We found striking flow-driven elongation of an isolated skyrmion and flow alignment of skyrmions in the many-skyrmion system, both of which are also observed in the experiments. In the simulations, particular attention was given to the dissipation rate and to the various dissipation channels for a single skyrmion under external flow. This analysis provides insight on the observed scaling regime of the elongation of isolated flowing skyrmions and revealed a surprising plastic response at very short times, which may be relevant in applications based on the alignment of soft structures such as liquid crystal skyrmions.

We investigate, using experiments and numerical simulations, the distortions and the alignment of skyrmions in liquid crystal under external flows for a range of average flow velocities. The simulations are based on the Landau-de Gennes Q tensor theory both for isolated as well as for systems with many skyrmions. We found striking flow driven elongation of an isolated skyrmion and flow alignment of skyrmions in the many-skyrmion system, both of which are also observed in the experiments. In the simulations, particular attention was given to the dissipation rate and to the various dissipation channels for a single skyrmion under external flow. This analysis provides insight on the observed scaling regime of the elongation of isolated flowing skyrmions and revealed a surprising plastic response at very short times, which may be relevant in applications based on the alignment of soft structures such as liquid crystal skyrmions.
Liquid crystals (LCs) are unique soft materials that exhibit facile responses to external electric and magnetic fields. Additionally, the LC molecular orientation field, the director, can be easily distorted by anchoring at solid boundaries. These properties of LCs are the cornerstone of modern display technologies [1]. In LC displays the director field undergoes continuous transitions between distinct topologically trivial configurations, i.e., configurations that can be morphed into a uniform structure continuously. Topological point or line defects may be spontaneously nucleated but these singular configurations are usually short-lived transient structures.
Recently, a new type of topologically non-trivial nonsingular director configurations was realised experimentally [2]. These solitonic structures are a soft matter analogue of magnetic skyrmions, which are two-dimensional solitons with localized winding of the magnetic moment observed in hard condensed matter systems [3,4], which are relevant for technological applications such as spintronics [5]. Likewise, LC skyrmions are localized topologically protected distortions of the director field, which in some realisations can even include knotted or linked field lines [6,7]. Additionally, skyrmions exhibit particlelike behaviour, such as long-ranged effective interactions between skyrmions and between skyrmions and colloids [8,9], self-assembly into spatially periodic structures [10], and may acquire self-propelled motion when subject to periodic time-dependent electric fields [11,12].
As topological structures, skyrmions may be indexed or classified in terms of the elements of homotopy groups [6]. Experiments and numerical simulations revealed electric and magnetic transitions between skyrmionic textures with distinct topological indices [13]. It has been suggested that this behaviour may be used to control three-dimensional orientational structures of anisotropic nanoparticles, aligned with the local LC director, and drive them through transitions between different topological configurations [14,15].
Despite the extensive body of research on LC skyrmions, their interaction with externally imposed material flow fields remains poorly understood. In our recent paper [16] we reported a numerical study, based on the Ericksen-Leslie nematohydrodynamics, of the effect of external flows on the structures and dynamics of LC skyrmions. We found a configurational transition driven by the flow, from skyrmionic distortions along the flow in weak flows to skyrmionic distortions in the direction perpendicular to the flow in strong flows.
Here, we extend that study and address explicitly the reversibility of the flow induced skyrmionic distortions. This analysis provides insight on the distortions observed in skyrmions under flow and reveals the mechanisms of their plastic response, which may be key in applications based on the elongation and alignment of soft structures arXiv:2303.12045v1 [cond-mat.soft] 21 Mar 2023 such as liquid crystal skyrmions.
We reformulate the problem in terms of the Landaude Gennes Q tensor theory and analyse the different contributions to the dissipated energy. Our numerical results for the flow-induced elongation of a skyrmion agree qualitatively with experimental observations. We also report on the collective behaviour of skyrmions driven by external flows observed both in the simulations and in experiments. The most significant many-body effects are the hindering of skyrmionic elongations by the presence of other skyrmions and the collective skyrmion alignment by the external flow. see Fig. 1 fl o w a b FIG. 1. Schematic illustration of the simulation setup. a) Initial configuration of the director field. b) Elongated skyrmion as a result of the externally applied mass flow. The arrows represent the director field while the colours stand for its z-component (black for nz = −1 and yellow for nz = 1).

Flow-induced elongation of isolated skyrmions: numerics
We start by describing the elongation of a single skyrmion due to the externally imposed material flow. In the simulations, the flow is driven by an external constant force density field. The flow field rapidly reaches a uniform steady state, due to the friction at the surfaces, which enables an accurate description of the skyrmion behaviour under controlled conditions.
Black corresponds to |nz| = 1 and yellow to |nz| = 0. The flow field is directed from the top to the bottom of the panels as is indicated by the white arrows. The scale bars correspond to one cholesteric pitch p.
average flow velocity u = 4636.15 µm/s (t ch ≈ 0.002s). Starting from a circular shape, Fig. 2(a), the skyrmion becomes axisymmetric with a small tail-like wake, Fig. 2(c), and undergoes continuous elongation in the flow direction (from the top to the bottom of the panel) beyond a characteristic time t ch = p/ u , where p is the cholesteric pitch. t ch is the time taken by the skyrmion to move by one cholesteric pitch, which is roughly the skyrmion size at t = 0. The skyrmion elongation regime is illustrated in Fig. 2(d)-(f), where the configurations are calculated at t = 100t ch , 300t ch , 500t ch , respectively. In Fig. 3(a), we plot the time evolution of the skyrmion aspect ratio or elongation ε, defined as the ratio of its length (parallel to the flow) and its width (perpendicular to the flow). We define the boundary of the skyrmion as the set of points where |n z | = 0.5. The four curves in Fig. 3(a), corresponding to different flow velocities u , collapse (approximately) onto a master curve when plotted against t/t ch . The data collapse is better at short times t 10t ch , which indicates that t ch is the relevant time scale in this regime. The circular skyrmion shape starts to distort at t 2t ch (see Fig 2(b), and configuration 2 in supplemental Figs. S1-S4). Significant elongations, with ε > 1.5, are observed at t 20t ch for the lowest u (orange diamonds in Fig. 3(a)), and at t 40t ch for the other velocities (blue circles, green triangles, and red squares in Fig. 3 We observe that at t 30t ch when ε 2.0 (compare the orange diamonds with the other symbols in Fig. 3(a)) the scaling of the elongation ε starts to break down. In particular, the data for sample 4 (orange diamonds) shifts systematically upwards relative to the master curve, i.e. the skyrmion elongation proceeds at a higher rate when t 30t ch . Sample 3 follows the master curve until t ≈ 70t ch when ε ≈ 2.3, and then shifts upwards similarly to the previous case. Finally, at t 150t ch , when ε ≈ 3.0, the data for sample 2 (green triangles) departs from the data for sample 1. This cascade of thresholds where the scaling behaviour of the elongation breaks down can be related to skyrmion shape changes as illustrated in supplemental Fig. S5(a)-(c), for samples 4-to-2 at times t = 30, 70 and 150t ch , respectively. All the configurations have a thin tail-like region on the left ( Fig. S5(a)-(c)), which after forming lags behind the head-like front on the right, leading to effectively faster skyrmion elongations. For comparison, we also show in supplemental Fig. S5(d) the configuration for sample 1 at t = 150t c h which lacks a well defined tail-like feature. As a result this configuration elongates more slowly.
With the aim of comparing the simulation results with experiments we depict in the inset of Fig. 3(a) the aspect ratio ε as a function of the skyrmion displacement r. The latter is defined as the distance between the skyrmion's centre of mass at a given time and its initial position. We recall, that the simulations stop for all samples at the same absolute time, corresponding to different reduced times and maximal skyrmion displacements. The curves exhibit a quasi-linear behaviour and good data collapse is observed for r 10p, in line with the scaling of ε(t/t ch ) discussed above.
The linear behaviour is clear at low velocities u , or large characteristic times (same pitch), see the orange diamonds in the inset of Fig. 3(a). At higher velocities, we find ε ≈ 1 + 0.009r/p for sample 2 (green triangles). The breakdown of scaling of the elongation at r 10p (t 10t ch ), and the observation that in this regime ε grows faster at lower u (compare the orange diamonds with the red squares in Fig. 3(a)) are related to differences in the director field configurations at a given t/t ch for different flow velocities u . Figures. 3(b)-(e) depict the skyrmion configurations under different flow velocities u at t ≈ 10t ch , while Figs. 3(f)-(i) depict the configurations at a later time, t ≈ 50t ch . The first four configurations exhibit similar shapes and are in the scaling regime revealed in Fig. 3(a). Significant shape changes, however, are observed at t ≈ 50t ch as illustrated by the last four configurations in Fig. 3(f)-(i). Thus, the slowest skyrmion ( Fig. 3(f)) is already in the elongation regime, with wellpronounced head-and tail-like regions. By contrast, the fastest skyrmion ( Fig. 3(i)) has a bullet-like shape with a flat head (recall that the flow direction is from top to bottom.) Furthermore, the skyrmion shape becomes fuzzier we discuss only the dissipation rate D(t) arising from the liquid crystal flow and the director reorientation in the bulk. As we will show below, D(t) sensitively probes the dynamics of the skyrmion distortions, and provides a better understanding of the failure of the scaling of the elongation ε(t) at t ≥ 10t ch . First, we note that the scalar order parameter S hardly changes in the course of the simulations, and thus we use the Ericksen-Leslie theory to calculate the dissipated energy. We follow Ref. [17] to relate the parameters of the Beris-Edwards to those of the Ericksen-Leslie equations. In the Ericksen-Leslie theory the dissipated energy by a flowing liquid crystal is given by [18]: where the α n 's are the Leslie viscosities, γ 1 = α 3 − α 2 , γ 2 = α 6 −α 5 and the co-rotational time flux of the director is defined as N α = ∂ t n α + u β ∂ β n α − W αβ n β . Diagonalisation of the Q tensor field yields the director field up to a sign. Therefore, we carried out a standard procedure of director vectorisation. The resulting vector field is continuous and can be used to calculate the dissipated energy using Eq. 1. We verified that upon applying the constant external force density, the fluid velocity away from the skyrmion relaxes almost instantaneously, as compared to the relaxation of the director field. Next, due to the 2D nature of this problem, this constant flow field away from the skyrmion cannot perturb the director field which has n z = 1 in that region of the simulation domain. Therefore, the time dependence of all the contributions to the total dissipation rate discussed below is due solely to the director dynamics in the skyrmion proximal region. We found that the most important contribution to D(t) is due to the rotational viscosity term γ 1 N α N α in Eq. 1. Supplemental Fig. S7, compares D(t) with the dissipation rate D 1−5 (t) of the other five terms ∝ α 1 , γ 2 , α 4 , (α 5 + α 6 ) in Eq. 1. The conclusion is that the contribution D 1−5 is an order of magnitude smaller than the dissipation arising from the γ 1 term. In the vicinity of the skyrmions the local flow differs from u (see Fig. S6(e)-(h)), which results in a weak temporal variation of the frictional dissipation rate as shown in supplemental Fig. S8. The quantity plotted is the reduced frictional dissipation rate D f (t) = S d 2 rχρ(u − u ) 2 , defined as the excess over the frictional dissipation for a uniform flow u . D f (t) is two orders of magnitude smaller than D(t), and thus can be safely ignored. Figure. 4(a) illustrates D(t) as a function of the reduced time t/t ch for several values of u . We find that the three curves corresponding to the lower values of u collapse approximately onto a single curve, while the curve corresponding to the largest u stands apart. This can be understood qualitatively by inspecting the average velocity fluctuations ∆u(r) = (u(r) − u ), which is at least three times larger than for the other three samples (compare the colour bar in supplemental Fig. S6(h) with those on Fig. S6(e)-(g)). Large velocity fluctuations |∆u| enhance significantly the local velocity gradients, which in turn couple to the director field and drive it locally away from the preferred orientation. This is one of the reasons why the dissipation rate for sample 1 (red squares in Fig. 4(a)) lies markedly above the other three curves.
The regular symmetric flow pattern shows that the system is not in the turbulent regime.
Despite this difference, the curves exhibit similar qualitative behaviour: , with t * /t ch ≈ 100, 20, 10, 2 for samples 1 to 4, respectively. Interestingly, these values of the threshold time t * are an order of magnitude smaller than the characteristic time scale γ 1 p 2 /L ≈ 1.3s (where p is the cholesteric pitch and L is the average elastic constant) for the relaxation of the director field.
We emphasize that this separation into three regimes is only qualitative. For instance, we note a more complex behaviour of D(t) within regime 2) for the three highest flow velocities, where we observed short plateau-like features for 3 t/t ch 5, see supplemental Figs. S1-S3. The emergence of such features highlights the complex skyrmion distortion dynamics even when the skyrmion shape remains almost constant at short times, see the skyrmion configurations labelled 2 and 3 in supplemental Figs. S1-S3, for example. Interestingly, this stagnation of the skyrmion elongation is accompanied by a burst in the dissipation rate due to the coupling of the flow gradients to the director field as shown in the bottom panel of supplemental Fig. S7. This panel reveals that the dynamical regime 1) t/t ch 1 is also observed in the relaxation of the dissipation term D 1−5 .
In regime 1) the dissipation rate D(t) is reduced by a factor of two, while the skyrmion maintains its circular shape with ε ≈ 1, as shown by the configuration depicted in Fig. 2(b) obtained at the end of this regime, see also the configurations labelled 1 in supplemental Figs. S1-S4. Next, in regime 2) the total dissipation rate D(t) is further reduced by approximately an order of magnitude, while the elongation is ε 1.5 for samples 2, 3 and 4 as shown in Fig. 2 At late times, in regime 3) the skyrmions start to elon-gate and the total dissipation rate D(t) takes an almost constant value (ca. 20 lower than the dissipation rate at early times). This behaviour is illustrated in supplemental Fig. S1 (panels 7, 8), Fig. S2 (panels 6, 7), Fig. S3 (panels 6, 7), and Fig. S4 (panels 5, 6). The elongation dynamics is more pronounced at the highest flow velocity (sample 1), panels 7 and 8 in supplemental Fig. S1. The skyrmion configuration 8 extends almost throughout the whole system, reaching the limiting value of the aspect ratio. This fact is reflected by the saturation of the Landau-de Gennes free energy F LdG (t), as shown in Fig. 4(b) (red squares for t 1000t ch ). At shorter times, F LdG (t/t ch ) exhibits a linear behaviour, and all the curves collapse when the time is scaled by the characteristic time t ch . A closer inspection of the curves, at the three highest flow velocities, reveals the existence of plateau-like behaviour at 10 t/t ch 20, where the skyrmion aspect ratio hardly changes, as illustrated in supplemental Figs. S1-S3 (panels 4 and 5). Again, this observation emphasizes a complex non-monotonic behaviour of the skyrmion elongation dynamics, which, among others, leads to the pearly string shapes shown in Fig 2(f).
We now proceed to discuss the reversibility of the skyrmion distortions induced by the flow. To this end, we calculate the Landau-de Gennes free energy as a function of time along two consecutive cycles when the flow field is switched on and off. Figure 5 illustrates the results for several values of the on-flow cycle. When the flow is on, after a short transient, the Landau-de Gennes free energy increases linearly with time (see open red squares in Fig. 5)-the work of the external force density is stored in the elastic distortions of the skyrmion. When the flow is switched off, a fraction of the stored elastic free energy is dissipated, and F LdG drops to a constant value, which depends on the duration of the preceding on flow cycle (see the open triangles, circles and diamonds in Fig. 5). Somewhat unexpectedly, the free energy does not relax to its pre-flow value even when the flow was applied for periods as short as 7t ch , as shown by the green triangles in Fig. 5. Correspondingly, the flow-elongated skyrmions (Figs. 5(b), (d), and (f)) do not return to their original circular shape (Figs. 5(a)) and remain in distorted metastable configurations (Figs. 5(c), (e), and (g)). This reveals that flow-induced skyrmion distortions exhibit surprising plastic behavior even at elongations as small as 10%. A natural question arises whether one can design spatio-temporal flow patterns which would result in some pre-defined skyrmion shape?  Figure 6(h) depicts the aspect ratio as a function of time for two toron realizations and flow intensities, while the inset shows the average toron velocity as a function of time. As the toron velocity is not constant, the results plotted in Fig. 6(h) can not be compared directly with those in Fig. 3(a). However, by plotting the aspect ratio as a function of the toron displacement, Fig. 6(i), we observe data collapse similar to that found in the simulations (inset of Fig. 3(a)). This observation suggests that, if the toron velocities were constant, there would be a characteristic time for which the curves for the elongation ε would collapse. It is noteworthy that even with the varying toron velocity, ε is approximately linear in the displacement. For the data shown in Fig. 6(i) we find ε ≈ 1+0.17r/p, where the slope is about 18 times larger than that found in simulations.
We speculate that this difference may be due to the 2D approximation used in the simulations. Nevertheless, it is encouraging that even a 2D model agrees qualitatively with the experimental observations.

Collective effects of flowing skyrmions/torons
When many skyrmions are subject to a flow field, their elongation may be hindered or even reversed as a result of effective interactions mediated by distortions of the LC director. In Ref. [16], two flowing skyrmions were simulated, and the resulting elongation was found to be smaller than that of an isolated skyrmion. The skyrmionskyrmion interaction becomes relevant at distances of the order of the cholesteric pitch.
We skyrmions fill the entire domain. We also observed that some skyrmions shrink as they are caged by neighbouring skyrmions, as shown by the arrows in Fig. 7(a)-(d).
We calculated the average skyrmion aspect ratio and compared it with the single skyrmion case in Fig. 7(e). Initially, the two curves evolve similarly, but at t ≈ 0.2s they start to deviate indicating that a skyrmion in the ensemble elongates less than an isolated skyrmion. This occurs as the interactions with neighbouring skyrmions hinder the elongation of an isolated skyrmion. The error bars that increase with time indicate that the skyrmion distortions become more heterogeneous, as illustrated in Fig. 7(a)-(d). Similar behaviour was observed in experiments with ensembles of flowing torons as shown in Fig. 7(f)-(i). The toron indicated by the arrow initially elongates and then shrinks back to an almost circular shape as the result of interactions with the elongating neighbours.
Clearly, the elongation of flowing skyrmions discussed in the previous section occurs only if they are isolated or if their elongation proceeds along a common direction. The time dependence of the free energy of the system with 30 skyrmions is similar to that of the average aspect ratio (see supplemental Fig. S9). At short times it is similar to that of an isolated skyrmion, flowing at the same average velocity, and it exhibits many-body effects beyond t ≈ 0.2s as the interactions between the elongated skyrmions become relevant. At late times the free energy increases more slowly in line with the suppression of the average elongation of the skyrmions.

DISCUSSION
We investigated, using experiments and numerical simulations, the distortions and the alignment of liquid crystal skyrmions under external flows for a range of average flow velocities.
We performed extensive simulations based on the Landau-de Gennes Q tensor theory both for isolated as well as for systems with many skyrmions. The most striking effects were the flow driven elongation of single skyrmions and their alignment along the flow direction in the many skyrmion system, both of which were also observed in the experiments.
The simulations were run for different flow velocities u and the results revealed a characteristic time scale t ch = p/ u , which controls the elongation of single skyrmions, leading to data collapse or scaling of ε as a function of t/t ch .
Specific features were observed depending on u : the larger the flow velocity u , the larger the extention of the scaling regime. This was related to a configurational transformation where the skyrmion takes a comet-like shape where the tail moves more slowly than the head. In other words the skyrmion is pulled apart by the flow, with distinct regions responding to the flow in different manners.
Furthermore, due to finite size effects, observed at high flow velocities and long times, the free energy flattens or saturates and remains almost static as the flow induced elongation is hindered by the image skyrmion that results from the periodic boundary conditions.
We have also calculated the dissipation rate and the various dissipation channels for a single skyrmion under external flows. This analysis provided insight on the observed scaling regime of the elongation of a single skyrmion and revealed a somewhat surprising plastic re-sponse at very short times. The plastic character of the skyrmion elongation occurs even at small strains, where one would expect an elastic response to the forcing. This plastic behaviour is expected to be present also in the many skyrmion system and may be relevant in applications based on the alignment of soft structures such as liquid crystal skyrmions.
The experimental observations of the flow induced distortions and alignment of single and many-skyrmion systems agree qualitatively with the numerical results, despite important differences between the model and the real system. We single out the fact that the simulations were carried out in 2D under controlled uniform flow conditions, while the experiments are in 3D and the flow was not uniform. Nevertheless the elongation and the alignment of the distorted skyrmions by external flows appear to be robust and are qualitatively similar. The plastic behaviour uncovered by our analysis of the dissipation may also have an impact on experimental observations and applications.
We finish by reporting an interesting experimental observation that is not reproduced by the 2D model. Under certain experimental conditions, skyrmions have been observed to de-stabilize and transform into cholesteric fingers, which elongate in random directions. Preliminary work suggests that applying a flow field for a short time results in orientational ordering of the elongated cholesteric fingers, which persists when the flow is switched off.
Previous experimental work reported that torons were stabilized by applying a weak voltage across the cell and were observed to transform into cholesteric fingers [19] when reducing this voltage. The fingers grow in random directions triggered by a decrease of the voltage below a certain threshold. Each toron (initially with a spherical shape in Fig. 8(a)) elongates in a random direction ( Fig. 8(b)) and continues to grow until filling the available space. Figures 8(c) shows growing fingers, but now under an applied external flow that is switched on for a short time just after the voltage reduction. We found that it was sufficient to maintain the flow for 3s in order to align all the fingers in the flow direction, as shown in Fig. 8(c). Moreover, after switching the flow off, the fingers continue to elongate along the same direction, filling the entire space. This suggests that an external flow may be effective to drive orientational ordering of cholesteric fingers. Interestingly, this ordered texture of fingers remains stable for days, although, some of the fingers shrink back to their original circular shape under an effective compressing action of the neighbouring fingers.
There exist several distinct types of cholesteric fingers [9]. To determine the type of the fingers observed here, we carried out a quasi-static numerical analysis of the 3D structure of the elongated torons, by minimizing the Frank-Oseen elastic free energy. Several, examples of the stretched toron configuration are presented in Figs. 8(d)-(f), using colour coded director preimages, i.e., those spatial regions where the nematic director n has a given fixed value. The vectorial representations of n(r), shown in Fig. 8(g) and (h), reveals that the topology of this extended structure remains unchanged, because the vectorized n(r) taken on the (x, z) cross section of this extended configuration covers the order parameter space S 2 once. This finding is also confirmed by directly verifying that the skyrmion number 1 4π n · ∂n ∂x × ∂n ∂z dxdz = 1. The corresponding skyrmion number surface density N sk in Fig. 8(i) highlights the existence of two nonsingular λ−disclinations, or fractional skyrmions [20], at the top and the bottom of the cross section. In fact, this configuration is characteristic of a cholesteric finger of the second type, CF-2, studied previously [21], and is often found as a metastable configuration occurring spontaneously. The director configuration in CF-2 shows how the fractional values of the skyrmion number, two half skyrmions in this case, can add up to unity in the finger, being embedded into the uniform far field background.
Figures 8(j)-(l) compares computer-simulated and experimental polarized optical microscopy (POM) as well as three-photon excitation fluorescence polarizing microscopy ((3PEF-PM) images for a system of aligned extended torons. We observe very good qualitative agreement, which demonstrates that the experimental elongated structures are cholesteric fingers of the second type. It is natural to assume that similar director configurations are associated with the torons stretched by the flow, being subject to stabilising voltage at the same time. More detailed numerical investigation of these 3D structures and their interaction with the external flow is beyond the scope of this study, which will be addressed in a future work.
Other extensions of this study include the effect of a time-dependent, e.g. step-like, flow field, the introduction of local changes in the flow direction or the effect of flows with non-zero vorticity.

Beris-Edwards model for flowing skyrmions
The model used to describe the liquid crystal skyrmions under mass flow is based on the Landau-de Gennes Q tensor theory of nematics [22,23]. For uniaxial ordered phases, the tensor order parameter is Q αβ = S(n α n β − δ αβ /3), where S is the scalar order parameter, which equals to zero in the isotropic phase and unity in phases with perfect alignment, and n α is the director field. The free energy of the liquid crystal is F = V d 3 r f , where the free energy density f has three contributions: The first is the bulk free energy and a, b and c are positive material constants. The second is the elastic free energy with two elastic constants L 1 and L 2 [24,25] and p is the equilibrium cholesteric pitch. The last contribution arises from the anchoring at the confining surfaces. As the simulations are effectively 2D, one needs to apply this anchoring everywhere (as a bulk term) to stabilize the skyrmions [16,26,27]. The parameter W controls the anchoring strength. The dynamics is governed by the Beris-Edwards, the continuity and the Navier-Stokes [28] equations: where ρ is the liquid crystal density and Γ is the system dependent rotational diffusivity. Equation 5 describes the time evolution of the order parameter Q αβ , which depends on the velocity field u α . The dynamics of the velocity field is given by Eq. 7, which depends on Q αβ . The first term on the right-hand side of Eq. 7 stands for the effective friction at the surfaces [16]. For Poiseuille flow of a fluid with absolute viscosity η = ρν between two surfaces separated by a distance L, the friction coefficient in 2D that yields the same average velocity of the 3D flow is χ = 12η/L 2 , where ν is the kinematic viscousity. The fluid moves driven by an external force f α . In the equations above, the shear rate and the vorticity are given by A αβ = (∂ α u β + ∂ β u α )/2 and W αβ = (∂ β u α − ∂ α u β )/2 and the co-rotational term reads: where ξ is the flow aligning parameter. The molecular field is: and the nematic stress tensor is given by: with P 0 the hydrostatic pressure.

Numerics
As discussed in Ref. [16], the time scales for changes in the director and the velocity fields differ by six orders of magnitude, rendering the simulations challenging. The simulations reported in Ref. [16] use finite differences for both fields (n α and u α ) and an adaptive time step to speed up the convergence of the velocity field. The dynamics of the director field was described by the Ericksen-Leslie theory. It was also assumed that the fluid relaxes instantaneously when compared to the director field. Here, we use the Beris-Edwards equation for the dynamics of the director field, which is solved using finite differences (predictor-corrector algorithm) while the lattice Boltzmann method is used for the velocity field. This approach results in slower but more reliable simulations, since the dynamics of the two fields are solved with the same time step. In addition, the lattice Boltzmann method is conservative while the method based on finite differences is not. Perhaps the major advantage of this approach is that the simulations are stable for parameters close to the experimental ones, which was not the case for the Ericksen-Leslie theory [16].
Due to the different time scales, which require long simulations, a compromise between the domain size and the simulation time is required. We used 2D simulation boxes, which provide useful insights into the dynamics of skyrmions in feasible simulation times. However, a 2D domain assumes invariance in the direction perpendicular to the plane, which may be a strong approximation for the velocity field. In addition, the skyrmion structure in the experiments is not invariant in the direction perpendicular to the plane. Thus, quantitative differences between the simulations and the experiments are inevitable although the qualitative behaviour is expected to be captured. For instance, the flow velocities in the simulations of stable skyrmions can be much higher than those in the experiments, as shear flows in 2D are rather weak by comparison to 3D where these strong flows destroy the skyrmions.
The parameters used in the simulations are provided in Table S1 (in the supplemental material) both in simulation and physical units. We assume the material parameters of 5CB at room temperature [11,29,30]. Some of these were converted from the parameters of the Ericksen-Leslie theory to the Beris-Edwards one using the expressions given in Ref. [17] and the equilibrium nematic order parameter S N = 0.65 (corresponding to the selected values of a, b and c). For this set of parameters, corresponding to the experimental parameters of 5CB, the simulations remain stable at the cost of choosing a very small time step. The simulation domain is L X × L Y = 112 × 56 for the simulations of a single skyrmion and L X × L Y = 200 × 200 for the simulations of 30 skyrmions, with periodic boundary conditions in both directions.
The initial configuration of each skyrmion is set as in Ref. [16] using an Ansatz for the director field: n x = sin(ã) sin mb + g n y = sin(ã) cos mb + g The parameter R controls the size of the skyrmion, B controls the width of the interface that separates the inner and the outer regions, m is the winding number of the skyrmion, g controls the direction of the skyrmion, r measures the distance from the skyrmion centre andb is the polar angle. We set the values of these parameters (in simulation units): m = 1, g = π/2, R = 0.7p, B = 0.5, C x = L X /2, C y = L Y /2. The velocity field is set to zero while the Ansatz is relaxed until it reaches the steady state.

Frank-Oseen free energy for quasi-static analysis
We utilize the Frank-Oseen free energy to study the 3D director structures of stretched torons and CF-2s. For chiral nematic LCs, the free energy can be expressed as where the Frank elastic constant K 11 , K 22 and K 33 determine the energy cost of splay, twist and bend deformations, respectively, and p is the cholesteric pitch. The electric field E is along z−axis ( E = U/d, with U being the voltage across the cell and d the cell thickness), ε 0 is the vacuum permittivity, and ∆ε is the dielectric anisotropy of the LC. Torons and CF-2s emerge as local or global minima of F F O , and a relaxation routine based on the variational method is used to identify energy-minimizing configurations. For elongated torons of elongated length, we prepare initial structures by stretching the point defects of the original toron and letting the whole structure to relax. The 3D computations are performed on 40×100×40 square grids with periodic boundary conditions in the x− and y−directions and with fixed boundary conditions at the bounding surfaces along the z−direction. For all calculations, the following values of the model parameters are used: d/p = 1.5, K 11 = 6.4 × 10 −12 N, K 22 = 3 × 10 −12 N, K 33 = 10 −12 N, U = 0.85V and ∆ε = 13.8.

Simulated polarized optical microscopy images
The POM image is simulated by the Jones-matrix method, using the energy-minimizing configurations of n(r) for the periodic CF-2s. The cell is split into 40 thin sublayers along the z direction, then we calculate the Jones matrix for each pixel in each sublayer by identifying the local optical axis and ordinary and extraordinary phase retardation. The phase retardation arises from the optical anisotropy of LC (n o = 1.58 and n e = 1.77 for 5CB), where the optical axis is aligned with the local molecular direction. The Jones matrix for the whole LC cell is obtained by multiplying all Jones matrices corresponding to each sublayer. We obtain the singlewavelength POM by the respective component of the product of the Jones matrix and the incident polarization. To properly reproduce the experiment POMs, we produced images separately for three different wavelengths spanning the entire visible spectrum (450, 550, and 650nm) and then superimposed them, according to light source intensities at the corresponding wavelengths.
Sample preparation and experimental methods Chiral LCs are prepared by mixing 4-Cyano-4'pentylbiphenyl (5CB, EM Chemicals) with a left-handed chiral additive, cholesterol pelargonate (Sigma-Aldrich). To define the pitch (p) of the ensuing chiral LCs, the weight fraction of the added chiral dopant is calculated using C dopant = 1/(h htp ), where the helical twisting power h htp = 6.25µm −1 for the cholesterol pelargonate additive.
The sample cells are assembled from indium-tinoxide (ITO)-coated glass slides treated with polyimide SE5661 (Nissan Chemicals) to obtain strong perpendicular (homeotropic) boundary conditions. The polyimide is applied to the surfaces by spin-coating at 2700rpm for 30s followed by baking (5 min at 90 • C and then 1 h at 180 • C).
The LC cell gap thickness is defined by silica spheres as spacers between two surfaces to be 7 µm with the cell gap to pitch ratio d/p = 1.25. The lateral sides of the LC cell are sealed by epoxy except for the two entrances. We use a homemade connection obtained by 3D-printing with the Formlabs Form 2 resin 3D printer (purchased from Formlab) to connect one entrance and tube, the other end of the tube being connected to a syringe. Then by pushing or pulling the plunger, the LC inside the cell channel is driven forwards or backwards. Metal wires were attached to ITO and connected to an external voltage supply (GFG-8216A, GW Instek) for electric control.
We utilize a ytterbium-doped fibre laser (YLR-10-1064, IPG Photonics, operating at 1064nm) to generate torons. The torons are three dimensional topological solitons where the skyrmion tube is embedded between the two substrates along the perpendicular axis, terminated by two point defects. We use a sinusoidal electric field with the frequency of 1000Hz. Initially, we adjust the voltage across the cell to be around U ≈ 3.5V where the LC is uniform and unwound. Then we melt the LC locally using a laser power around 30mW and switch off the laser tweezers, with torons being spontaneously generated after the LC quenching back. As we drive the LC flow forwards or backwards, the torons follow the flow.
Polarizing optical microscopy images are obtained with a multi-modal imaging setup built around an IX-81 Olympus inverted microscope and charge-coupled device cameras (Grasshopper, Point Grey Research). Olympus objectives 20x and 10x with numerical aperture NA=0.4 and 0.1 are used. The speed, length, and trajectories of the torons are analysed by using ImageJ (freeware from NIH).
Three-dimensional nonlinear optical imaging 3D nonlinear optical imaging of periodic CF-2 structures is performed using the three-photon excitation fluorescence polarizing microscopy (3PEF-PM) setup built around the IX-81 Olympus inverted optical microscope integrated with the ytterbium-doped fibre laser [31]. We use a Ti-Sapphire oscillator (Chameleon Ultra II; Coherent) operating at 870nm with 140f s pulses at an 80 MHz repetition rate as the source of the laser excitation light. An oil-immersion 40× objective (NA = 0.75) is used to collect the fluorescence signal. The fluorescence signal is detected by a photomultiplier tube (H5784-20, Hamamatsu) after a 417/60nm bandpass filter. The 3PEF-PM imaging involves a third-order nonlinear process, 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.