Multifrequency Bessel beams with adjustable group velocity and longitudinal acceleration in free space

The group velocity of an optical beam in free space is usually considered to be equal to the speed of light in vacuum. However, it has been recently realized that, by structuring the beam’s angular and temporal spectra, one can achieve well pronounced and controlled subluminal and superluminal propagation. In this work, we consider multifrequency Bessel beams that are known to propagate without divergence and show a variety of possibilities to adjust the group velocity of the beam by means of designed angular dispersion. We present several examples of multifrequency Bessel beams with negative and arbitrary positive group velocities, as well as longitudinally accelerating beams and beams with periodically oscillating local group velocities. The results of these studies can be of interest to scientists working in the fields of optical beam engineering, light amplitude and intensity interferometry, ultrafast optics, and optical tweezers.


Introduction
Structured light with controlled degrees of freedom has wide-range applications in science and technology [1,2]. For example, optical beams exhibiting orbital angular momentum, specific polarization distributions, or propagation invariance are used in optical microscopy, information optics, and laser micromachining [1][2][3]. Recently, the possibility to control the group velocity of structured cw and pulsed optical beams in free space has attracted considerable attention in the optical research community. It has been realized that the group velocity depends not only on the material dispersion of the propagation medium, but also on the angular dispersion, i.e., the relation between the propagation direction and the frequency of the plane-wave components of the field in question. Hence, by tailoring the angular dispersion of an optical beam, the group velocity can be controlled also in free space [4]. There are many demonstrations of structured light fields with specifically adjusted angular dispersion possessing group velocities that are higher (superluminal) or lower (subluminal) than the speed of light in vacuum c, or even negative (see, e.g., [5][6][7][8]).
We refer to the group velocity as the velocity of the intensity peaks along the beam axis [9][10][11]. It is worth noting that the group velocity could also refer to the expected value of the group velocity operator, the weak (local) value of which can be superluminal or negative for inhomogeneous monochromatic superoscillating and vortex beams [12,13]. Superluminal velocities can also be related to the Wigner time delays of light pulses transmitted through a dispersive barrier [14][15][16]. Recently, such time delays corresponding to both subluminal and superluminal propagation were realized by reflecting and refracting optical pulses even without dispersion, through coupling of spatial and temporal degrees of freedom of optical vortex pulses [17]. In general, these group velocity concepts differ from those of the signal and energy-transport velocities that are limited by the speed of light in vacuum [10,[18][19][20]. So does also the velocity of the intensity centroid in accelerating beams that, following the Ehrenfest theorem, has to be subluminal [21].
Bessel beams belong to a distinctive family of self-healing and non-spreading optical beams [22]. Examples of other such beams are elliptical Mathieu beams [23] and parabolic Weber beams [24]. Bessel beams are composed of plane waves with a common angle of propagation with respect to the beam axis. Therefore, they form a convenient basis for constructing multifrequency beams with controllable group velocity [25]. Pulsed multifrequency Bessel beams, which neither diffract nor disperse, such as X-waves [26] and focus wave modes [27,28], have also been studied extensively. So far, pulsed Bessel beams that propagate at subluminal [29,30], superluminal [31][32][33][34][35], and accelerating group velocities [36,37] have been demonstrated. Alternatively, the angular dispersion of pulsed Bessel beams can be used to compensate for the material dispersion of the propagation medium, resulting in both diffraction-free and dispersion-free propagation [38]. In addition, Bessel beam pulses (coined as X-wave bullets) with negative group velocity due to extremely tilted pulse fronts of the beam's plane-wave components have been proposed and demonstrated by calculations [39].
In this work, we consider in general terms multifrequency Bessel beams and show how their group velocity can be controlled in free space via the designed angular dispersion. Not only pulsed beams that contain fully correlated frequency components, but also cw beams with uncorrelated components are studied. We show that cw beams in the form of a superposition of a few independent quasimonochromatic Bessel beams can have a well-defined adjustable group velocity. Exhibiting the wave beating phenomenon, they behave essentially deterministically if the beating period is much shorter than the beating harmonicity time that is mainly determined by the spectral width of the interfering components [40,41]. Both pulsed and cw beams can have arbitrary positive and negative as well as spatially and temporally varying group velocities, which is demonstrated by examples. In general, the group velocity depends on the location of the intensity peak, and the local group velocities for a cw multifrequency Bessel beam can be found by solving a certain initial value problem, the general form of which is derived in this work. We show that the local group velocity can not only change upon the beam propagation, but also oscillate periodically both in space and time. We anticipate that our findings will lead to new applications, e.g., in the fields of optical interferometry, ultrafast optics, and optical tweezers.
The paper is organized as follows. In section 2, we give a brief theoretical description of multifrequency Bessel beams. Sections 3 and 4 show how the group velocity can be calculated for pulsed and cw beams, respectively, presenting also some examples of such beams. The conclusions drawn from our results are introduced in section 5.

Multifrequency Bessel beams
Bessel beams are non-diffractive solutions of the Helmholtz equation in cylindrical coordinates [42]. An accurate description of an electromagnetic Bessel beam requires a vectorial model. The vector Helmholtz equation has many non-diffracting Bessel beam solutions [43]. In this work, we consider a Bessel beam obtained by illuminating an axicon at normal incidence with a plane-wave field propagating in the z-direction and polarized along the x-direction. The resulting Bessel beam has the following electric field components [44]: where E 0 , ω, and k 0 are the complex amplitude, angular frequency, and the wave number of the incident plane-wave, respectively, and α is the angle of propagation with respect to the z-axis for the refracted plane-wave components of the beam. In equation (1), the wave and position vectors are given in cylindrical coordinates (ρ, ϕ, z). Multifrequency Bessel beams contain multiple monochromatic frequency components with the field distribution given by equation (1). The frequency components interfere at the beam axis forming an intensity pattern that propagates with a certain longitudinal group velocity. This velocity depends on the wavefront velocities of the frequency components along the beam axis, i.e., the longitudinal phase velocities. For each frequency component, the longitudinal phase velocity depends on the propagation angle α. Therefore, the group velocity of a Bessel beam is determined by the angular dispersion of the beam.
We consider two types of multifrequency Bessel beams: (1) broadband pulsed beams with correlated frequency components and (2) cw beams composed of multiple narrowband beams with coinciding beam axes. For the pulsed beams, the intensity distribution and its evolution are deterministic [45]. On the contrary, the frequency components of the (considered) cw beams do not correlate and the intensity evolution of the beam is not fully predictable. Nonetheless, the intensity beating of the beam can be considered harmonic within a certain time interval called the beating harmonicity time τ bh [40,41], which is defined as the time it takes for the envelope of a normalized intensity correlation function to decrease to the value of 1/2. References [40,41] present an approximate expression for the beating harmonicity time for two beating fields with Gaussian spectra and Gaussian statistics, which can be generalized for N fields as ( Here,Ī l is the average intensity and Δω l is the spectral width of the beam with index l. Equation (2) shows that for sufficiently narrow spectra, the beating harmonicity time is long with respect to the beating period, resulting in an intensity distribution with a well predictable evolution.

Group velocity of pulsed Bessel beams
A pulsed Bessel beam can be generated by transmitting a laser pulse through an axicon. The angular dispersion of the generated pulse is determined by the refraction introduced by the axicon. For a conventional axicon, the refraction at the conical surface is shown in figure 1. The refraction angle α is given by Snell's law as where γ is the apex angle and n is the refractive index of the axicon. Due to the material dispersion of the axicon, α is a function of frequency, which leads to an angular dispersion of the transmitted pulse found by differentiating the refraction angle with respect to frequency ( This equation shows that the angular dispersion is linearly proportional to the material dispersion of the axicon (∂ ω α ∝ ∂ ω n). Alternatively, Bessel beams can be generated by using diffractive (holographic) axicons [46,47], where the diffraction introduces anomalous angular dispersion (i.e., ∂ ω α < 0), or by using an annular slit and the Fourier transform property of a lens [48], in which case the angular dispersion is introduced by chromatic aberration of the lens. The longitudinal group velocity of an angularly dispersed pulse is given by the projection of the group velocity vector onto the beam axis, which in free space yields The group velocity of a pulse as a function of angular dispersion ∂ ω α is shown in figure 2. A pulse with sufficiently high anomalous dispersion ∂ ω α < (cos α − 1)/(ω sin α) has subluminal group velocity, meaning that v g < c. In the range (cos α − 1)/(ω sin α) < ∂ ω α < cot α/ω, the group velocity is superluminal (v g > c). For dispersion ∂ ω α = cot α/ω, the group velocity function has a pole. If the angular dispersion of the pulse is greater than that at the pole, the group velocity is negative, meaning that the  intensity peaks move towards the light source. A similar relation has been found in [6] between the group velocity of a two-dimensional optical space-time packet and the tilt angle of a spectral hyperplane intersecting the light cone.
In the literature on pulsed Bessel beams, angular dispersion is often neglected and the beams are considered to propagate at a group velocity equal to the longitudinal phase velocity c/cos α that is close to c. Literature discussing such pulsed Bessel beams is vast (see, e.g., [35], and references therein). In general, the angular dispersion can lead to pulsed Bessel beams with pronounced subluminal, negative, and superluminal velocities. Recently, papers discussing subluminal pulses have attracted much attention [29,30,49,50]. On the other hand, pulsed Bessel beams with a negative group velocity, while being proposed more than a decade ago [39], have gone unnoticed.
Negative group velocity, as shown in figure 2, requires high normal dispersion, and therefore, an axicon made of a material with high material dispersion could be used. A good candidate for the (axicon) material would be zinc selenide (ZnSe). Alongside the material dispersion, the refractive index of ZnSe is high at visible wavelengths, which further increases the angular dispersion, as shown by equation (4). Furthermore, ZnSe is a common material for optical elements due to its low absorption in the infrared. Figure 3 shows the group velocity as function of the central wavelength for a pulse transmitted through a ZnSe axicon with an apex angle of 136 • , calculated using the Sellmeier model for the refractive index of ZnSe taken from [51].  For the chosen apex angle, beams with wavelengths shorter than 550 nm experience total internal reflection at the conical surface of the axicon. In addition, the extinction coefficient of ZnSe becomes significant for these short wavelengths (<550 nm) [52]. However, at wavelengths between 550 and 600 nm, the on-axis group velocity is negative, so that the pulse should move towards the source. To verify this result, we consider a pulse with a Gaussian spectrum, for which the central wavelength is 580 nm and full width at half maximum (FWHM) is 5 nm, transmitted by the axicon. The intensity distribution at any instance of time t can be calculated from with the electric field components given by equation (1). Figure 4 shows a time series of the calculated intensity distributions in the xz-plane. The pulse is seen to propagate in the negative z-direction with group velocity v g ≈ −2c, in agreement with figure 3. The negative group velocity originates from the refraction of the pulse by the axicon. After the axicon, the angular spectrum of the pulse is localized at α = 59 • that is the propagation angle of the plane-wave components of the pulse. On the other hand, due to the axicon-induced angular dispersion, the pulse has an X-shaped intensity distribution in any longitudinal plane (see figure 5, where the intensity is plotted in logarithmic scale to make low-intensity regions visible). The blue and green dashed lines in figure 5 mark the X-shaped profile of the pulse, and the blue and green arrows show the corresponding propagation directions of these lines. The tilt angle β of the X is 106.5 • in this case (see the figure). Obviously, the line intersection point in the X-profile of the pulse moves backwards upon the pulse propagation, i.e., in the negative z-direction. Similar X-wave bullets with a negative group velocity have been described in [39]. However, the method to generate them based on the angular dispersion by refraction at a single axicon is remarkably simple compared to that in [39].
Not only are pulsed Bessel beams able to have arbitrary group velocities, but it is also possible for such beams to accelerate or decelerate as they propagate. This peculiar behavior is obtained, if the propagation angles of the plane waves contributing to the formation of an intensity peak on the beam axis change during the propagation. Accelerating/decelerating Bessel beams were investigated numerically in [36], and later realized in [37], by placing a lens behind the axicon forming the Bessel beam. The lens makes the angle α of light refraction depend on the radial distance from the optical axis. In the paraxial approximation, the refraction angle introduced by the lens can be added to the refraction angle introduced by the axicon. This results in a radially dependent propagation direction and, therefore, in a Bessel beam accelerating along z. The propagation angle as a function of the longitudinal coordinate z is given by where f is the focal length of the lens and α is the refraction angle near the center of the lens, the same as α in equation (3). Equation (7) is derived in appendix A using the paraxial approximation. The method works for z considerably smaller than f, and the achievable acceleration is moderate. We suggest that a considerably larger acceleration can be achieved by using a highly dispersive (chromatically aberrative) lens, such as a diffractive lens, because then, not only α, but also the tilt angle β of the X-profile will be a function of z. In the next section, we consider a strongly accelerating dual-frequency Bessel beam with only one of the components being focused.

Local group velocity of beating cw beams
Multiple narrowband laser beams aligned on top of each other can form a beating intensity pattern with intensity peaks that move along the common beam axis. In this section, an equation for the velocity of such intensity peaks is derived. We refer to the velocity as the local group velocity, as the velocity is defined only on the beam axis and can depend on the location of the peak. Consider a group of N monochromatic beams with coinciding beam axes. For simplicity we assume that the beams are paraxial and their polarization is the same. In this case, the complex amplitude of the total field at the z-axis is given by where E l , ω l , k zl (ω l ), and φ l are the amplitude, frequency, longitudinal wave number, and phase of the beam with index l, respectively. The longitudinal wave number k zl (ω l ) depends on frequency ω l . Although we are interested in angular dispersion, the following derivation works equally well for material dispersion or the combination of both. The intensity of the total field is proportional to From equation (9), it can be seen that the intensity pattern is a linear combination of a constant term and terms describing beating pairs of the beams. In order to calculate the velocity of an intensity peak in the pattern, we have to determine its location. An intensity peak is a local maximum of an intensity distribution at a fixed instant of time. Therefore, the position of an intensity peak as a function of time, ζ(t), has to satisfy the equation at each instance of time. Furthermore, to ensure that the function ζ(t) locates a maximum instead of a minimum, we require that In an effort to find the function ζ(t) and the corresponding local group velocity v g (t, ζ(t)) = ∂ t ζ, equation (10) is turned into an initial value problem. An initial value ζ 0 is found by solving equations (10) and (11) at t = 0. The ordinary differential equation (ODE) is obtained by differentiating equation (10) with respect to time and solving for the local group velocity, yielding Equation (12) is a general expression for the local, on-axis, group velocity. Let us consider a special case of two beating frequency components of the beam, for which the local group velocity simplifies to a familiar form [10] v Denoting the group velocity of a pair of beating components with indices n and m as Equation (12) can be rewritten as Equation (15) highlights how the local group velocity is determined by the beating of the pairs of the beam components. Notably, if the beating of a pair propagates at a velocity that deviates from the local group velocity of the beam (v nm g = v g ), its contribution to the local group velocity oscillates as a function of time, and so does the local group velocity. Therefore, the local group velocity is constant if and only if v nm where {v 12 g , v 13 g , . . . , v (N−1)N g } is the set of group velocities of the pairs of the frequency components of the beam. In free space, a constant local group velocity is obtained when all of the frequency components have the same convergence angle α of their plane-wave components. In this case, the resulting local group velocity is v g = c/cos α, just like for a nondispersive pulse.
The propagation of the intensity peaks originating at different locations can be analyzed by expressing the local group velocity as a function of the initial peak position: This expression shows that the initial peak position contributes to the local group velocity by introducing constant terms (k zm − k zn )ζ 0 to the expression. Consequently, at any instance of time, different peaks can have different velocities. However, as the group velocities of the pairs v nm g are constant, the velocity of a peak will be periodically varying in time about a constant value, which is independent of the initial position of the peak. Similarly, the phases φ l affect only the initial peak position and instantaneous velocity.
In order for the peaks to move at different average velocities at different intervals along the beam axis, the group velocities v nm g have to depend on the peak location. This can be realized in practice by making the frequencies of the beams vary with the longitudinal coordinate z. A Bessel beam with longitudinally modulated frequency can be obtained by transmitting a broadband optical beam through a spectral filter with a radially modulated transmission frequency followed by an axicon. In the presence of the modulation, the complex amplitude of the total field is where the central frequency ω l of the beam with the index l is modulated by a function h l (z), which results in a spatially varying phase φ l (z). The local group velocity of the field is derived in appendix B. If the frequency changes slowly with z, so that {∂ z h l , ∂ zz h l , ∂ z φ l , ∂ zz φ l } k zl , the local group velocity is given by where each v nm g (ζ) depends on the peak position in accordance with v nm Alternatively, the group velocities v nm g can be made to depend on the peak position, by making the angle α for one or both of the components to be a function of z. In this case, the group velocities of the beating pairs are given by v nm A way to realize a longitudinally changing propagation angle α by using a lens-axicon doublet is proposed in section 3. The technique results in a modulated longitudinal wave number where α(z) is given by equation (7). The initial value problem for the group velocity can be solved numerically. In this work, we do this using an algorithm described in appendix C.

Examples of cw multifrequency Bessel beams
In this section, we apply the derived equations to several particular examples of cw multifrequency Bessel beams. The results are then verified by calculating the time series of the corresponding intensity distributions in the xz-plane. The intensity distributions are calculated to be proportional to where the electric field components are of the form given by equation (1). In all of the examples, we assume that the beams are narrowband and, for simplicity, set the phases φ l to be zero. The amplitudes of the plane waves incident on the axicon are assumed to be E 0l = (cos α l + 1) −1 V m −1 so that the electric fields polarized along the x-direction are equal on the beam axis, resulting in strong interference.

Two-component Bessel beam with a fixed group velocity
For a two-component Bessel beam, the on-axis group velocity is given by v g = (ω 2 − ω 1 )c ω 2 cos α 2 − ω 1 cos α 1 .
As an example, let us assume that the beam components have wavelengths λ 1 = 610 nm and λ 2 = 600 nm, and the plane-wave propagation angle of the first component, α 1 , is fixed at 10 • . The local group velocity of the beam as a function of angle α 2 is plotted in figure 6. The result shows that, by properly selecting the  plane-wave propagation angle α 2 of the second component, one can obtain subluminal, superluminal, and negative group velocities. Because the beam components in a cw beam are not correlated, the phases φ l will fluctuate independently, affecting the instantaneous value of the group velocity. However, for narrowband components, the relative phase drift will be slow and the group velocity will change insignificantly. Indeed, the deviation Δv g of the group velocity from the value given by equation (24) is determined by the beating period L b along the z-axis and the harmonicity time τ bh as if the fields have the same bandwidth Δω. It can be seen that, if Δω becomes comparable to ω 2 − ω 1 , the deviation Δv g becomes comparable to v g , and the value of the group velocity becomes unpredictable. However, if Δω ω 2 − ω 1 , the fluctuations of v g can be ignored.

Longitudinally accelerating two-component Bessel beam
In this section, we consider another type of a two-component Bessel beam. While the first component is collimated in front of the axicon, the second one is slightly focused. We choose the first component to have the same parameters as in the example above, i.e., λ 1 = 610 nm and α 1 = 10 • . For the second component, the wavelength is λ 2 = 600 nm and the plane-wave propagation angle is where f = 5 cm is the distance from the axicon to the focal plane of the second beam. The calculated peak positions ζ(t) and the local group velocities v g (t) for the two peaks of the beam with initial positions ζ 0 of 7.5 and 10 mm are shown in figure 7. The two peaks are seen to accelerate and approach infinite positive and negative group velocities when ζ(t) approaches a point ζ c ≈ 8.8 mm, where the peaks collide and abruptly reduce their group velocities to 0. The acceleration and the subsequent collision of the two peaks can also be seen in figure 8, in which the time series of the intensity patterns of the two-component beam are shown. Such beams could be used as optical tweezers to trap particles in high intensity regions by a dipole force and move the particles towards a fixed point in space [53].

Frequency modulated two-component Bessel beam
In the following example, we consider a two-component Bessel beam that has a periodically spatially modulated group velocity along the beam axis. The first component of the beam is assumed to be the same as in the two examples above, but the second component is obtained by transmitting a collimated beam through a spectral filter with radially modulated transmission frequency. This results in a longitudinal on-axis frequency modulation of the second beam component. We choose the central frequency of this component to coincide with the frequency of the first component. The angle α 2 is chosen to be 5 • and the frequency modulation to follow equation (18) with h 2 (z) = 1 + 0.001 cos(πz/Λ), where Λ = 2 mm. Let us assume for a moment that all frequency components of the beam are initially in phase. Because the modulation amplitude is small, we can also neglect the phase variation φ 2 (z) of the modulated component. In this case, the frequency modulation divides the z-axis into periodically alternating regions of positive and negative local group velocities. For example, for z ∈ [−1 mm, 1 mm], the local group velocity is positive, whereas for z ∈ [1 mm, 3 mm], it is negative. The peaks originating in the neighboring regions of positive and negative local group velocity propagate towards the borders of the regions, such as the one at z = 1 mm. Figure 9 shows the peak positions ζ(t) and local group velocities v g (t) for two peaks with initial values ζ 0 of 0.91 and 1.07 mm, and figure 10 presents the corresponding time series of the longitudinal intensity patterns. Unlike in the previous example, the peaks do not collide, but decelerate towards a stationary interference pattern, as can be seen in figure 10. Obviously, the process cannot continue indefinitely. Indeed, the frequency components can be either correlated, which results in light pulses, or uncorrelated, which makes the interference pattern randomly change in time and space, preventing any indefinite decrease of the pattern period. Accordingly, in the latter case, the group velocity will on average decrease towards Δv g given by equation (25), eventually becoming unpredictable. The period of the intensity variation along z will also become random with the average value of L b = λ 1 λ 2 /(λ 1 cos α 2 − λ 2 cos α 1 ) = 53 μm, independently of the initial intensity distribution.  Figure 11 shows the location ζ(t) and group velocity v g (t) of the intensity peak with an initial position ζ 0 = 0. The corresponding time series of the intensity pattern of the beam is shown in figure 12. It can be seen that the peaks at different locations exhibit the same oscillating behavior, but shifted in time. If the components are narrowband, this periodic behaviour will be essentially deterministic.

Conclusion
The main conclusion of this work is that a variety of multifrequency pulsed and cw Bessel beams can be created to obtain essentially any fixed or varying group velocity in free space. We have introduced the conditions for pulsed Bessel beams to have subluminal, superluminal, and negative group velocities, as well as methods to fulfill these conditions. We also proposed simple methods to realize backward propagating and longitudinally accelerating optical pulses in free space using a controlled angular dispersion. For cw Bessel beams containing multiple frequency components, we showed that a well-defined group velocity at the beam axis can be found, if the components are sufficiently narrowband. We introduced a general calculation technique, formulated in terms of an initial value problem, to find the local group velocity, both in space and time, and demonstrated that the on-axis group velocity of a two-component Bessel beam can have essentially arbitrary values. Furthermore, if one of the components is focused, the group velocity changes upon propagation resulting in a longitudinally accelerating or decelerating cw beam. We also showed that Bessel-like beams with periodically alternating regions of positive and negative, as well as locally oscillating on-axis local group velocities can be created.
The presented analysis of the group velocity for multifrequency Bessel beams in terms of the beam angular dispersion offers new possibilities in engineering of optical beams with extraordinary properties. Such beams can find applications in science and technology, including the fields of optical interferometry, ultrafast optics, and optical tweezers.
For the complex amplitude given by equation (18) we find At an instance of time t, the function for the location of a peak, ζ(t), satisfies the equation where h l = ∂ z h l | z=ζ(t) and φ l = ∂ z φ l | z=ζ(t) . Equation (B.2) is turned into an ODE by differentiating it with respect to time, i.e., where v g h l = ∂ t h l (ζ), v g h l = ∂ t h l (ζ), v g φ l = ∂ t φ l (ζ), and v g φ l = ∂ t φ l (ζ). Equation (B.3) can be rearranged to solve for the local group velocity:

Appendix C. Algorithm for solving local group velocity initial value problems
The process flow chart of the algorithm used for solving the initial value problem set by equations (10)- (12) or (B.4) is shown in figure C1. For a start, an initial guess for a peak location is made, which is used to find a root of equation (10) with Newton's method. Next, the algorithm checks whether or not the root corresponds to a local maximum. If not, the root is removed from the function using Maehly's procedure, so that the same root cannot be found again. The algorithm keeps finding roots of the function by alternating between Newton's method and Maehly's procedure until a local maximum is found. With the local maximum (i.e., the initial value) found, the initial value problem with ODE given by equation (12) or (B.4) can be solved. Here, the classic fourth-order Runge-Kutta method (RK4) is used. Figure C1. The process flow chart of the algorithm used to solve the initial value problem for the local group velocity.