Generation and evolution of different terahertz singular beams from long gas-plasma filaments

We theoretically and numerically investigate the generation and evolution of different pulsed terahertz (THz) singular beams with an ultrabroad bandwidth (0.1–40 THz) in long gas-plasma filaments induced by a shaped two-color laser field, i.e., a vortex fundamental pulse (ω0) and a Gaussian second harmonic pulse (2ω0). Based on the unidirectional propagation model under group-velocity moving reference frame, the simulating results demonstrate that three different THz singular beams, including the THz necklace beams with a π-stepwise phase profile, the THz angular accelerating vortex beams (AAVBs) with nonlinear phase profile, and the THz vortex beams with linear phase profile, are generated. The THz necklace beams are generated first at millimeter-scale length. Then, with the increase of the filament length, THz AAVBs and THz vortex beams appear in turn almost periodically. Our calculations confirm that all these different THz singular beams result from the coherent superposition of the two collinear THz vortex beams with variable relative amplitudes and conjugated topological charges (TCs), i.e., +2 and −2. These two THz vortex beams could come from the two four-wave mixing (FWM) processes, respectively, i.e., ω0+ω0−2ω0→ωTHz and –(ω0+ω0)+ 2ω0→ωTHz. The evolution of the different THz singular beams depends on the combined effect of the pump ω0−2ω0 time delay and the separate, periodical, and helical plasma channels. And the TC sign of the generated THz singular beams can be easily controlled by changing the sign of the ω0−2ω0 time delay. We believe that these results will deepen the understanding of the THz singular beam generation mechanism and orbital angular momentum (OAM) conversion in laser induced gas-filamentation. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
In the last decade, the generation of many terahertz (THz) singular beams have created widespread interest across a wide range of fields, especially THz beams with orbital angular momentum (OAM) [1][2][3][4]. These demands for numerous applications have drawn considerable attention to different THz OAM beam generation for years [1,2]. However, there still have been very fewer methods for the generation of different THz OAM beams [1]. Very recently, much attention has been paid to THz OAM beam generation via vortex-laser-induced gas-plasma [5,6]. In specific, these are intensity modulated THz vortex waves, which were also called THz angular accelerating vortex beams (AAVBs) due to their phases evolving nonlinearly with azimuthal angle θ [7]. Such plasma-based THz emitters can avoid the limitation to material ionization threshold. However, in these studies mentioned above, propagation effects are basically neglected because of the short gas-plasma lengths and only one type of THz singular beams is generated and investigated. While the propagation effects have significant influences on the filament dynamics and the emitted THz spatial characteristics [8]. So, if propagation effects are taken into consideration in long gas-plasma filaments, what changes will we get?
In this paper, theoretical and numerical investigations on the generation and evolution of different pulsed THz singular beams with ultrabroad bandwidth (0.1-40 THz) in long gas-plasma filaments are carried out via performing a unidirectional propagation model under group-velocity moving reference frame. The THz radiation is excited via long gas-plasma filaments induced by a shaped two-color laser field, which includes a Laguerre-Gaussian (LG) fundamental pulse (ω 0 ) and Gaussian second harmonic pulse (2ω 0 ). Surprisingly, it is found that, as the filament length increases, three different THz singular beams occur, including the THz necklace beams (π-stepwise phase profile), the THz AAVBs (nonlinear phase profile) and the THz vortex beams (linear phase profile). At a millimeter-scale filament length, the THz necklace beams are generated firstly. Then, with the increase of the filament length, THz AAVBs and THz vortex beams appear in turn almost periodically. It is confirmed that the two four-wave mixing (FWM) processes, i.e., ω 0 +ω 0 −2ω 0 →ω THz and -(ω 0 +ω 0 ) + 2ω 0 →ω THz , take place with different probability amplitudes at different filament lengths, which contributes to the generation of the THz vortex beams with the conjugated topological charges (TCs), respectively. They interfere with each other to synthetize the different THz singular beams. Besides, it is demonstrated that the combined effect of the pump ω 0 −2ω 0 time delay and the separate, periodical and helical plasma channels plays a key role in the generation efficiencies of the two FWM processes. In addition, the sign of the time delay determines the dominated FWM process, thereby determining the TC sign of the generated THz singular beams. This work may provide a practical and efficient strategy to control local OAM in the THz region, thus to generate different THz singular beams. These findings will provide some theoretical supports for the OAM (de)multiplexing ultra-broadband THz communication, particle manipulation, plasma control, optically-driven micro-machines, trapping in-parallel molecular or cell assays, material processing and nonlinear optics [9,10].

Physical Model
In our numerical experiments, we consider a shaped two-color field consisting of a linearly polarized 800 nm LG fundamental pulse (ω 0 ) with its TC of l=1 and its Gaussian second harmonic pulse (2ω 0 ) with the same linear polarization focused by a lens with a focal length of f in nitrogen (N 2 ) under the standard atmospheric pressure. Without loss of generality, the original phase of Gaussian 2ω 0 is assumed as zero. With this regard, the phase difference between ω 0 and 2ω 0 on the azimuthal angle θ has the form of ∆φ=2lθ=2θ.
A numerical simulation of the laser field in gas-plasma filaments based on the unidirectional propagation model [11,12] is carried out to reveal the generation and evolution characteristics of the different THz singular beams. The spectral-domain unidirectional propagation equation in the group-velocity moving reference frame can be expressed as [11][12][13] whereÊ is the spectra of electric field. ∇ 2 ⊥ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the transverse Laplacian, (k (0) /n!)(ω − ω 0 ) n describes the high-order dispersion with k (0) =n 0 ω 0 /c and k (n) =[d n k(ω)/dω n ] ω=ω0 , while µ 0 is the vacuum permeability and n 0 is linear refractive index of the laser.F NL =P NL + iĴ e /ω + iĴ loss /ω is the nonlinearity term including the third-order Kerr nonlinear polarizationP NL , the current densityĴ e and a ionization energy loss termĴ loss . For clarity, here we have dropped the (x, y, z, ω) dependence of (Ê,F NL ,P NL ,Ĵ e ,Ĵ loss ) and the (x, y, z, t) dependence of (E, F NL , P NL , J e , J loss , N e ). And P NL =2ε 0 n 0 n 2 |E| 2 E, J e =(e 2 /m e )[(v c +iω)/(v c 2 +ω 2 )]N e E, J loss =R 1 (N at -N e )U i /E, where ε 0 and n 2 are the dielectric constant in vacuum and nonlinear refractive index of the laser, respectively. Generally, in an early self-focusing stage, self-focusing stage FWM is dominant; but when plasma occurs, photo-ionization prevails [14]. The FWM mechanism leads to a cos(∆φ) dependency, while the photocurrent (PC) mechanism leads to a sin(∆φ) dependency [14,15]. Basically, these dependencies are based on point-like source without propagation and dispersion influences. However, as for centimeter-scale long plasma filaments, propagation effects have to be taken into consideration, then it would not follow these simple forms of dependency relationship [16]. Recently, both experimental and numerical results have proved that, for long plasma filaments, both Kerr nonlinearity (FWM mechanism) and plasma current (PC mechanism) contribute to THz generation [8,17]. However, they usually work at different situations and contribute to different THz frequency ranges [8]. It is noteworthy that the critical power of a Gaussian pulse is given by P crit (0)≈λ 0 2 /[4πn(ω 0 )n 2 ] [18]. So the critical power for an 800 nm Gaussian pulse is about 1.27 GW. In contrast, an 800 nm LG pulse with a TC of l=1 has 4 times larger critical power as that of a Gaussian pulse, i.e., P crit (1) = 5.09 GW, which is much bigger than that of a Gaussian one. According to the Ref. [8], it is found that the Kerr nonlinearity is dominated over plasma current at moderate pulse intensity <100 TW/cm 2 (z <=4.5 cm), in which the generated THz spectrum predominates at tens of THz and is depleted near zero frequencies. Whereas the plasma current overwhelm the contributions from Kerr nonlinearity when the intensity reaches clamping intensity (>=150 TW/cm 2 ) (z >=5.7 cm). These results are under the condition of a pump power of P≈1.2 P crit (0) for a Gaussian-shaped case with a duration of 20 fs in argon (Ar). Whereas under the conditions of our case in which the pump field consisting of a LG ω 0 and a Gaussian 2ω 0 , the clamping intensity would be bigger than 100 TW/cm 2 . In our numerical simulation, we choose a moderate pump level (10 TW/cm 2 <I pump <100 TW/cm 2 ). Note that N 2 has similar ionization potential as Ar, and the filament length in our simulation is less than 4.5 cm. Basically, it can be almost equivalent to the early stage of a Gaussian-shaped case with a duration of 20 fs aforementioned, in which Kerr nonlinearity dominates over plasma current [8]. Then the FWM model is responsible for the ultra-broadband THz radiation generation. Given that complete 3 + 1 dimensional numerical simulation is very time consuming, we will also take the strategy in Ref. [11], namely, dividing the propagation into linear and nonlinear propagation. Please refer to Appendix for more details.

Numerical Simulation and Discussion
After the linear propagation, the transverse electric field distribution of the two-color field at t=0.5T (T is the time cycle) is shown in Fig. 1(a). It seems like an asymmetric Taichi-like profile due to the asymmetry of the two-color field, the focusing effect via the lens and the linear propagation. The time-integrated intensity distribution keeps a doughnut pattern [ Fig. 1(b)]. However, the center is with non-zero intensity due to Gaussian 2ω 0 , which can be seen from cross section intensity profile curve, as shown with the white line in Fig. 1(b). It is worth noting that a small temporal offset will accumulate with distance. Over the linear propagation zone, we find that, the ω 0 −2ω 0 time delay is about +20.01 fs in temporal domain, where the sign '+' means the ω 0 is ahead of the 2ω 0 .
Generally, based on the FWM model, there are two processes taking place simultaneously to generate THz emission in the two-color gas-plasma filaments, namely, ω 0 +ω 0 −2ω 0 →ω THz and -(ω 0 +ω 0 ) + 2ω 0 →ω THz [19,20]. It has been revealed that the generation efficiencies of the two FWM processes are extremely sensitive to the temporal walk-off between the ω 0 and 2ω 0 pulses, namely, the ω 0 −2ω 0 time delay [19,21]. For a short plasma filament (i.e., millimeter-scale), Kerr nonlinear effect is relatively weak, so the influence of the ω 0 −2ω 0 time delay on the two FWM processes can be very weak, too. When the filament length increases (i.e., centimeter-scale), the ω 0 −2ω 0 time delay starts to take effect. As a result, one of the FWM processes is growing while the other is diminishing as the ω 0 −2ω 0 time delay increases within a period of ∼0.67 fs [21,22]. In the case of this shaped two-color field, the OAM conservation laws of the two FWM processes of ω 0 +ω 0 −2ω 0 →ω THz and -(ω 0 +ω 0 ) + 2ω 0 →ω THz can be expressed by 1 -h+1h−0 -h=2h and -(1 -h+1h) + 0 -h=−2h, respectively, in which two THz vortex pulses with TC of +2 and −2 can be produced. As is known, an interference vortex beam [23] was produced by the superposition of two collinear vortex beams with variable relative amplitudes and conjugated TCs. Accordingly, when they interfere with each other, different THz interference vortex beams are predictable to be generated [23].
In the low-intensity linear propagation zone, the Kerr self-focusing and plasma generation still can occur slightly. A reference frame is established where the laser pulse propagates along the Z axis direction and the focus of the lens (f =0.3 m) is selected as the origin of coordinates. For simplicity, we assumed that the plasma channels start at z=−0.05 m, namely, at the end of the linear propagation zone, and then to evaluate numerically Eq. (1) for the nonlinear propagation characteristics in the gas-plasma filaments. Figure 1(c) shows the calculated spectra of the output electric field at the filament lengths of 7.4 and 28.6 mm. Compared the spectrum at 28.6 mm with that at 7.4 mm, it broadens slightly due to self-phase modulation. The inset shows the ultrabroad bandwidth ranging from 0.1 to 40 THz, peaked at ∼16 THz. In other words, the generated THz spectrum predominates at tens of THz and is depleted near zero frequencies, which are the typical characteristics for THz generation induced by Kerr nonlinearity.
Before discussing the results in our condition, it would be worth recalling the phase variance at different frequencies when nonlinear interactions of Gaussian ω 0 and vortex 2ω 0 pulses with underdense point-like plasma occur in Ref. [5], which is calculated by Fourier-transformed THz source term ∂ t J. To be specific, the larger the frequency, the smoother the azimuthal intensity modulations are, and the more linear the phase variance is. In fact, our partial results follows this rule, too. In the following, we will focus on the spatial mode characteristics of the generated ultra-broadband (0.1-40 THz) THz radiation in our parameter condition. In order to navigate more easily and comfortably, we put the simulation results in the form of a table as Fig. 2. First of all, for a short plasma filament (several millimeters, for example, 7.4 mm or shorter), the FWM scheme indicates that the amplitudes of the generated THz waves change as a cosine function of the relative phase between ω 0 and 2ω 0 , i.e, E THz ∼cos(∆φ)=cos(2θ) in our shaped two-color field, so the generated THz pulse is a ring-shaped beam with an azimuthally-modulated intensity profile and π-stepwise phase profile along azimuthal angle θ, resembling a necklace with four petals [see the first figure at the first row in Fig. 2]. The white circle in the intensity distributions indicates the circle with the radius of the maximum THz electric field distribution, referred to by r max . We call it THz necklace beams [24,25]. Its π-stepwise phase variance is similar to that of a type of interference vortex beams which is generated by two vortex beams with same amplitudes and conjugated TCs [23], shown in the second figure at the first row in Fig.  2. In specific, the θ-dependent phase keeps constant crossing each lobe but jumps when crossing the boundaries of the lobes, which indicates the polarity reversal in the THz waveforms from one lobe to another. This π-stepwise phase profile with the strongest nonlinearity contribution shows that it carries no OAM [23]. The temporal THz waveforms along azimuthal angle θ are illustrated in the third figure at the first row in Fig. 2. When θ=π/4, 3π/4, 5π/4 and 7π/4, the amplitudes of the THz wave almost vanish, which is consistent with FWM scheme [20,24]. Clearly, the four representative waveforms at different θ [marked as vertical white dash lines and donated by A, B, C, and D in the third figure at the first row in Fig. 2] are plotted in the fourth figure at the first row in Fig. 2. The amplitudes of the THz pulses at the azimuthal angles A and C almost keep at the same value with an opposite polarities of the THz fields, so do the waveforms at the positions B and D. However, the peak amplitude at B is much smaller than that at A. In the light of the aforementioned analysis, the THz necklace beams with the strongest phase nonlinearity are obtained.
As the filament length increases slowly, we found that the THz phase tends to a nonlinear function of θ and the nonlinearity of the phase variation becomes weaker compared to the case of the millimeter-scale filament. It is similar to another general type of interference vortex beams which is generated by two vortex beams with different relative amplitudes and conjugated TCs [23]. It is referred to as a THz AAVB [7]. We take the filament length of 10.8 mm for an example. The THz spatial intensity profile still keeps a necklace pattern [see the first figure at the second row in Fig. 2], while its phase increases nonlinearly and smoothly by about 4π when θ varies from 0 to 2π [see the second figure at the second row in Fig. 2]. Here the evolution of the THz waveforms vs θ becomes a little more complicated. The cos(∆φ)-dependency of THz amplitude will not be followed any more, i.e., its minimum does not go to zero [see the third figure at the second row in Fig. 2]. Correspondingly, the four representative waveforms are plotted in the fourth figure at the second row in Fig. 2. The waveforms of A and C, or B and D have the same changing rule as the case at 7.4 mm, but the amplitudes of B and D are higher.
Interestingly, upon further propagation, the nonlinearity of THz phase variation becomes weakest, i.e., its phase changes almost linearly. Correspondingly, it is referred that the THz vortex beams are generated. In specific, this filament length is about 13.1 mm. The intensity distribution [see the first figure at the third row in Fig. 2] shows the center singularity, but its shape is a little different from a LG pulse, which may be caused by the diffraction and dispersion in the filaments [26]. Compared the phase curve at 13.1 mm with that at 10.8 mm, the THz phase increases linearly with a slope of +2 basically, which indicates THz radiation is a vortex pulse with its TC of +2. Obviously, the waveforms typically vary continuously along the θ [see the second figure at the third row in Fig. 2]. Intuitively, their positive or negative electric field peaks with almost same peak amplitude in the image exhibit downward inclined stripes [see the third figure at the third row in Fig. 2]. Clearly, the four typical waveforms plotted in the fourth figure at the third row in Fig. 2 are with different phases.
After that, we found that the THz phase nonlinearity begins to increase and reaches to a relative stronger one, and then decreases again and almost reaches to the weakest, which indicates that THz AAVBs and THz vortex beams appear in turn almost periodically. We take three cases (i.e., ∼15.4, ∼19.3, and ∼28.6 mm) as examples to observe the evolution rule of the generated THz radiation. Their results are shown at the fourth to sixth row in Fig. 2, respectively. When the filament length is ∼15.4 mm, an intensity profile along θ with the intensity modulation depth of smaller than 1 is presented [see the first figure at the fourth row in Fig. 2]. The phase along θ also shows nonlinear variance [see the second figure at the fourth row in Fig. 2], similar to the case of 10.8 mm [see the second figure at the second row in Fig. 2], but with weaker nonlinearity. In a similar vein, the positive or negative field peaks in the image exhibit almost nearly-inclined stripes [see the third figure at the fourth row in Fig. 2], a little similar to the case of 13.1 mm. However, color depth variance is shown along the nearly-inclined stripes, in which the peak amplitudes are not equal. Particularly, the ratio of the peak amplitudes at A and B is about 2.1, as shown in the fourth figure at the fourth row in Fig. 2. Furthermore, when the filament length reaches to ∼19.3 mm, its intensity distribution turns back to be necklace-shaped again [see the first figure at the fifth row in Fig. 2], and the phase nonlinearity becomes stronger compared with the case of 15.4 mm, as shown in the second figure at the fifth row in Fig. 2. The THz waveforms along the whole azimuthal range [see the third and fourth figures at the fifth row in Fig. 2] are with opposite polarities as those in the case of 10.8 mm. In addition, it is possible to get the THz vortex beams again when the filament length reaches to ∼28.6 mm [see the sixth row in Fig. 2]. However, compared to the case of 13.1 mm [see the third row in Fig. 2], diffraction, dispersion effects, spatiotemporal instabilities and etc. would be more obvious. As a result, the intensity distribution of the generated THz beam here is a little different from that in the case of 10.8 mm, as shown in the first figure at the sixth row in Fig. 2. However, it can be still regarded as a THz vortex pulse due to the nearly-linear phase variance along θ [see the second figure at the sixth row in Fig. 2] and the waveforms with nearly-inclined stripes [see the third figure at the sixth row in Fig. 2]. By and large, it seems that, for a filament with the millimeter-scale length, the generated THz radiation performs with the THz necklace beams; with the increase of the filament length, THz AAVBs and THz vortex beams appear in turn almost periodically.
Obviously, necklace beams, AAVBs, and even vortex beams are three different types of interference vortex beams [23]. Actually vortex beams can be considered as a special case of AAVBs, and both of them carry OAM. Accordingly, we may reasonably assume that here the two processes of ω 0 +ω 0 −2ω 0 →ω THz and -(ω 0 +ω 0 ) + 2ω 0 →ω THz do occur but with variable probability amplitudes in the process of THz generation. Based on the OAM conservation law, the former process will generate THz emission with a TC value of +2, and that of the later will be with −2. Consequently, when these two THz vortex beams with same amplitude and conjugated TCs of +2 and −2 interfere with each other, the THz necklace beams with no OAM are produced, which is consistent with the case of the millimeter-scale filament length [see the first row in Fig.  2]. In contrast, when they are with different amplitudes, the THz AAVBs with nonlinear phase profile occur. This is in accordance with the most cases like ∼10.8 mm [see the second row in Fig. 2], ∼15.4 mm [see the fourth row in Fig. 2] and so on. When one process is much more dominated over the other, or one amplitude is much larger than the other, then the THz vortex beams are acquired. This is agreed with the cases of ∼13.1 mm [see the third row in Fig. 2] and ∼28.6 mm [see the sixth row in Fig. 2].
According to the qualitative discussion above, the generated THz radiation can be regarded as the superposition of the two THz components with conjugated TC values (i.e., +2 and −2). Next, we begin to give a quantitative analysis by calculating the mode purity of the generated THz beams via expanding the target electric field into the superposition of a series of LG p,l modes [27]. The purity a p,l of LG p,l in the target electric field E T for a coherent superposition can be expressed as a p,l =(<E T |LG p,l ><LG p,l |E T >)/(<E T |E T ><LG p,l |LG p,l >), where the inner product < B|A > of the two electric fields A and B is defined as ∫∫ B * (x, y)A(x, y)dxdy. The calculation results of the mode purity at different filament lengths are shown in Fig. 3. The insets are the corresponding intensity images in the CCD plane [28,29]. The results at different filament lengths indicate that the generated THz emission does have two TC components, i.e., +2 and −2. For the filament with a length of ∼7.4 mm, the purities for TC of +2 and −2 are almost the same [ Fig. 3(a)], which indicates two FWM processes with almost same possibility amplitude. As a result, the interference between the two FWM processes is strongest, which results in the generation of the THz necklace beams with four intensity maxima but without OAM [see the first row in Fig. 2]. As the filament length increases gradually, it is found that the purity for the TC of +2 increases, while that of −2 decreases [compare Figs. 3(a), 3(b) and 3(c)]. Correspondingly, the THz AAVBs with nonlinear phase profiles are generated by the coherent superposition of the two THz vortex beams with different amplitudes [ Fig. 3(b)]. The purity for TC of +2 can even reach to more than 0.9, while the other one decreases to less than 0.1, as shown in Fig. 3(c), namely, the latter can be negligible. Here the THz vortex beams with TC of +2 are obtained. Subsequently, the purity for TC of +2 begins to decline, and then grow again [ Figs. 3(d)-3(f)]. In other words, the variance trend of the THz phase nonlinearity experiences decreases, increases and then decreases again, which indicates that after the THz necklace beams, the THz AAVBs and the THz vortex beams appear in turn almost periodically. According to the topological competition [30], the "winning" TC of an interference vortex beam is the one with the larger absolute value of the amplitude. Accordingly, the TCs of the THz AAVBs and the THz vortex beams in these cases are all equal to +2 [ Figs. 3(b)-3(f)]. However, here comes to two questions: First, why does the possibility amplitude of one FWM process (i.e., ω 0 +ω 0 −2ω 0 →ω THz ) grow and decline almost periodically as the filament length increases? Second, why is this FWM process (i.e., ω 0 +ω 0 −2ω 0 →ω THz ) always dominated over the other one with the increase of the filament length? For the first question, physically, we found that it is the combined effect of the ω 0 −2ω 0 time delay [ Fig. 4(a)] and the separate, periodical and helical plasma channels [ Fig. 4(b)] that plays a vital role in the generation efficiencies of the two processes. In detail, on the one hand, it has been demonstrated that [19,21], for a general two-color laser field, the generation efficiencies of the two FWM processes are very sensitive to the ω 0 −2ω 0 time delay. The intensity changes with the time period of ∼0.67 fs, corresponding to a distance of ∼200 nm [19]. Based on this, we checked the time delays of the ω 0 and 2ω 0 pulses at different positions. It is found that the ω 0 −2ω 0 time delays at filament lengths of 0, 7. The relationship between ∆t r and the filament lengths is shown in Fig. 4(a). It seems that, the closer the time delay is to an integer multiple of the period, the weaker the interference is, which means only one FWM process is more dominant, namely, ω 0 +ω 0 −2ω 0 →ω THz in this case. For the cases at 13.1 and 28.6 mm, the interference almost vanishes. For the other cases whose time delays are relative far away from an integer multiple of the period, their phase distributions show nonlinear variances, which indicates that the other process of -(ω 0 +ω 0 ) +2ω 0 →ω THz plays a role. The THz spatial modes depend strongly on the time delay, which, of course, also affects the yield of the THz radiation [21,22,31]. On the other hand, generally, during gas-filamentation process, which always considers Kerr self-focusing and plasma effects, the vortex beams are unstable to the transverse perturbations and would break up randomly [32,33]. This is referred to as modulation instability [34]. However, in the case of the shaped two-color field, according to the Eq. (10) of plasma density in the Appendix, it seems that the plasma density is proportional to |cos(∆φ)|= |cos(2θ)|, which leads to four separate plasma spots in transverse plane. That's to say, whether they are point-like sources (without considering Kerr nonlinearity) or long plasma filaments (considering Kerr and plasma effects), there always are four separate plasma spots at a certain position in the transverse plane. Furthermore, as the filaments increase, the dispersion causes the change of the phase difference between ω 0 and 2ω 0 . Thus, the four separate plasma spots would rotate, thereby leading to the generation of helical plasma channels. As a result, the shaped two-color field results filamentation with four separate and helical plasma channels in a regular and predictable manner.
Moreover, for clarity, it is necessary to present global curves of the peak intensity and the corresponding electron density vs filament length. As shown in Fig. 4(c), one can see that, at the beginning of the filaments, the peak intensity of pump field is about 31.25 TW/cm 2 , which corresponds to the power of P=6.26 GW, satisfying P≈1.2 P crit (1). The peak intensity decreases slightly. Correspondingly, the electron density decreases [see the inset in Fig. 4(c)]. Subsequently, due to self-focusing effect, peak intensity increases intensively to about 50.86 TW/cm 2 at ∼20.2 mm. During this time, more ionizing electrons are generated at first. However, with the growth of the electron density, defocusing effect and ionization loss also increase, which could cause the slight decline of peak intensity. The electron density decreases again. Later, it seems that the electron density varies repeatedly with higher values. Specifically, there are actually three bursts of plasma, i.e., 1.16×10 13 cm −3 at ∼0.3 mm, 1.99×10 14 cm −3 at ∼13.5 mm, and 1.45×10 15 cm −3 at ∼25.5 mm, respectively, whose levels are basically consistent with those in the Refs. [35,36], i.e., when the peak intensity is about 30 TW/cm 2 , a weakly ionized plasma with the electron density of around 10 13 cm −3 is generated typically [35]. Actually, this periodic variation can be the result of a cyclic process of self-focusing and plasma defocusing via dynamic spatial replenishment [37].
In addition, we compared the temporal waveforms and transverse patterns of the full electric field at four different filament lengths (i.e., ∼7.4, ∼13.1, ∼19.3 and ∼28.6 mm), shown in Fig.  4(d). As the filament length increases, the waveforms of full electric field vary, especially the envelopes experience deformation due to THz and Brunel harmonics generation [38] [see the first row in Fig. 4(d)]. And the intensity distribution at 0.5 T becomes more and more asymmetry, at the same time, it shows the intensity increase and maxima rotation [see the second row in Fig.  4(d)]. However, the time-integrated intensity distributions almost keeps the same because the generated THz Brunel harmonics can be negligible compared with pump ω 0 and 2ω 0 [see the third row in Fig. 4(d)]. Also, we calculated the separated parts that correspond to the ω 0 and 2ω 0 pulses [not show here]. We found the temporal waveforms and time-integrated transverse intensity patterns of the separated ω 0 and 2ω 0 pulses are almost the same as those at the beginning of nonlinear propagation, except for the change of ω 0 −2ω 0 time delay in time domain.
It is found that, when the plasma density is almost at peak (i.e., the filament length is around 13.1 mm), and meanwhile, the ω 0 −2ω 0 time delay is much closer to integer cycles of the time delay period (i.e., ∼0.67 fs), almost only one FWM process occurs, i.e., ω 0 +ω 0 −2ω 0 →ω THz in this case. Accordingly, the THz vortex beams are generated. For the second question, we found that the sign of the ω 0 −2ω 0 time delay determines the dominated process, thereby determining the TC sign of the generated THz AAVBs and THz vortex beams. When it is positive (i.e., ω 0 is ahead of 2ω 0 ), the FWM process of ω 0 +ω 0 −2ω 0 →ω THz is dominant. The aforementioned cases are all based on this. As expected, it turns out that, when the ω 0 −2ω 0 time delay is adjusted to be negative (i.e., 2ω 0 is ahead of ω 0 ), which can be realized by inserting a phase compensator in the practical experiments [21,22], the other FWM process of -(ω 0 +ω 0 ) +2ω 0 →ω THz becomes dominated. The corresponding numerical results at the different filament lengths are shown in Fig. 5. Basically, the evolution is much similar to the cases discussed above. To avoid repetition, the detailed discussion will not be covered here.

Conclusion
In conclusion, we have demonstrated an effective way to generate three different THz singular beams with ultrabroad bandwidth ranging from 0.1 to 40 THz via long gas-plasma filaments excited with a shaped two-color laser field. It is revealed that the three different THz singular beams can be obtained at different filament lengths, including the THz necklace beams, the THz AAVBs and the THz vortex beams. The THz necklace beams occur at millimeter-scale length; after that, the THz AAVBs and the THz vortex beams appear in turn almost periodically. It is found that the two FWM processes with different probability amplitudes interfere with each other at different filament lengths and would be responsible for the different THz spatial modes. Physically, the combined effect of the ω 0 −2ω 0 time delay and the separate, periodical and helical plasma channels is critical for the generation efficiencies of the two processes. Besides, the sign of the ω 0 −2ω 0 time delay determines which process is dominant, thereby determining the TC sign of the generated THz singular beams. We believe that this simulation-based study will prompts the development for exploitation of the special THz singular light beams and investigation on nonlinear dynamics in the process of filamentation.

Appendix: Laser Parameters and Numerical Approach of Unidirectional Propagation Equation
In the beginning in our simulation, the space-time electric field E 0 ≡E 0 (x,y,t;z=-f ), including the original electrical fields (E ω0 and E 2ω0 ) of ω 0 and 2ω 0 , before the lens can be described as It should be noted that the LG ω 0 field is expressed by superposition two Hermite-Gaussian beams. Then the electric field in space-frequency domainÊ 0 ≡Ê 0 (x, y, ω; z = −f ) at the plane immediately behind the lens shall bê here i is imaginary unit, r=(x 2 +y 2 ) 1/2 and θ=arg(x,y) are the coordinates in polar notation. A 01 and A 02 are the field amplitudes, τ 1 and τ 2 are the pulse durations, while a 1 and a 2 are the spot sizes of ω 0 and 2ω 0 , respectively. We set a 2 =2a 1 in order to increase their transverse overlapping area. The wave number is k(ω)=n(ω)ω/c, where n(ω) is the frequency (ω)-dependent refractive index, and c is the light speed in vacuum.
In the regions near the back plane of the lens and far away from the focus point where the laser intensity is relatively low compared with that in the filaments, we consider it as linear propagation. Consequently, in the linear propagation zone, we will first solve the linear propagation equation by neglecting the nonlinear terms in Eq. (1), i.e., which can be evaluated analytically. In terms of Eqs. (3) and (4), the solution in the position z can be expressed analytically aŝ E(x, y, z, ω) = √ πτ 1 A 01 8a 1 (z)p 2 1 q 2 1 (x + iy)exp with p 1 = 1/a 2 1 (z) + ik(ω)/2f , p 2 = 1/a 2 2 (z) + ik(ω)/2f , Here, we set f =0.3 m, z=−0.05 m, a 1 (z) = 160 µm, a 2 (z) = 2a 1 (z), τ 1 =τ 2 = 30 fs, I 01 =20 TW/cm 2 and I 02 =0.2I 01 , A 01 =(2ηI 01 ) 1/2 , A 02 =(2ηI 02 ) 1/2 , where η=ε 0 n 0 c is the optical impedance in N 2 . It is worth noting that a peak intensity of 20 TW/cm 2 corresponds to the power of 4.02 GW here, which is much higher than the Gaussian critical power (P crit (0) = 1.27 GW for a 800 nm Gaussian pulse), but lower than the LG critical power (P crit (1) = 4 P crit (0) = 5.09 GW for a 800 nm LG pulse with TC of l=1) [18]. Usually, Kerr nonlinearity-induced self-focusing prevails over diffraction when the power of the beam exceeds the corresponding critical power [39]. Thus, in the regions near the back plane of the lens and far away from the focus point, nonlinear effects can be ignored, and a linear numerical solver is employed. After that, the plasma channels are formed, then it will be taken as nonlinear propagation. Alternating direction implicit (ADI) Peaceman-Rachford finite difference scheme [13] is used to solve Eq. (1) numerically.
Nowadays, a number of numerical algorithms have been employed to solve the unidirectional propagation equation. For example, the Crank-Nicolson method [40] is a popular implicit method for solving partial differential equations with second-order accuracy in both time and space [41,42], and is unconditionally stable and strictly non-dissipative [43]. However, the approximate solutions can still contain (decaying) spurious oscillations if the ratio of propagation step to the square of space step is large (typically larger than 1/2) [41,42]. Also, the high cost is seriously obstacle to using the Crank-Nicolson method in practice. Whereas the ADI scheme [44] was show to be accurate and stable to simulate the propagation in comoving frame and paraxial approximation [45]. Particularly, the ADI Peaceman-Rachford scheme shows to be much suitable and high efficient for nonlinear interaction of two co-propagating laser beams with underdense plasmas, which gives additional savings in computer CPU time over the Crank-Nicolson method [11]. That is why we choose ADI Peaceman-Rachford scheme.
The transient plasma dynamics are described by free electron density [46], i.e., ∂N e ∂t = R 1 (N at − N e ) + R 2 N e , (10) where N at =2.7×10 19 cm −3 is the density of neutral atoms, with R 1 and R 2 being the laser field and the avalanche ionization rates, respectively. R 1 can be obtained from Keldysh formula, while R 2 =σ(ω 0 )E 2 /U i . Here, σ(ω 0 )= (e 2 /m e )[v c /(v c 2 +ω 2 )] is the inverse Bremsstrahlung cross-section. U i is the atomic ionization potential of N 2 , while e, m e and v c denote electron charge, electron mass, and electron-ion collision rate, respectively.