Kinetics of motile solitons in nematic liquid crystals

The generation of spatially localized, soliton-like hydrodynamic disturbances in microscale fluidic systems is an intriguing challenge. Herein, we introduce nonequilibrium solitons in nematic liquid crystals stimulated by an electric field. These dynamic solitons are robust as long as the electric field is maintained. Interestingly, their kinetic behaviours depend on the field condition—Tuning of the amplitude and frequency of the applied electric field alters the solitons to self-assemble into lattice ordering like physical particles or to command them to various dynamic states. Our key property to the realisation is the electrohydrodynamic instability due to the coupling between the fluid elasticity and the background convection. This paper describes a new mechanism for realising dynamic solitons in fluid systems on the basis of the electrohydrodynamic phenomena. Electric field induced convective instabilities in liquid crystal slabs can assume a localized shape. The authors show how to generate, manipulate and shepard these dissipative solitary excitations that do not require topological constraints for stabilization.

S olitons are conventionally known as spatially localised, shape-preserving travelling wave packets. Since the first observations of the solitonic behaviour of water waves by John S. Russell in 1834 and the introduction of the term of 'soliton' by Zabusky and Kruskal 1 , solitons have been identified in various natural phenomena. In theoretical frameworks such as the Korteweg-de Vries (KdV) or the nonlinear Schrödinger equation 2 , solitons appear as particle-like wave packets that travel without shape change and dissipation by balancing dispersive and nonlinear effects. This ideal behaviour of the solitons would enable the distortion-free long-distance transport of waves or structures, which has attracted considerable interest from both fundamental and technological viewpoints in many branches of science, including atmospheric circulation 3,4 , fluid 5 and mechanical waves [6][7][8] , living organisms 9,10 , Bose-Einstein condensates [11][12][13] , light confinement, and propagation [14][15][16] .
Very recently, there have been several groundbreaking studies observing solitonic analogues in polar materials, such as spin systems [17][18][19] , ferromagnetic fluids 20 , and liquid crystals [21][22][23][24][25][26] . In most cases, these solitons possess topological singularities in spin-like vector fields that are distinct from the far-field background. In other words, the structures of the solitons are topologically protected. Mochizuki demonstrated a unidirectional thermally driven ratchet motion of magnetic skyrmions in the presence of a temperature gradient 17 . Smalyukh et al. experimentally and theoretically revealed that the twisted orientational order of chiral nematic liquid crystals under confinement induces assorted topological solitonic structures that can squirm in a specific direction under an electric field [20][21][22] . These topologically protected structures are static in thermal equilibrium but motile under an inhomogeneous electric field. Of course, since there are lots of examples of static topological structures known in many physical systems, these experimental facts do not guarantee all topologically protected structures to become motile. On the other hand, there have been reported intrinsically dynamic non-topologically protected solitons, where the KdV theory is applicable. As an example of such nontrivial solitons, Li et al. have recently observed travelling solitonic waves, which proposes a novel formation mechanism based on the flexoelectric effect in a non-chiral nematic liquid crystal 23,25,26 . While these previous works, including the ones published by Li et al. 25 and Sohn et al. 26 during the submission of the present paper, considered the mechanisms mainly based on the localised solitons without the electrohydrodynamic instability, still it is well known that electrohydrodynamic instability plays crucial roles in the formation mechanisms and dynamic behaviours of localised structures and patterns in many nonequilibrium fluid systems [27][28][29] . Therefore, employing the electrohydrodynamics is a fascinating strategy for exploration of dynamic solitons.
In this article, we show how solitons can arise in a parameter regime where the Carr-Helfrich electrohydrodynamic instabilities would be expected to occur. We even tune the anisotropies of both the dielectricity and the conductivity to generate a smooth path from standard flexoelectric or electrohydrodynamic states to dynamic solitonic states. The localised non-equilibrium excitation assumes the character of quasi-particles and enables various kinetics, e.g. directional translation, collision, reflection, and proliferation. The study establishes the fundamental physics and strategy of manipulating the kinetics of solitons in more general experimental settings and permits the observation of their time evolution, facilitating further experimental and theoretical studies of dynamic soliton excitations in soft-matter systems.

Results
Soliton induction. In contrast to refs. 23,25 , we consider more general conditions accompanied by electrohydrodynamics, where dielectricity and ionic conductivity play important roles. In order to systematically examine the effects of the material parameters, we use the frustrative mixtures comprising two nematics, E7 (a mixture of cyanobiphenyl and cyanoterphenyl components) and 4ʹ-butyl-4-heptyl-bicyclohexyl-4-carbonitrile (CCN47) (Fig. 1a) from the beginning. By altering the concentration of E7 or adding excesses of the ionic species tetrabutylammonium benzoate (TBABE) to the CCN47, we can facilitate a continuous transition of the anisotropic dielectric/ conductivity states from (−−) to (++), where the pairs of the signs, (ΔεΔσ) (the so-called de Gennes' notation 28 ), are given corresponding to the signs of the dielectric and conductivity anisotropies, Δε ¼ ε k À ε ? and Δσ ¼ σ k À σ ? . The subscripts || and ? indicate the directional components parallel and perpendicular to the average molecular orientation, denoted by the nematic director n, respectively. In particular, we pass the (−+) state, characteristic for the electrohydrodynamic instability of the Carr-Helfrich phenomenology.
Through extensive experiments, we find that solitons can be created in the frequency range of 2-60 Hz and for mixtures of E7-CCN47 with 15-18 wt% E7 (see also the state diagrams in Fig. 2f-h). As seen in the generalised scheme of the procedural conditions for the solitons (Fig. 2g), stable solitons are observed only in (−+) at low electric fields (≤10 V), that is, a typical condition in which the Carr-Helfrich electrohydrodynamics becomes active [27][28][29] . In this case, a convective flow occurs parallel to the background nematic director n, while in contrast, no observable flow occurs for (−−) and (++) accompanied with the flexoelectric effect and Freederikz effect-driven uniform reorientation, respectively. The present mixture system does not exhibit (+−), but it is known in the conventional systems that the Freederikz effect dominates in (+−).
As expected, conductivity is vital in determining the stability of the solitons in the present system, i.e. the solitons appear in the limited range of the moderate conductivity, 8 10 À9 < σ < 4 10 À8 Ω À1 m À1 (Fig. 2h). However, this conductivity value is much lower than those required for the conventional Carr-Helfrich electrohydrodynamic phenomena. For example, even in the (+−) state of almost the same E7-CCN47 mixture with approximately 16 wt% E7, there emerges only the global flow pattern without isolated solitons if too many ions are added and the conductivity exceeds the above-mentioned range. This means, the conventional convection regime upon the Carr-Helfrich electrohydrodynamic effect occurs everywhere in the material. On the other hand, in the case of too few ions and hence too low conductivity, the electrohydrodynamic effect is suppressed and the Freederikz effect dominates over the pattern formation (see Supplementary Note 11). The similar low-conductivity condition has been tested by Li et al. for the (−−) state, where the flexoelectric effect works instead of both the Carr-Helfrich electrohydrodynamic and the Freederikz effects 23 . One possible physical situation deduced from these facts is that some localised effect on the top of the Carr-Helfrich electrohydrodynamic mechanism is happening for the emergence of the present solitons: the solitons are a superposition of the electrohydrodynamic rolls, which can coexist with the global electrohydrodynamic roll patterns (states (II), (III), and (V)) as discussed later (see also Supplementary Note 8 and Supplementary Movie 12). Indeed, as evidenced in Supplementary Notes 1 and 2, and Supplementary Figs. 1 and 2, the ions are considered to be localised to produce independent solitons, compared to the relatively uniform spatial distribution of the ions for the global flow pattern based on the conventional Carr-Helfrich electrohydrodynamics. In such ion-rich areas, the electrohydrodynamic flow can easily occur due to the larger hydrodynamic torque 29 . The soliton radius, R, is defined as the radius of the primary area undergoing elastic deformation. The space-filling ratio X is defined as the ratio of the total area of N solitons to the sample area, A. Image width, 65 µm. g Counterpart to f taken with a tint plate. h Monochromatic textures with a λ/4 plate and schematics of the alternating director field between two orientational states, taken with a high-frame-rate camera. Image width, 30 µm.
In Fig. 1b-d and Supplementary Movie 1, we show timedependent microscopic texture images acquired during soliton creation in a mixture containing 16 wt% E7 for a 10 V electric field at 12 Hz at the specific times designated in Fig. 1e. Independent cruciform solitons are randomly excited in space and time (Fig. 1b, c), while they disappear when the field is absent ( Supplementary Fig. 1). These solitons fluctuate around and remain relatively near their induction positions with a space-filling ratio of X = πR 2 /A < 0.6, where R and A are the soliton radius and the total area of the region, respectively. Since the time dependence of X reflects how the involved ions localise, it is fitted by the Kolmogorov-Johnson-Mehl-Avrami expression, which might give insights into the procedure of the migration of the ions during the soliton generation (Supplementary Notes 1 and 2, where X, t, K, and m are the space occupancy of the directordeformed region, the time, the temperature-dependent Avrami coefficient, and the exponent, respectively [30][31][32] . The value obtained for m is 2.9, suggesting the involvement of two-dimensional homogeneous nucleation and growth processes 33  packing structures occur in some self-driven systems 23,25,34,35 . Next, we investigate the director structure of each soliton. Under an exposure time sufficiently longer than the oscillatory period of the applied electric field, the dynamic structural properties of the solitons cannot be captured and only a time-averaged structure is clarified. Figure 1f, g show magnified views of the time-averaged solitonic structures obtained with an exposure time of 200 ms under crossed polarisers and with a standard tint waveplate (retardation equal to 530 nm) taken by a conventional camera. An array of the concentric patterns of the time-averaged solitons pretends as if s = +1 topological defects ordering in the midplane of the nematic film. In order to elucidate the dynamic structure, we perform high-frame-rate polarising microscopy at 100 frames per second. Figure 1h shows a time-dependent oscillation of the structure of the solitons. It becomes clear that the director field of the solitons transforms between two defectless structures, i.e. biconvex and biconcave structures, synchronised with the applied electric field. It is noted that, by analysing the director fields, the biconvex structure has a larger deviation angle of the director in the midplane of the nematic film than the biconcave structure (time-resolved structures in Supplementary Fig. 5). This explains why the time-averaged structure looks like a concentric domain with a s = +1 defect.
As the soliton-filling ratio increases, mutual repulsion, separating the solitons by distances of about twice the soliton radius, is observed (Supplementary Movie 1). As a result of this repulsion, the solitons become resistant to the external stress without coalescence, confirming the pseudoparticle nature of the solitons.
Conditions of soliton creation. The orientational pattern that includes the solitons varies depending on the voltage and frequency ( Fig. 2a-e). On the basis of the director deformations and flow profiles, we classify the frustrated states into five categories: (I) ground planar alignment (G-state), (II) ||-roll pattern (||-state), (III) rectangular pattern (R-state), (IV) soliton state (Soliton-state), and (V) ?-roll pattern (?-state). Figure 2f shows a state diagram for a mixture containing 16 wt% E7, which exhibits (−+) anisotropy.
In the G-state, neither elastic deformation nor hydrodynamic effect occurs; therefore, the initial uniform uniaxial alignment is maintained. In the ||-state, which appears in the lower-frequency region <9 Hz, only an in-plane splay elastic deformation of the director occurs because of the coupling with the flexion-induced polarisation, the so-called flexoelectric effect 27,36 Because the strong surface-anchoring effect forces fix the surface directors in the initial state, twist deformation along the electric field direction (z direction) is induced, yielding a one-dimensional periodic pattern. Particle tracking confirms the absence of flow (see 'Methods'). An increase in the voltage in the ||-state causes a fieldinduced hydrodynamic convective flow parallel to the initial director ( n 0 ). The coupling between the flow and the flexoelectric effect induces a two-dimensional version of the ||-state, i.e. the Rstate. The ?-state occurs in a dielectric regime in which ions are effectively immobile, but the nematic director oscillates with the applied electric field. The convective flow field parallel to n 0 arising from an electrohydrodynamic effect (Carr-Helfrich instability) yields a stripe pattern in which the director has a wavy modulation in space, as observed in nematics with (−+) anisotropy 29 . The appearances of the respective roll patterns in the ||-state and ?-state significantly differ in terms of the roll direction: rolls in the ||-state run parallel to the director; those in the ?-state are perpendicular. The Soliton-state occurs in the frequency range between the ||-/R-state and ?-state and spans a conductive regime in which ionic convection dominates the dielectric oscillation of the director. Accordingly, the soliton system has features of both the ||-state and ?-state with a compresence of elastic and convective effects, as discussed below. The way the solitons are observed to move in our system depends on frequency as well as on the amplitude of the voltage applied. For instance, at constant frequency but with the increased voltage of the electric field (and thus the director oscillation), a transition from a motional pattern coined 'swimming' to 'proliferation' is observed. Worth noting that voltage or frequency contributes differently to the kinetics of the solitons as discussed below.
Collision and elastic reflection in the swimming regime. An increase in the voltage or frequency in the oscillation regime (the transition from the orange to purple-blue regions in Supplementary Fig. 7) produces steady directional motion of the solitons in which the initial time-averaged cruciform configuration changes to a time-averaged bug-eye-like pattern because of asymmetric elastic deformation (Fig. 3a, b). In this regime, the director structure of the solitons during the motion is also dynamically oscillating with time ( Supplementary Fig. 8). The direction of the motion can be tuned by changing either the amplitude or frequency of the electric field.  Fig. 3c), the solitons move directionally parallel to n 0 in the xy plane. Notably, they randomly move to the left or right because of the previously discussed dielectric oscillation. As plotted in Fig. 3i, at V > 8.5 V, the directional angle of motion with respect to n 0 continuously varies as the voltage is increased (Fig. 3d-h, Supplementary Movies 3-8), corresponding to the increasingly asymmetric deformation of the bug-eye-like pattern (Fig. 3b). High-speed camera observations reveal that this movement is driven by the oscillation of the director field between distorted versions of two orientational states, as shown in Fig. 1h (Supplementary  Movie 9). The directional motion ends at about 10 V, after which the soliton motion becomes somewhat chaotic, leading to a fractalised proliferation process that will be discussed later. The average velocity of the motion of the solitons is measured through a trajectory tracking as shown in Fig. 3c-h and plotted as a function of voltage and frequency in Fig. 3j, k. Although the velocity is nearly constant with respect to the voltage because the background flow field is mostly unaffected in the limited voltage range, a significant frequency dependence is observed. The solitons are resonantly active at about 20 Hz, as shown in Fig. 3k. This frequency corresponds to a critical frequency that separates the conductive and dielectric regimes 29 , as described by the competition between the relaxation process of mobile charge in space and the dielectric relaxation of the director, Given σ ⊥ = 10 −9 Ω −1 m −1 and ε ? = 4, as measured experimentally, we obtain the critical frequency of 28 Hz by calculation, consistent with the experimentally observed value. The dynamics of the collisions between solitons and their reflection by obstacles can be understood in terms of their pseudoparticle nature. Figure 4a-c and Supplementary Movie 10 show that their head-on collisional behaviour varies depending on the degree of offset, δ, between soliton cores. When R < δ < 2R, interference occurs between the peripheries of soliton pairs with the same orientation. Consequently, the solitons repel each other at a small post-collision angle; at 0 < δ < R, elastic deformation becomes more significant because of the increased mismatch of the director field between solitons, leading to vertical repulsion between the solitons. Following this collision, the solitons separate such that elastic deformation is avoided and then continue to move along their original directions of motion. Interestingly, the solitons are elastically reflected at an air-nematic interface in the manner of real particles colliding with a wall (Fig. 4d).
Proliferation regime. Figure 5a, b show the fractalised proliferation trajectory of destabilised solitons for a 10 V field at 20 Hz (Supplementary Movie 11). Soliton fractalisation produces clones for movement perpendicular to n 0 . Upon fractalisation, the solitons bifurcate and move obliquely for a short distance. Their dynamic structure is preserved during this process, and the emerged two solitons have an opposite phase in their orientational structure. The fractalisation process corresponds to the continuous growth in soliton size. Considering that an optimum size exists for the solitons at a fixed film thickness ( Supplementary Fig. 11), the fractalisation process results in the continuous accumulation of excess elastic energy. The solitons fractalise when the energy penalty from elastic deformation cannot be compensated by the effects of surface anchoring and dielectric interaction. This is illustrated in Fig. 5c, which shows the soliton diameter, 2R, as a function of time; fractalisation occurs when the soliton size exceeds a specific value. Notably, the fractalisation simultaneously occurs for all solitons (Fig. 5a, b). This implies that a synchronised growth mechanism is involved for the solitons, resulting in a solitonic trajectory with a Hausdorff fractal dimension of D = 1.

Discussion
The localised solitary oscillating director structures can also exist under conditions for global flow patterns to occur; there is a coexistence of states (IV), Soliton-state, and (V), the ?-state (see consist of two oblique waves travelling at an angle θ with respect to n 0 , as introduced elsewhere 37-39 . As a first-order approximation, the amplitude expansion of the flow field is, where A and B are the respective spatiotemporal variable amplitudes of the two assumed travelling waves; k A , k B , and k C are the wavenumbers; r is the spatial coordinate; t is time; and ω is the angular frequency. For small amplitudes where Ω is the angular frequency, the description of the flow waves in the Ginzburg-Landau framework is valid, yielding and where u A,B are the group velocities of the waves, μ is the bifurcation parameter, and c and g are coupling coefficients. Forming a symmetric superposition mode of the two waves in the xy plane with a finite Ω proportional to A 0 j j 2 produces a two-dimensional localised travelling rectangle solution corresponding to a solitonic state. If only the flow field is under action, two types of solitonic motion are possible: a directional motion parallel to n 0 and a twodimensional random walk within the xy plane. However, in our results, variable angles of solitonic motion with respect to n 0 are achieved. This motivates us to consider the elastic deformation of the bulk directors by introducing a spatially modulated imbalance of amplitudes between A and B. Although convective flow is considered in ref. 37 , our solitons have splay deformation in the xy plane and twist-bend deformation in the xz plane (Fig. 2d) in addition to a convective flow. Owing to the flow-alignment property of the director, its orientation tends to be parallel to the flow field. Since the orientation of the director oscillates with time, the local flow also changes its direction accordingly through a self-regulating procedure ( Supplementary Fig. 12). Such processes can trigger the local motion of topological defects 40 , thereby providing an additional travelling mode out of the xz plane. Indeed, the electrohydrodynamic convection flow can be modulated by twist deformation to produce a spiral flow 41 . Because the elastic deformation in the xy plane weakens in the transition from the ||-state through the Soliton-state to the ?-roll state as the frequency is increased or the voltage is decreased, the inclination of the oblique solitonic movement decreases with respect to n 0 ; with neither flow nor elastic deformation in the xy plane, no two-dimensional motion is observed in the studied parameters (Fig. 2a-e). Thus the local coupling between the elastic deformation and the flow field in an anisotropic medium permits effective dictation of the solitonic motion directionality, which cannot be achieved through a flow field alone.
Next, we discuss the pseudoparticle properties of solitons in anisotropic flows. For characterisation, we analyse a radial distribution function of packed solitons, g(r cc ) ( Supplementary  Fig. 13) and the thermal energy-activated solitonic Brownian diffusion. Figure 6a shows the pairwise potential V(r cc ) as a function of the centre-to-centre distance between solitons, r cc . At the dilute limit, V(r cc ) is directly calculated as Àk B T log g r cc ð Þ ½ . The first potential closely corresponds to the mean nearestneighbour distance between solitons, r cc,o = 14.2 µm. The depth of the potential well slightly exceeds the thermal energy, k B T, which is consistent with the observed minimal soliton fluctuation beyond r cc,o . The solitonic Brownian diffusion process is calculated as the two-dimensional mean square displacement of the solitons over time t (Fig. 6b). The effective diffusion constant of the solitons is 1.4 µm 2 s −1 , corresponding to a mean solitonic fluctuation velocity of about 1 µm s −1 . The effective mass of the solitons, m eff , is expressed as where U is the kinetic energy and v the mean velocity of the solitonic fluctuation. By balancing the thermal energy, k B T, and the kinetic energy of the solitons, m eff is estimated to be 1 picogram. In contrast, the real mass of molecules constituting each soliton is estimated to be ρπR 2 h ≈ 1000 picograms, where ρ and h are the material density and nematic film thickness, respectively. The effective mass of the soliton is significantly smaller than the real mass of the materials constituting the solitons. Notably, m eff is independent of the film thickness irrespective of the change in the real mass ( Supplementary Fig. 14). Since the present soliton is a localised deformation of the director field as shown in Fig. 2d, such a time-averaged three-dimensional state is approximated using a spherical distortion, where R and d are radius of the solitons and the film thickness with a relationship of R ∝ d as shown in Supplementary Fig. 11. The elastic energy of the solitons is calculated by considering the Frank elasticity: with different film thicknesses, by using experimentally determined splay (K 11 = 13 pN), twist (K 22 = 6.2 pN), and bend elastic constants (K 33 = 15 pN). Therefore, the decrease of elastic energy density upon increasing the film thickness is compensated by the increase of the volume of the soliton, leading to the invariance of the total elastic energy ( Supplementary Fig. 15). This implies that the constant m eff reflects the elastic energy stored by the elastic deformation in the director field.
In summary, we have demonstrated the generation of nonequilibrium solitons and the control of their kinetics in nematic liquid crystals. Unlike the well-known Russell's solitons that propagate in isolation after their excitation, the present nonequilibrium solitons are structurally robust upon application of an electric field and dissipate after the field is turned off. It is revealed that the stabilisation of the solitons is rooted in the socalled Carr-Helfrich electrohydrodynamic mechanism, where the interactions between background flow and elastic deformation play important roles. Interestingly, these solitons show unique dynamic behaviours, i.e. controllable directional swimming and proliferation. The results suggest a new possible pathway that  solitons can be made and manipulated in nonequilibrium conditions based on the electrohydrodynamic mechanism in anisotropic fluid systems.

Methods
Sample preparation. The frustrative nematic mixtures were formed by mixing E7 (Wako, Aldrich, and made in house) with CCN47 (Nematel GmbH) by weight. These molecules form a uniaxial nematic phase in which the average molecular orientation is characterised by a vectorial director, n. Because of the macroscopic and nonpolar nature of the nematic phase, the head and tail of n are equivalent. The liquid crystals are uniformly aligned between two glass plates with a controlled thickness in the range of 2-50 µm. In the manuscript, if not specified, we show solitonic behaviours in a 5.2-µm-thick film. Although E7 has positive anisotropies (ΔεΔσ) = (++), those of CCN47 are negative (ΔεΔσ) = (−−). TBABE was added to the E7-CCN47 mixture to tune its net strength of conductivity. Films comprising these mixtures were sandwiched between pairs of substrates coated with rubbed polyimide on indium-tin-oxide electrode layers, with the film thickness adjusted by adding silica particles to obtain a 2-12-µm-thick film (micromer, Micromod) or by adding polyester film spacers to obtain a 20-50-µm-thick film (Mylar, Dupont). To create solitons, a sinusoidal or rectangular AC electric field is applied normal to the plates (Fig. 1a). In many other liquid crystalline systems, including Schiff bases and cyanobiphenyl mesogens, solitons were also found.
Conductivity measurement. The electrical and dielectric responses in the liquid crystal mixtures were characterised via dielectric spectroscopy (Solartron Analytical, impedance/gain-phase analyser 1260A with dielectric interface system 1296A). An AC electric field at 0.05 V rms sweeping at a frequency of 10 −2 -10 4 Hz was applied to the films, and the resulting impedance response was recorded and analysed ( Supplementary Fig. 16).
Director mapping. The orientation of the directors was determined optically using polarising microscopy. The highest birefringence of the mixtures was sufficiently low (0.03 < Δn = n e − n o < 0.08, where n e and n o are the extraordinary and ordinary refractive indices, respectively) to carry out this task. On the basis of the measured transmittance of light through the films, we calculated the orientation of the directors. Our method was the same as that presented in ref. 20 . We also conducted fluorescence confocal polarising microscopy (TCS SP8 STED, Leica) to probe the three-dimensional distribution of the soliton directors. Although volume scanning was performed during soliton movement, owing to their fast fluctuation that was synchronised with the electric field, the solitons were not clear enough to allow their analysis. Instead, we measured the emission intensities at single planes on the surfaces of the slab and in the middle of the cell. The intensity on the surfaces was nearly constant during the application of the electric field, confirming a strong surface-anchoring effect. The signal was maximum at the middle of the slab, corresponding to elastic deformation of nearly zero on the surface that reached a maximum in the middle of the slab.
Characterisation of flow field. The particle tracking method was used to visualise the effective flow in the mixtures in the presence of an electric field. In this process, a small number of 2-µm particles were dispersed into the nematic mixtures, and their motions were tracked by both bright-field and polarising microscopy.
Time-dependent trajectory analysis. The recorded dynamics of the solitons were analysed using custom-made ImageJ macros. For the precise image analysis needed to extract sharp shapes and calculate the centres of gravity of the solitons, a dark but slightly structured time-dependent background was averaged over time.

Data availability
All data that support the findings in this study are available in the article and in Supplementary Materials. Additional information is available from the corresponding author upon reasonable request.