Dynamics of necklace beams in nonlinear colloidal suspensions

Recently, we have predicted that the modulation instability of optical vortex solitons propagating in nonlinear colloidal suspensions with exponential saturable nonlinearity leads to formation of necklace beams (NBs). Here, we investigate the dynamics of NB formation and propagation, and demonstrate a variety of optical beam structures emerging upon vortex beam propagation in engineered nonlinear colloidal medium. In particular, we show that the distance at which the NB is formed depends on the input power of the vortex beam. Moreover, we show that the NB trajectories are not necessarily tangent to the initial vortex ring, and that their velocities have components stemming both from the beam diffraction and from the beam orbital angular momentum. We also demonstrate the generation of elliptical rotating solitons and analyze the influence of losses on their propagation. Finally, we investigate the conservation of the orbital angular momentum in necklace and elliptical rotating beams. Our studies, performed in ideal lossless media and in realistic colloidal suspensions with losses, provide a detailed description of NB dynamics, and may be useful in analysis of light propagation in highly scattering colloids and biological samples.

charge and the radius fixed. In the first part, we show how the distance at which the instability onsets depends on the initial vortex power. Moreover, we study the trajectories of the solitons forming the NB and analyze the evolution of the NB profile for fixed power, depending on the charge of the vortex and its radius. We demonstrate the generation of elliptical rotating solitons resulting from the fusion of two solitons resulting from an interference of two vortices. Finally, we investigate the OAM conservation and the role played by losses in all the phenomena described above.
We study the propagation of beams carrying the OAM in a CS built of air bubbles with refractive index n p = 1 uniformly distributed in water with the refractive index n b = 1.33, resulting in a negative polarizability of the bubble. We note that in the laboratory experiments, one can use polytetrafluoroethylene (PTFE) particles dispersed in a glycerin-water solution to realize a stable negative polarizability CS. The volumetric filing fraction of this solution is f = 1‰ and the radius of the air bubbles is taken as r p = 50 nm. The suspension is assumed to be at the room temperature T = 293 K and the incoming free-space wavelength of light is λ 0 = 2π/k 0 = 532 nm. The nonlinearity in such a system has exponentially saturable character 40 and the propagation of the slowly varying envelope of the electric field φ can be described by: is the volume of the particle, and k B is the Boltzmann constant. The scattering cross section responsible for loss in the Rayleigh regime is given by σ π where ε 0 denotes the vacuum permittivity.
First, we study the dependence of the distance at which the MI onsets as a function of the power of the vortex with the profile of a soliton originating from the Laguerre-Gaussian profile with a fixed charge m = 2 and radius R = 20 μm (measured from the center of the vortex to the location of the peak intensity). Typical iso-intensity surfaces obtained during a 40-mm-long propagation are shown in Fig. 1(a)-(c). For all the iso-intensity plots in this Letter, the iso-intensity surfaces are drawn at the value of = I max{I(x, y, z)}/8 iso . Figure 1(b) shows the case where the input beam is a stable vortex solution. In the first 15 mm of propagation, we can see a stationary propagation of the vortex, during which the transverse profile does not change. After that, the MI onsets due to the presence of the numerical (discretization) noise, and four solitons are generated which travel away from the vortex ring. At the power levels lower than that for the stable solution (<8 W), the diffraction of the vortex beam dominates. As it can be seen in Fig. 1(a), for the power level of 6 W, in the first stage of the propagation, the vortex beam diffracts but the MI still manages to create the NB in a later stage of the propagation. Below this value of power, the MI is not observed and the vortex experiences only linear diffraction.
At the power levels higher than 8 W, the MI onsets at a longer distance than for the stable vortex soliton because of the saturation of the nonlinear effects in the CS. For a perturbation with a given modulation amplitude, the MI-mediated growth is slower in the regime close to the saturation than below the saturation regime, due to different slopes of the intensity dependent refractive index curve, n(I), in these two regimes (for more details, see the Supplementary Materials). Figure 1(c) illustrates the evolution of the vortex beam at a high-power level, close to the nonlinearity saturation. At first, due to the excess of power, the radius of the vortex soliton oscillates, which resembles the behavior of higher-order solitons. Then, the MI creates a NB built of four separate solitons. However, due to the excess of power in each of the solitons and the presence of the OAM in the beam, the created solitons have an internal twisted structure, that reflects the presence of two closely spaced solitons rotating around a common center of mass.
The plot in Fig. 1(d) summarizes the results of the study of the distance at which the MI onsets in function of the input power. The studies in which no noise was added to the input beam show a clear dependence of the MI-onset distance on the input power. The MI onsets at the shortest distance for the power corresponding to the stable vortex solution. Deviation from this power in any direction results in increase of the MI-onset distance either due to competition with diffraction (lower powers) or saturation effects (higher powers). Addition of a white noise with amplitude of 5% to the input beam accelerates the MI-mediated formation of NBs and leads to a weaker dependence of the MI-onset distance on the input power. Nonetheless, the character of this dependence is the same as in the noiseless case, as it can be seen in Fig. 1(d).
The presence of losses in the system prevents the NB formation at low-power levels but does not dramatically change the MI-onset distance for high-power levels. In Fig. 2, we present a comparison of iso-intensity surfaces for the high-power beams propagating in systems with and without loss. Loss strongly limits the propagation distance of the solitons created via the MI. Increase of the input power results in increased propagation distance of the created solitons due to the nonlinear self-induced transparency in the CS studied here 43 . With the increase of power, the loss experienced by the beam decreases, as the scattering particles with negative polarizability are expelled from the high-intensity region.
The number of beams into which the vortex splits due to the MI is not affected by the input power, for moderate levels of power, as it can be seen in Fig. 1(a)-(c). For powers of 20 W and 30 W, corresponding to more than 2.5 times the power of the stable soliton, the number of beams is reduced to three, and these beams have a strongly twisted character, as illustrated in Fig. 2(a) and (b).
Next, we investigate the trajectories and escape velocities of the solitons forming the NBs. Theoretical models have predicted that the escape velocity of a soliton v ⊥ , defined as the distance traveled in the transverse plane (x-y) divided by the distance traveled along the propagation axis (z), is proportional to the absolute value of the charge of the vortex and inversely proportional to the input vortex radius: 27 . This prediction is valid for vortex beams with radius R much larger than the width of the vortex donut [defined by the parameter ξ in eq. 2] 23 . Additionally, up to now, the trajectories have always been reported to be tangent to the NB ring 26,27 . Here, we analyze the validity of this analytical prediction for various radii of the vortex beams, and show that the trajectories can form an arbitrary angle with the line tangent to the NB ring.
Depending on the radius, the power level, and the charge of a vortex, the MI leads to formation of NBs with a different number of maxima [27][28][29][30] . Here, we analyze beams where each of the beam parameters can be controlled independently. To achieve that, instead of analyzing the MI-mediated breakup of the input vortex, we start from a NB with a predetermined number of maxima and fixed charge, radius, and input power. Such NBs can be generated using two methods (i) interference of two vortices with different charges but the same radii 22-24 resulting in a continuous distribution of phase in the azimuthal direction or (ii) arrangement of N solitons into a circular pattern where the initial phase of the jth soliton is given by 2πmj/N. Here, m denotes the NB charge and N is the number of solitons building the NB. This method results in a step-wise distribution of the phase in the azimuthal direction. Despite the fact that the latter method of NB generation creates more stable beams [50][51][52][53] , it is more cumbersome from the experimental perspective, as it requires generation of N independent solitons with appropriate phases. Therefore, here we use the former technique, which requires only two independent beams regardless of the number of beams building the NB. The input profile with N = 8 maxima and charge m is described in cylindrical coordinates (r, θ, z) by:   Fig. 3(d)], the contribution of diffraction is negligible compared to the OAM-induced velocity, and the beams travel along the trajectories quasi-tangent to the NB ring. We have performed similar studies for NBs with radii R = 100 μm and R = 500 μm, and the results are summarized in Fig. 3(e). As expected, for beams with larger radius, the contribution of diffraction is smaller and it is visible only for low charges m of the beam. For radius R = 500 μm the contribution of diffraction is negligible for the whole range of charges studied and the theoretically predicted relation | | =  they propagate. Generation of stable spiraling solitons was theoretically and experimentally demonstrated in materials with nonlocal nonlinearities 8,9,[11][12][13][15][16][17][18][19] . Since the nonlinearity in our model is assumed to be local and instantaneous, the spiraling beams are unstable in such a medium 54 . Above the critical radius, the two beams escape from each other; and below this radius, the beams fuse. This latter phenomenon generates elliptical rotating beams: corkscrew-like beams where two tightly spaced solitons rotate around the common center of mass.  Figure 4 shows the propagation of elliptical rotating beams in lossless and lossy media for various input-power levels. In the lossless case, at low power, the elliptical rotating beam is formed at the beginning of the propagation, but after approximately 20 mm, the rotation is no longer visible and the beam propagates as a regular soliton. The excess of energy, which is a result of the fusion of the two beams, is expelled away from the main fundamental soliton beam and only a single stable soliton remains. For higher input energies, the twisted nature of the beam is preserved for a longer propagation distance, as it takes more time to expel the energy excess. Increase of the energy does not allow for formation of higher-order solitons because the required nonlinear index modification exceeds the maximum attainable index change in the saturable CS studied here. In the case of propagation in the lossy medium, the elliptical rotating beam maintains the twisted character and propagates for longer distance as the power is increased, due to the self-induced transparency effect in the negative polarizability CSs 43 .
Finally, we have studied the conservation of the OAM of the NBs and the elliptical rotating beams. The total OAM is calculated as: . Figure 5 presents plots of the total OAM as a function of the propagation distance for four beams studied above: two NBs and two elliptical rotating beams. Blue and red curves show M z ( ) tot for NBs without and with losses, respectively. In the lossless case, for which the iso-intensity surface is shown in Fig. 1(b), the OAM is conserved both when the vortex beam propagates and transforms into the NB, and during the propagation of the generated solitons. In the case of the NB generation in a lossy medium [shown in Fig. 2(d)], the OAM is not fully conserved. It is conserved during the propagation of the vortex beam, but when the NB is formed, part of the OAM is lost. After the solitons lose energy and diffract, the OAM does not decrease any more. The loss of the OAM in the middle stage of the propagation is caused by fact that the three solitons forming the NB do not have equal energies, as they grow from random noise due to the MI. In our model, loss is intensity dependent, and the solitons with the lower power lose the energy faster than the solitons with higher power. This results in a non-symmetric distribution of the energy and leads to changes in the OAM carried by the NB. When beam intensities decrease, all the beams decay at the same rate (corresponding to the loss level in the linear regime), and the OAM of the linearly diffracting field becomes conserved again.
For the elliptical rotating beams resulting from the fusion of two solitons, the OAM is not conserved, as shown by the orange and violet curves Fig. 5. Here, the OAM is expelled from the elliptical rotating high-intensity beam and it is carried by the low-intensity field (not visible in iso-intensity plots) away from the central beam. The confirmation of this hypothesis can be seen comparing the plots of the OAM and the power evolution, shown in the videos provided in the Supplementary Materials. During the lossless propagation shown in Fig. 4(c), the low-intensity field carrying the OAM is absorbed when it reaches the simulation domain boundaries. Surprisingly, the twisted character of the beam is present, despite the fact that the beam does no longer carry the OAM. In the system with losses [see Fig. 4(e) and the corresponding violet curve in Fig. 5], the low-intensity field that carries the OAM is rapidly scattered by the CS. As a result, the OAM in this system decreases to zero in the first 10 mm of propagation.
In conclusion, we have investigated the dynamics of vortices in nonlinear colloidal suspensions and showed that the distance at which the modulation instability transforms a vortex beam into a necklace beam depends on the input power, and that the minimum distance is achieved for the power corresponding to the stationary vortex soliton. Moreover, we have demonstrated the influence of diffraction on the trajectories of necklace beams and the value of the escape velocity. The solitons forming the necklace beam have been shown to follow linear trajectories that are not necessarily tangent to the necklace beam ring, and the angle of escape can be controlled by the interplay of the necklace beam radius and its charge. We have presented the possibility of elliptical rotating soliton formation using necklace beams with two maxima. Finally, we have analyzed the influence of losses and the conservation of orbital angular momentum in the necklace beams and elliptical rotating beams. These results may be of interest for the formation of complex photonic structures in liquids as well as in the broader field of light filamentation. Data Availability. No datasets were generated or analyzed during the current study. Figure 5. The total OAM as a function of propagation distance for necklace beams (blue and red curves) and elliptical rotating beams (orange and violet curves) shown in the plots indicated in the legend. OAM conservation during propagation in colloidal media without loss (blue and orange curves), and with loss (red and violet curves) is analyzed.