Optical wave-packet with nearly-programmable group velocities

During the process of Bessel beam generation in free space, spatiotemporal optical wave-packets with tunable group velocities and accelerations can be created by deforming pulse-fronts of injected pulsed beams. So far, only one determined motion form (superluminal or luminal or subluminal for the case of group velocity; and accelerating or uniform-motion or decelerating for the case of acceleration) could be achieved in a single propagation path. Here we show that deformed pulse-fronts with well-designed axisymmetric distributions (unlike conical and spherical pulse-fronts used in previous studies) allow us to obtain nearly-programmable group velocities with several different motion forms in a single propagation path. Our simulation shows that this unusual optical wave-packet can propagate at alternating superluminal and subluminal group velocities along a straight-line trajectory with corresponding instantaneous accelerations that vary periodically between positive (acceleration) and negative (deceleration) values, almost encompassing all motion forms of the group velocity in a single propagation path. Such unusual optical wave-packets with nearly-programmable group velocities may offer new opportunities for optical and physical applications. Unprecedented control of optical beams has allowed demonstration of superluminal or subluminal beams by controlling the phase of the wave-packet. Here, nearly-programmable control of a beam’s group velocity results in alternating acceleration and deceleration regimes in a single propagation path.

R ecent studies of non-diffraction beams include the Bessel, X, Y, Airy, and parabolic waves in linear mediums [1][2][3][4][5][6] , the self-trapping, self-focusing, and nonlinear X wave in nonlinear mediums [7][8][9][10][11] , and multidimensional solitons in complex systems [12][13][14][15] . Based on these studies, three-dimensional (3D) localized spatiotemporal optical wave-packets with long propagation distances and nearly-invariant intensity profiles have been widely demonstrated in both linear [16][17][18] and nonlinear [19][20][21][22] mediums. Applications of such wave-packets range from particlemanipulation to bio-imaging and plasma-physics. Beside the well-known properties of long-distance self-similarity and/or selfhealing of these localized spatiotemporal optical wave-packets, tunable group velocity is another degree of freedom relevant to some novel applications. In nonlinear optics, the group velocity υ g = c/n g (n g is the group refractive index and c is the speed of light) of an optical wave-packet can be well controlled by crafting the wavelength-dependent refractive index n λ 23 . However, some very special materials or systems are required [24][25][26][27][28][29][30] , challenging the application of this principle to linear systems where the controllability of n g is very limited.
The Bessel beam is a famous family of diffraction-free beams resulting from conical superposition, which usually propagates with a constant superluminal group velocity in free space [31][32][33][34][35][36][37][38][39] . By phase-modulating the incident cross-section, it is possible to produce self-accelerating Bessel-like beams having arbitrary curved trajectories [40][41][42][43] . In this case, the combination of a Bessel beam and a pulse can produce a superluminal and/or selfaccelerating optical wave-packet 18,44,45 . For example, directly combining a Bessel beam with an Airy pulse can create a selfaccelerating optical wave-packet in transmission materials 18 ; however, the accelerating value is limited within the short pulse duration range. The X wave, demonstrated in both linear and nonlinear mediums, is another important limited diffraction beam [46][47][48] . In linear optics, broadband superposition of slightly distorted Bessel beams can create Bessel-X spatiotemporal optical wave-packets 49,50 , and in nonlinear optics pulsed beams with nonlinear material responses can generate X-shaped light bullets 51 . However, the controllability of group velocities of these Xshaped optical wave-packets is not high. The Airy beam also is a well-studied diffraction-free beam in which main intensity maxima and lobes tend to propagate [52][53][54][55][56][57] . By combing an Airy beam with a pulse, the resulting self-accelerating optical wave-packets can propagate at superluminal group velocities along parabolic trajectories in free space. In the above methods, group velocityvariable optical wave-packets in free space usually correspond to curved or bended propagation trajectories, allowing for novel applications that include particles guiding/trapping along curved paths and self-bending plasma channels generation [58][59][60] . In some other applications instead, straight-line propagation optical wavepackets with tunable group velocities have irreplaceable advantages.
In linear optics, it has recently been shown that spatiotemporal coupling permits high controllability of the group velocity and/or the acceleration of spatiotemporal optical wave-packets. Saari et al. and Abouraddy et al. invented a new type 2D spatiotemporal optical wave-packet by manipulating the spatial and temporal degrees of freedom jointly, where both diffractive spreading and pulse broadening are eliminated 61,62 . This optical wave-packet can be described by a spectral trajectory resulting from the intersection of the light-cone (k x , k z , ω/c) with a tilt plane (k z , ω/c) in spectral-space, where k x and k z are the transverse and longitudinal wavenumbers, respectively, x and z are the transverse and longitudinal coordinates, respectively, and ω is the angular frequency. By adjusting the tilt angle of the plane (k z , ω/c), the group velocity of the optical wave-packet in free space or in transmission materials can be controlled, including all motion forms, i.e. superluminal, subluminal, accelerating, decelerating, and backward-propagation [63][64][65][66][67][68] . Another spatiotemporal coupling method is to control the group velocity of the intensity peak of a focused ultra-short pulse within the extended Rayleigh length (named as sliding focus or flying focus) by combining temporal chirp and longitudinal chromatism. This method was independently demonstrated by Quéré et al. in theory 69,70 and Froula et al. in experiments 71 . In this method, the longitudinal chromatism separates wavelength-dependent focuses along the propagation axis and the temporal chirp controls the appearance times of these focuses, so that the sliding/flying focus possesses a tunable effective group velocity, also achieving all motion forms (superluminal, subluminal, accelerating, decelerating, and backward-propagation). More recently, we theoretically demonstrated a third spatiotemporal coupling method to generate group velocity and acceleration tunable 3D optical wave-packets in free space 72 . The pulsed beam used for the Bessel beam generation is deformed to have an axisymmetric pulse-front which deviates from its plane phase-front. In the generation of the Gauss−Bessel optical wave-packet, the plane phase-front determines the Bessel beam generation in space, while the deformed pulse-front determines the optical wave-packet propagation in time. Consequently, the group velocity and acceleration of the optical wavepacket can be adjusted by changing the pulse-front deformation, also including all motion forms of superluminal, subluminal, accelerating, decelerating, and backward-propagation group velocities. In all these spatiotemporal coupling methods however, by controlling one degree of freedom it is possible to achieve only one determined motion form (superluminal or luminal or subluminal for the case of group velocity; and accelerating or uniform-motion or decelerating for the case of acceleration) in a single propagation path.
In this article, by combining well-designed complex axisymmetric pulse-front deformations with our recently reported method 72 , we achieve a nearly programmable group velocity with several different motion forms in a single propagation path. The created optical wave-packet can fly with superluminal and subluminal group velocities alternately, and the corresponding instantaneous acceleration varies between positive (accelerating) and negative (decelerating) values. For example, when periodically distributed axisymmetric pulse-front deformations are introduced, the group velocity and acceleration of optical wavepackets display periodical variations during propagation, showing an alternate appearance of superluminal−subliminal group velocities and accelerating−decelerating accelerations. This unusual optical wave-packet presenting different group velocity motion forms in a single propagation path may provide new opportunities for applications.

Results
Setup of the method. The schematic diagram of the method is shown in Fig. 1a. A collimated pulsed beam is reflected by a deformable mirror (DM) (or a free-surface mirror) to shape its pulse-and phase-fronts from a flat surface to a required axisymmetric surface. A transmission spatial light modulator (SLM) is positioned just behind DM to correct the phase-front back to the original flat surface for long-distance propagation while keeping the shaped axisymmetric pulse-front unchanged 73 . A beam splitter (BS) samples the shaped pulsed beam into a parabola telescope, which is used for three purposes: first, to image the shaped pulsed beam into the Bessel beam generation region formed by an axicon for suppressing propagation diffraction; second, to enhance the spatial resolution of the phase-front correction that is limited by the pixel size of SLM by beam reduction (ten times in this article); and third, to increase the instantaneous pulse-front variation across the beam aperture also by beam reduction. Finally, an ideal thin axicon is used to generate a Bessel beam in the conical superposition region.
Using the Simulation model and parameters given in the "Methods" section, the optical fields at different locations are calculated. Because the optical fields are always axisymmetric about the propagation axis, only the 2D distributions in the lateral plane of the r−z plane, where r and z are the transverse and longitudinal coordinates, are given. Figure 1b shows the input pulsed beam, and the carrier frequency is multiplied by 0.1 to avoid too fast oscillations for observation. Figure 1c shows the optical field, when only DM (or a free-surface mirror) is considered, both pulse-and phase-fronts are deformed, and Fig. 1f gives the corresponding surface of DM. Figure 1d shows the optical field, when SLM is considered, which has a deformed pulse-front and an unchanged (plane) phase-front, and Fig. 1g gives the corresponding phase-front correction of SLM. Figure 1e shows the optical field after 5-m free propagation, and some diffraction distortions can be found, showing the necessity of the image relay telescope. In Fig. 1c-e, g, a small range is enlarged for observation. The continuous surface of DM won't change the smooth pulse-and phase-fronts in Fig. 1c, while the segmented phase-front-correction of SLM (limited by the SLM pixel size shown in Fig. 1g) generates fine structure modulations in Fig. 1d. The pixel−pixel gap of SLM introduces net-like spatial amplitude modulation, and the pixel size of SLM keeps a very small phasefront tilt within each segment spatial range. In real cases, these two high spatial-frequency modulations can be reduced by suitable propagation diffraction in Fig. 1e. All in all, SLM introduces optimized (−π, π) phase-corrections across the beam Fig. 1 Setup for variable-group velocity optical wave-packet. a In the setup, input Gaussian pulsed beam experiences reflection by a deformable mirror (DM) (or free-surface mirror), phase-front correction by a spatial light modulator (SLM), reflection by a beam splitter (BS), image relay and beam reduction by a parabola telescope, and finally Gauss−Bessel optical wave-packet (OWP) generation by an ideal thin axicon. b Optical field of the input pulsed beam has plane phase-and pulse-fronts. c Optical field with DM (without SLM) has deformed phase-and pulse-fronts. d Optical field with SLM has deformed pulse-front but near-plane phase-front. e Optical field after 5-m free propagation shows slight diffraction distortions. Insets in (c-e) illustrate fine structures in a small regions of 200 μm × 60 fs. Inset in (c) shows smooth phase-and pulse-fronts, inset in (d) shows spatial amplitude modulation and residual small phase-front tilt induced by SLM pixel-pixel gap and pixel size, respectively, and inset in (e) shows smearing effect by propagation diffraction. f Surface deformation of DM. g Phase-front correction of SLM, and a small range is enlarged for observing the pixel effect. For observation, the carrier frequency in (b-e) is multiplied by 0.1 to avoid too fast oscillations. aperture (like a Fresnel lens) and then collimates propagation directions (or wave vectors) at different transverse positions. In this case, after the ideal thin axicon, a Gauss−Bessel optical wavepacket can be created in the superposition region. Due to the pulse-front deformation, the created optical wave-packet does not possess a constant superluminal group velocity governed by υ b = c/cos α anymore 39 , and the detail is going to be introduced in the following section.
Variable-group velocity and acceleration. Figure 2 shows the simulated intensity distributions (illustrated by red−yellow distributions) in 2D x−z/t and 3D x-y−z/t spatiotemporal domains. Because of the axisymmetric profiles about the propagation axis, here we only discuss the results in the 2D lateral plane x−z/t containing the propagation axis. The axicon spatially divides the input pulsed beam into two and changes their traveling directions (illustrated by green dash-line arrows) with (half) conical angles of α = ±0.5°. The plane phase-fronts (illustrated by white solid lines) of two pulsed beams are always perpendicular to the traveling directions. In the superposition region, a Gauss−Bessel spatiotemporal optical wave-packet is generated at the core position resulting from conical superposition. We keep the same definitions as that in our previous article 72 : in space, the propagation axis of the optical wave-packet (or Bessel beam) is defined as the z-axis and the geometrical center of the superposition region is defined as the spatial origin of (x = y = z = 0) (see Fig. 1a); in time, the moment when the intersection of two phasefronts (the intersection of two white solid lines) arriving at z = 0 is defined as the temporal origin of t = 0.
When the deformed pulse-front has a cosine-function-like axisymmetric distribution, Fig. 2a shows at the very beginning the generated optical wave-packet (the core of red−yellow distributions) has the same location with the intersection of the phasefronts (the intersection of two white solid lines); however, during Propagation of variable-group velocity optical wave-packet. a After the ideal thin axicon, in 2D lateral plane the divided two half beams have symmetrical propagation directions (green dash-line arrows) with an angle of 2α. When deformed pulse-front (red−yellow distributions) has an axisymmetric cosine-function-like distribution, longitudinal gap Δz between the intersection of phase-fronts (white solid lines) and optical wave-packet (core of red-yellow distributions) increases and decreases in the firstand second-half propagations, showing subluminal and superluminal group velocities, respectively. b Isosurface plots (10 and 80% of the maximum intensity) show dynamics of optical wave-packet in the entire pulsed beam at different propagating times of −100 ps, −50 ps and 0 in (a). c When deformed pulse-front has an opposite axisymmetric cosine-function-like distribution, longitudinal gap Δz between the intersection of phase-fronts and optical wave-packet decreases and increases in the firstand second-half propagations, showing superluminal and subluminal group velocities, respectively. The (half) conical angle generated by the axicon is α = 0.5°, and then the intersection of phase-fronts travels at a constant velocity of 1.00004c. the propagation, the optical wave-packet will be temporally delayed along the z-axis. From left to right in Fig. 2a, the detailed distributions of the pulse-fronts (red−yellow distributions), the phase-fronts (white solid lines), and the optical wave-packet (the core of red−yellow distributions) at different propagating times of t = −150, −100, −50, 0, 50, 100 and 150 ps are illustrated. The intersection of the phase-fronts travels at a constant velocity of 1.00004c governed by υ g = c/cos α, while the optical wave-packet has a variable-group velocity. In the first half propagation from t = −150 ps to t = 0, the longitudinal gap Δz between the intersection of the phase-fronts and the optical wave-packet increases from zero to the maximum with the propagating time t, while in the second-half propagation from t = 0 to t = 150 ps, it gradually decreases from the maximum to zero again. Figure 3a shows the variation of the longitudinal gap Δz with the propagating time t (or position z), which has a cosine-functionlike distribution about the propagating time t (or position z). The blue curve in Fig. 3b shows the instantaneously variable subluminal group velocities in the first half propagation, reaching the minimum at around t = −100 ps, and the red curve illustrates the instantaneously variable superluminal group velocities in the second-half propagation, reaching the maximum at around t = 100 ps. Figure 3c shows the variation of the instantaneous acceleration during the entire propagation, and the optical wavepacket experiences decelerating, then accelerating and finally decelerating motions in three temporal periods of from the appearance time to around t = −100 ps, from around t = −100 ps to around t = 100 ps, and from around t = 100 ps to the disappearance time, respectively. Figure 2b illustrates three isosurface plots at 10% and 80% of the maximum intensity showing the dynamics at three different propagating times of t = −100 ps, −50 ps, and 0. It is clear that at different propagating times the relative temporal (or z-axis) locations of the optical wave-packet in the entire pulsed beam are different. At the moment of the perfect overlap (i.e. t = 0), a Bessel beam with a central core and a series of concentric rings appears; however, while comparing with the traditional Bessel beam, the central core and the concentric rings have different temporal (or z-axis) locations, which is dominated by the pulse-front deformation.
Next, when the deformed pulse-front has an opposite cosinefunction-like axisymmetric distribution (or with a π phase-shift), Fig. 2c shows the simulated distributions of the pulse-fronts, the phase-fronts, and the optical wave-packet at different propagating times of t = −150, −100, −50, 0, 50, 100 and 150 ps. At the very beginning, the optical wave-packet has the maximum temporal delay along the z-axis relative to the intersection of the phasefronts, i.e., the maximum longitudinal gap Δz. In the first half propagation from t = −150 ps to t = 0, the longitudinal gap Δz decreases to zero with the propagating time t, and in the secondhalf propagation from t = 0 to t = 150 ps, it increases back to the maximum again. Figure 3d shows the variation of the longitudinal gap Δz with the propagating time t (or position z). The red curve in Fig. 3e shows the instantaneously variable superluminal group velocities in the first half propagation, reaching the maximum at around t = −100 ps, and the blue curve shows the instantaneously variable subluminal group velocities in the second-half propagation, reaching the minimum at around t = 100 ps. Figure 3f shows the variation of the instantaneous acceleration, and the optical wave-packet experiences accelerating, then decelerating and finally accelerating motions in three temporal periods of from the appearance time to around t = −100 ps, from around t = −100 ps to around t = 100 ps, and from around t = 100 ps to the disappearance time, respectively.
By comparing Figs. 2a with 2c or 3a-c with 3d-f, the motions of the optical wave-packet in two cases are opposite, which is determined by the two opposite pulse-front deformations. Figure 3a, d shows the variation value of the longitudinal gap Δz (between the intersection of the phase-fronts and the optical wave-packet) is dominated by that of the pulse-front deformation. Because the intersection of the phase-fronts has a constant velocity of c/cos α (for an ideal thin axicon), the instantaneous group velocity of the optical wave-packet can be calculated by the derivative operation of υ b = c/cos α − d(Δz)/dt, and the instantaneous acceleration of the optical wave-packet can be obtained by a = dυ b /dt. In this case, by controlling the value of the pulse-front deformation, the instantaneous longitudinal gap Δz, and accordingly the instantaneous group velocity and acceleration of the optical wave-packet, in theory, can be well controlled. Furthermore, by introducing some unusual shapes of the pulse-front deformations, optical wave-packets with unusual motion forms (e.g., periodically variable-group velocity and-acceleration here) can also be created.
Moreover, in Fig. 2a, c, the best image relay position by the parabola telescope is at z = 0 (or t = 0), and when the pulsed beam propagates at z = ±45 mm (or t = ±150 ps) the propagation-diffraction-induced slight distortions can be found. In this case, when reducing the conical angle α to increase the propagation distance, the propagation diffraction should be considered, and the details are going to be discussed in the "Discussion" section.
Controllability of variable-group velocity. The above result shows that the shape and the value of the pulse-front deformation dominates the motion form and the group velocity value of the optical wave-packet, respectively. In our previous article 72 with the simplest form of pulse-front tilt (PFT), i.e., a linearly tilt pulse-front, the optical wave-packet has a constant group velocity governed by where α is the (half) conical angle formed by the ideal thin axicon, and β is the tilt angle between pulse-and phase-fronts. For a complex pulse-front deformation, for example the cosine-function-like profile here, Fig. 4a shows the instantaneous PFT can be obtained by calculating the tangent angle β of the tangent line at the intersection of two pulse-fronts. Figure 4b shows the dependence of the optical wave-packet group velocity υ b on PFT for three different (half) conical angles α. The tunable (or variation) range of the group velocity υ b increases with increasing PFT, and the capability can be further enhanced by choosing a larger (half) conical angle α. If PFT is enlarged, the difficulty of generating the required phase-correction by SLM will increase. Figure 4c shows when the pulse-front deformations have the same cosine-function-like shape but different peak-valley (PV) values, the required phase-corrections by SLM are quite different. When the PV value is enlarged, the instantaneous PFT changes drastically along the transverse axis, accordingly the variable range of the group velocity υ b increases; however, a high spatial-frequency phase-correction is also required, challenging the SLM resolution. Figure 4d shows the number of the phasecorrection periods (−π, π) within 1 mm increases linearly with increasing PFT. If one modulation period (−π, π) contains at least five pixels, the required SLM pixel size has a negative exponential function distribution with respect to PFT. For example, with reference to an available 4 μm pixel size, when the absolute value of PFT is larger than 280 fs mm −1 , beam reduction by the image relay telescope becomes necessary. This is another reason why a parabola telescope with ten times beam reduction is used in Fig. 1a. And the third reason is to reduce the beam size to increase the instantaneous PFT across the beam aperture and, accordingly, increase the variation range of the group velocity υ b .
The variation range of the (propagating time/position dependent) group velocity υ b is dominated by that of PFT across the deformed pulse-front. For the case of a cosine-function-like pulse-front deformation used in this article, the pulse-front is z = L/2·cos(r/D·2π), where L and D are the longitudinal PV value and the transverse period, respectively. The tangent angle can be obtained β = −L/D·π·sin(r/D·2π), and the extremum is β max/min = ±L/D·π. Based on an available commercial DM with a 100 μm PV value and a 20 mm diameter 74 , i.e., L = 100 μm and D = 10 mm for an axisymmetric cosine-function-like pulse-front deformation, the extremum is β max/min ≈ ±31.41 mrad (or PFT ≈ ±104.7 fs mm −1 ). When the ten times beam reduction by the parabola telescope is considered, the extremum is increased to β max/min ≈ ±314.1 mrad (or PFT ≈ ±1047 fs mm −1 ). By the substitution of Eq. (1) with the calculated extremum, for different (half) conical angles of α = 0.5°, 5°, and 10°, the variation range of the group velocity υ b is (0.997c, 1.003c), (0.976c, 1.033c), and (0.96c, 1.077c), respectively. We can find that the overall variation range of the group velocity is limited, and if required, which can be slightly increased by replacing DM with a free-surface mirror and/or increasing the magnification of the beam reduction telescope.

Discussion
The pulse-front deformation would shorten the propagation distance of a pulsed beam due to propagation diffraction, although the phase-front is corrected. When keeping all simulation parameters given in the "Methods" section and the setup shown in Fig. 1a unchanged, we remove the axicon and simulate the propagation of the pulsed beam around the image relay position. Figure 5a shows the optical fields at different positions of z = 0, 0.005Z R , 0.01Z R , and 0.015Z R (where z = 0 is the image relay position and Z R = 4 m is the Rayleigh length of the corresponding monochromatic Gaussian beam) and no serious diffraction distortion is found. When the Gaussian pulse bandwidth (FWHM) is increased from 10 to 20 nm, Fig. 5b shows after z = 0.005Z R the propagation diffraction seriously distorts the pulsed beam. Keeping the 10 nm bandwidth unchanged, when the PV value of the pulse-front deformation is enlarged from 300 to 600 fs, Fig. 5c shows that also after z = 0.005Z R the propagation diffraction seriously distorts the pulsed beam. In this case, for a broadband pulsed beam or a large-value pulse-front deformation, the propagation distance is limited by propagation diffraction. In this article, the length of the superposition region is around 115 mm (i.e., the beam waist 1 mm divided by the (half) conical angle 0.5°). If the pulse-front deformed pulsed beam is imaged into the center of the superposition region of z = 0, as shown in Fig. 5a, within the propagation range of (−0.015Z R , 0.015Z R ) the spatiotemporal distribution has no serious distortion, which can cover the whole superposition range. However, if the pulse bandwidth or the pulse-front deformation is enhanced, the flying length of the optical wave-packet is reduced. If the (half) conical angle α is enlarged to increase the variation range of the group velocity υ g (see Fig. 4b), the length of the superposition region is dramatically reduced, which relaxes the limitation induced by the propagation diffraction. We should also emphasize that in theory, the ideal Bessel beam can propagate over an infinite distance without any spread. However, in a real experiment, due to a finite beam aperture, the propagation invariant length is restricted. In this article, the pulse-front deformation is added onto the input pulsed beam, which enhances the propagation diffraction and accordingly shortens the propagation invariant length. In this case, the optical wave-packet slightly diverges during propagation, although both the quantitative analysis in this paragraph and the simulations in Figs. 2a, c and 5a show that the diffraction distortion within the conical superposition region is small.
Using the spatiotemporal coupling to control the group velocity and acceleration of an optical wave-packet or a focused intensity peak recently is an interesting and valuable technology. The spatiotemporal spectrum method invented by Abouraddy et al. 65 can change the group velocity of the optical wave-packet in a very large range by changing the tilt angle of the plane (k z , ω/c) with respect to the light-cone (k x , k z , ω/c) in spectral-space. Group velocities varying from −4c (in the backward direction) to 30c (in the forward direction) were measured, and in theory, arbitrary group velocities can be generated. The spatiotemporal dispersion method simultaneously demonstrated by Quéré et al. and Froula et al. can also adjust the group velocity of the focused intensity peak within a large range by changing the longitudinal chromatism and the temporal chirp [69][70][71] , and −0.09c to 39c flying focuses were measured in experiments. This work of the pulsefront deformed Bessel beam generation provides a third spatiotemporal coupling method to control the group velocity and acceleration 72 . Compared with the previous two methods, it provides an opportunity to precisely control the variations of both the group velocity and the acceleration of the optical wavepacket, although the tunable range of the group velocity is limited by the amount of the pulse-front deformation. Because the diversity of the pulse-front deformation makes the diversity of the group velocity (also the acceleration) possible, we can create some optical wave-packets with unusual motion forms, for example, the optical wave-packet with a nearly programmable group velocity (or acceleration) theoretically demonstrated in this article.
In conclusion, we have theoretically demonstrated an optical wave-packet having a nearly programmable group velocity by introducing a complex axisymmetric pulse-front deformation into the traditional Bessel beam generation. Different from the previous results of optical wave-packets displaying only a single motion form (superluminal or luminal or subluminal for the case of group velocity; and accelerating or uniform-motion or decelerating for the case of acceleration), the optical wave-packets here can propagate with nearly programmable motion forms during a single propagation path (e.g., superluminal followed by subluminal for the case of group velocity; and accelerating followed by decelerating for the case of acceleration). In this article, due to Fig. 4 Controllability of instantaneous group velocity and required SLM. a Pulse-front tilt (PFT) or tangent angle β at the intersection of two pulse-fronts during conical superposition. b Dependence of instantaneous group velocity υ b on PFT for (half) conical angles of α = 0.5°, 5°and 10°, and group velocity υ b is normalized by c. c Three different pulse-fronts of DM and required phase-corrections of SLM. d Dependence of modulation period number (within 1 mm) and required pixel size (five pixels within one modulation period) of SLM phase-correction on PFT. For a 4 μm SLM pixel size, when PFT is larger than 280 fs mm −1 , beam reduction by telescope is required. DM deformable mirror and SLM spatial light modulator. a periodically distributed pulse-front deformation along the transverse axis, the optical wave-packet propagates with superluminal and subluminal group velocities periodically along the longitudinal axis, and the corresponding instantaneous acceleration also varies between negative and positive values periodically. In this case, the propagating time-dependent motions can be well controlled by carefully optimizing the shape of the pulse-front deformation, creating optical wave-packets with unusual motion forms. As regards applications, we believe it can be used in some propagating velocity matched experiments, such as bio-imaging, particle-manipulation, particle acceleration, and radiation generation [75][76][77] , and the high spatiotemporal controllability could also offer new opportunities for fundamental studies in optics and physics.

Methods
Simulation model. The ideal thin axicon transfers a plane wave into a conical wave and generates the Bessel beam in the superposition region. Due to the axisymmetric distribution, the description is carried out in the 2D lateral plane containing the propagation axis for simplification, for example the x−z or y−z plane shown in Fig. 1a, and the divided half beams have individual traveling directions symmetrically bout the propagation axis. Both the input beam and the generated Bessel beam are described in the coordinate system of r−z, where r is the radial axis and z is the propagation axis. The divided half beam after the axicon is described in its own propagation (rotated) coordinate system of r α −z α , where r α is the radial axis and z α is the propagation axis. The origins of two coordinate systems have the same location at the geometrical center of the superposition region. When the clockwise rotation of the (half) conical angle α induced by the axicon is defined as the positive, and two coordinate systems of r−z and r α −z α satisfy the rotation relationship Under the paraxial approximation and the plane wave approximation, at the center of the superposition region (or the coordinate origin), the spectral optical field of a divided half beam is described in its rotated coordinate system of r α −z α and given by where A(r α ) and A(ω) are the spatial and spectral profiles of the amplitude, respectively, A SLM (r α ) is the spatial amplitude modulation by SLM, ϕ DM (r α ) and ϕ SLM (r α ) are the phase-modulations by DM and SLM, respectively. In this article, because of the image relay by the parabola telescope, we assume the shaped pulsed beam appears at the image relay position of z α = 0. The spatial amplitude modulation A SLM (r α ) by SLM is due to the net-like pixel−pixel gaps, which is described as where N is integer, rect() is the rectangular-function, d is the SLM pixel−pixel gap, p is the overall size of the SLM pixel (including the gap d), and then p−d is the effective size of the SLM pixel. Because DM has a continuous surface, it has no spatial amplitude modulation and only introduces a continuous phase modulation across the beam aperture. In this article, two axisymmetric cosine-function-like phase-modulations ϕ DM (r α ) by DM are used, respectively, and given by where k is the wavenumber, and L and D are the longitudinal PV value and the transverse period of the modulation, respectively. SLM in theory needs to introduce a conjugated [0, 2π) phase-correction for restoring a plane phase-front. However, the influence of the spatial resolution limited by the pixel size should be considered, and then the phase-front correction ϕ SLM (r α ) by SLM is described as where ⌊⌋ is the floor-function that gives as output the greatest integer less than or equal to the input, mod() is the modulo-function to return the positive remainder of a division, and k 0 is the wavenumber for the center wavelength (corresponding to SLM). The floor-function describes the spatially discrete phase induced by the SLM pixel size p, the modulo-function describes the [0, 2π) phase-variation by SLM, and finally the conjugated phase-correction is produced and moved to the region of (−π, π]. The angular spectrum method is used for modeling the propagation diffraction, and the optical field after z α propagation is given by where, at the initial position of z α = 0, the plane-wave angular spectrum and the corresponding optical field satisfy the Fourier-transform relationship After propagation, by using Eq. (2), the spectral optical fields of two divided half beams are described in the non-rotated coordinate system of r−z, and the coherent superposition (or interference) is given by E OWP r; z; ω ð Þ¼E r; z; ω; α ð ÞþE r; z; ω; Àα ð Þ : The temporal optical field is obtained by the Fourier-transform in spectrum and given by Finally, because the pulsed beam is axisymmetric about the propagation axis (z-axis), the 3D distribution is achieved by rotating the 2D result about the z-axis with r 2 = x 2 + y 2 .
Simulation parameters. Throughout this article, the parameters of the setup shown in Fig. 1a are as follows: the input Gaussian pulse has a 800 nm center wavelength and a 10 nm (FWHM) bandwidth; the beam diameter before and after the parabola telescope is 20 and 2 mm, respectively; DM introduces a cosinefunction-like axisymmetric deformation with a 10 mm period and a 90 μm (300 fs) PV value; SLM has [−π, π) (corresponding to the center wavelength) phasecorrection capability and a 40 μm pixel size including a 5 μm pixel−pixel gap; the parabola telescope introduces perfect ten times beam reduction; and the ideal thin axicon introduces a α = ±0.5°(half) conical angle. DM and SLM are positioned much closed to each other, and the beam divergence after DM is neglected. The parabola telescope images the pulse-front deformed pulsed beam into the geometrical center of the superposition region formed by the ideal thin axicon, where each divided half pulsed beam is assumed to have a six-order super-Gaussian beam profile in the lateral plane. All simulations are accomplished in the time domain with 1.2 ps window size and 1 fs accuracy (and the corresponding spectrum domain satisfies the Fourier relationship) and in the space domain with 100 mm window size and 2 μm accuracy (and the corresponding angular spectrum domain satisfies the Fourier relationship).
Group velocity and acceleration calculation. Under the approximation with a distortion-free pulse-front in a finite propagation length, the equations for the instantaneous group velocity and acceleration of the optical wave-packet are derived. Figure 6 shows, from the propagating time t 0 to t, the propagation distance along the z α -axis in the rotated coordinate system of r α −z α is given by where t 0 is the initial time when the phase-fronts arrive at the back-surface of the ideal thin axicon, and the propagation distance of the intersection of the phasefronts (or the intersection of the phase-front and the z-axis) along the z-axis in the non-rotated coordinate system of z−r is given by The propagation distance of the optical wave-packet along the z-axis in the coordinate system of r−z is given by where Δz(t) is the instantaneous longitudinal gap between the intersection of the phase-fronts and the optical wave-packet, satisfying The function z = f t0 (r) is the distribution of the deformed pulse-front in the coordinate system of r−z at the initial propagating time t 0 . In this article, the initial pulse-front distribution in the rotated coordinate system of r α −z α can be obtained by dividing Eq. (5) with the wavenumber k, and then, using Eq. (2), which in the non-rotated coordinate system of r−z can be obtained conveniently.
The instantaneous group velocity of the optical wave-packet is given by and by the substitution of Eq. (15) with Eq. (13), the instantaneous group velocity can also be described as Equation (16) shows the instantaneous group velocity is relevant to two terms: the first term is the constant velocity of the intersection of the phase-fronts; the second term is the variable velocity related to the change of the longitudinal gap Δz (between the intersection of the phase-fronts and the optical wave-packet), which is eventually dominated by the pulse-front deformation. This indicates why the pulsefront deformation can change the group velocity of the optical wave-packet and conforms to the phenomenon shown in Fig. 2. By the substitution of Eq. (16) with Eq. (14), the instantaneous group velocity can be re-described as where f t0 ′ (r) is the first-order derivative of the pulse-front function f t0 (r). The instantaneous acceleration of the optical wave-packet satisfies and by the substitution of Eq. (18) with Eq. (17), it is described as where f t0 ″(r) is the second-order derivative of the pulse-front function f t0 (r).
Equations (17) and (19) show, for a certain pulse-front deformation of z = f t0 (r), the instantaneous group velocity and acceleration of the optical wave-packet produced by this method can be directly calculated and, more importantly, well designed.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. Fig. 6 Propagation of phase-and pulse-fronts. Phase-and pulse-fronts propagate along z α -axis in the rotated coordinate system of r α −z α at the light speed in vacuum of c, and from propagating time t 0 to t, propagation distance is L 1 . Optical wave-packet propagates along z-axis in the coordinate system of r−z. Due to deviation between pulse-and phasefronts, from propagating time t 0 to t, propagation distance of the intersection of phase-fronts (or the intersection of phase-front and z-axis) L 2 and that of optical wave-packet (or the intersection of pulse-front and zaxis) L 3 are different, showing a variable-group velocity. t 0 is the initial propagating time when phase-fronts arrive at the back-surface of the axicon, and t is an arbitrary propagating time.