Manipulation and control of 3-D caustic beams over an arbitrary trajectory

We present an algorithm for manipulating and controlling 3-D field patterns, with energy confined to the narrow vicinity of predefined 3-D trajectories in free-space, which are of arbitrary curvature and torsion. This is done by setting the aperture field’s phase to form smooth caustic surfaces that include the desired trajectory. The aperture amplitude distribution is constructed to manipulate both the on-axis intensity profile and the off-axis beam-width, and is updated iteratively. Once the aperture distribution is calculated, the radiation from a finite sampled aperture is computed numerically using a Fast Fourier Transform-based scheme. This allows for both verification of the design and examination of its sensitivity to parameters of realistic discrete implementation. The algorithm is demonstrated for the cases of an Airy beam of a planar trajectory, as well as for helical and conical-helical trajectory beams. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

Several methods have been proposed, aiming at tailoring the on-axis intensity profile of curved beams. In [17], analytical solutions of Maxwell's equations with losses, propagating in absorbing media, were presented. These wave objects maintain their peak intensity along the beam-axis for large distances. Arbitrarily tailoring of the intensity distribution along the beam-axis of two-dimensional (2-D) surface-plasmon beams was implemented in [28]. The generation of such field patterns requires the design of source field distributions, often defined over wide apertures. For example, strong intensity profiles can be obtained at focal points of a bundle of geometrical optics (GO) rays (Bessel-type beams) [29,30] or along a smooth caustic of GO rays [6,7,[31][32][33]. For these cases, aperture field distributions (phase and amplitude) that give rise to the desired pattern can be calculated analytically.
Several methods have been proposed for shaping the accelerating beam trajectory and intensity profile. In [34], non-paraxial accelerating beams over arbitrary convex trajectories were designed. Some control over the intensity profiles along the caustic was obtained by changing the aperture field amplitude. In [35], independent control of both the trajectory and the maximum amplitude along it was achieved for 2-D caustic beams (CBs) using catastrophe theory. A practical algorithm for designing 2-D caustic beams, that meander along trajectories, consisting of both convex and concave segments, was introduced in [32]. The implementation was based on the generation of a caustic at the desired trajectory. This was done by "back-tracing" rays from the predefined beam trajectory to the source aperture, to determine the source distribution's phase. The aperture amplitude distribution was set to form a uniform intensity of the produced CB along its trajectory. These field design approaches are limited to particular 2-D cases with limited control over the beams' parameters. Yet, the principles used in these works inspire the design of more general field patterns that are of intensity confined to the vicinity of arbitrary three-dimensional (3-D) trajectories. Specifically, we are interested in the design of aperture distributions that produce patterns with 3-D trajectories of various curvature and torsion parameters. Moreover, along these trajectories, we would like to control the intensity profile and the beam-width.
In this work, we present a 3-D extension of the algorithm in [32], for the design of beam field patterns, along desired beam-axis trajectories with arbitrary curvature and torsion (as defined in section 2). In the proposed scheme, confining the field to the vicinity of the predefined beam-axis trajectory is done by designing the aperture field's phase, to give rise to a smooth caustic surface (section 3.1) that formed around the beam-axis "skeleton". Following the 2-D methodology (and using the definitions in section 3.2), we design the aperture amplitude distribution, with the aim of manipulating the on-axis intensity profile (section 3.3). In the 3-D case, another controllable feature of the pattern is the off-axis beam-width. It can be manipulated via the corresponding aperture distribution's effective off-axis width (section 3.4). The proposed technique is validated via numerical experimentation; Once the aperture distribution is calculated, the radiation from a finite sampled aperture is computed numerically using a Fast Fourier Transform (FFT)-based scheme (section 4). This technique also enables the examination of the design's sensitivity to parameters of realistic discrete implementation, in an experimental set-up, for example by using a finite emitters array. The algorithm is demonstrated for the cases of Airy beam planar trajectory, as well as for helical and conical-helical trajectories, in sections 5.1, 5.2, and 5.3, respectively.

Problem definition
Given a trajectory in a 3-D domain with σ ∈ [σ min , σ max ] denoting the arc length along the beam-axis trajectory, we aim at designing a source (aperture field) in the z = 0 plane that produces a 3-D beam-field, with a beam-axis that follows this curved trajectory, in free-space. The trajectory can be intricate, with arbitrary curvature and torsion. In this work, we assume that the characteristics of the trajectory are large on a wavelength scale and limit the discussion to trajectories with tangent vectors of positive z-component. The beam-field's desired transverse width is denoted W b (σ) (see Fig. 1) and I b (σ) denotes the desired intensity profile along the beam-axis. Our design approach relies on the observation that the radiated field's intensity can be confined to the vicinity of a given trajectory r b (σ), if the aperture phase, asymptotically, gives rise to a caustic surface that contains r b (σ). In other words, the approach produces a 3-D CB along the trajectory. To that end, the designed aperture field, denoted u 0 (x , y ), is sought in the form of a ray-field where k = ω/c is the wavenumber, with c being the free-space wave velocity. Also in (2), (x , y ) are coordinates in the z = 0 plane, and A and S are identified as the aperture field's amplitude and phase (Eikonal) distributions, respectively. In what follows, we present a method for designing the two functions, A and S, given the desired beam-axis trajectory, r b (σ), intensity, I b (σ), and beam-width, W b (σ).

Aperture field design
The design of an aperture that gives rise to a CB of desired properties consists of three parts: first, an aperture phase distribution is designed, in accordance with the desired CB's beam-axis trajectory. It is worth noting that, for planar trajectories, this part reduces to application of the 2-D technique in [32], for computing the phase along the trajectory's projection onto the z = 0 plane, followed by an extension of the phase in the normal direction. However, for treatment of arbitrary 3-D trajectories, both steps of forming the phase distribution must be generalized. Then, the aperture amplitude distribution is designed to control (i) the CB's on-axis intensity profile and (ii) the CB's width. While the phase distribution is computed only once per given trajectory, the aperture amplitude distribution can be further adjusted, iteratively.

Beam-axis manipulations
The proposed method attempts to form a caustic on a surface containing the beam-axis r b (σ). This is done by setting the aperture phase to give rise, at each point on the aperture, to a ray that participates in forming the caustic surface. First, we construct the phase that forms the beam-axis, which can be though of as the caustic "skeleton". Refereing to Fig. 1, from each point r b (σ) along the beam-axis, we follow the tangent ray back to its emanating point on the aperture, denoted r a (σ). The ray angles formed between the ray and the corresponding axes, θ x,y , are obtained via the unit vector tangent to the parametric trajectory in (1). i.e., In (5), the over dot denotes a derivative with respect to the argument, i.e., z(σ) = dz(σ)/dσ. Referring to Fig. 1, we identify the tangent ray trajectory emanating from point (x , y ) on r a (σ) as This procedure maps points on the axis to points on a curved trajectory r a (σ) in the aperture as depicted in Fig. 1.
Following this, we define the local curve coordinates, (σ , n ), on the aperture plane, where σ is the arc length along the curve r a (σ ) and n denotes the distance of a point (x , y ) from r a (σ ) (in the normal direction), with corresponding unit vectors,σ andn (see Fig. 2).
Here and henceforth, we use hat over bold fonts to denote unit vectors. Thus, the parameter σ ∈ [σ min , σ max ] is mapped to σ ∈ [0, σ max ] via (4). This mapping is denoted by σ (σ). From each point (x , y ) on the curve r a (σ ), the ray that forms the caustic over the beam-axis trajectory emanates in the direction of the unit vector where cos θ z (σ ) = 1 − cos 2 θ x (σ ) − cos 2 θ y (σ ). we project the ray directionŝ onto the unit vectorsσ andn namely, the tangent and normal to the curve r a . The aperture phase distribution is first constructed along r a via (7) and extended via (9) to other points with in the marked beam width domain w t in (12).
Thus, the phase along the curve r a (σ ) can be reconstructed via Here, ∇ S is evaluated by substituting θ x,y in (3) into (6) and dσ = dσ σ (σ being the tangent unit vector along r a (σ ) (see Fig. 2)). Next, the phase is extended to other points on the aperture. The extension is not unique and it dictates the caustic surface's shape. First we consider the phase extension for the special case of planar beam-axis where y(σ) = 0 in (1). In this case, cos θ y (σ ) = 0 and therefore, using (6), the normal phase derivative ∂S/∂n = ∂S/∂y = 0 for all points along r a (σ ). Any extension of the phase S(σ , n ) is subject to this boundary condition. The simplest form of extension is the linear one in which ∂S/∂n = 0 for all points over the aperture. This choice will reproduce the caustic curve for rays in any constant y plane, forming a ∂/∂ y = 0 caustic surface (within the beam width domain), see Fig. 3 (though our method of phase construction from tangent rays guarantees the formation of a caustic surface, this can be easily verified by evaluating the ray Jacobian using (3) and setting its determinant to zero). Other continuous extension, subject to the boundary condition, are possible but yield intricate caustic surfaces (a special interesting limit case is the conventional 3-D Airy beam where the caustic takes the form of a cusped double-layered caustic -see [9]). The phase extension for the general (non-planar) beam-axis is performed in a manner similar to that of the planar beam extension: since the integration in (7) reconstructs the phase along the curve r a (σ ), it sets the directional derivative along the curve to dS/dσ = ∇ S ·σ . In order to satisfy (6) for points over r a (σ ), the directional phase derivative normal to the curve should satisfy the boundary condition withŝ as defined in (5) (see Fig. 2). By using (8), and following the planar beam-axis phase extension, we extend the phase S(σ , 0) in (7) for small n values via a linear extension of the form S(σ , n ) = S(σ , 0) + n ŝ(σ ) ·n (σ ).
The phase in (7) and (9) forms a smooth caustic surface that passes through the beam-axis. The normal extension in (9) is carried out for n values for which the aperture field amplitude is not null (the amplitude reconstruction is discussed in the next sub-section). Here, it is assumed that the phase at each contributing point can be assigned uniquely, in accordance with the point's distance to the aperture curve. This requirement is satisfied for n values smaller than the local radius of curvature of r a (σ ). This implies that the transverse aperture amplitude distribution must be sufficiently localized (i.e., with energy confined to a finite width around the axis, denoted w t in (12)). Figure 3(b) shows the caustic surface for a 3-D Helical trajectory, similar to that described in detail in the results section. The surface was evaluated by setting the determinant of the ray Jacobian to zero for all rays that are emanating from the aperture within the beam domain (in the vicinity of r a (σ )).

Local beam coordinates
In order to manipulate and control the on-axis intensity and the beam-width, we apply the local beam coordinates that are associated with the beam-axis in (1). Denoting σ the arc length along the beam-axis trajectory, the beam local coordinates are defined by the unit vectorst(σ),n(σ), andn b =t ×n, being the tangent, normal, and bi-normal to r b directions, respectively at a point r b (σ) on r b . They are related by the Frenet-Serrer equations where K(σ) is the curvature of r b and κ(σ) is its torsion. Note that the beam coordinates are non-orthogonal for κ 0. A locally orthogonal coordinate system along the beam may be obtained by transverse rotation of the unit vectors [37]. Points near the beam may now be expressed as By using these definitions, we can discuss the controlling of the beam intensity profile along each of the three local coordinates, σ, n and n b . Since r b is a trajectory over a caustic surface, we identify the mechanism of enhanced intensity along the local coordinate n as the caustic local structure. This implies that, asymptotically (and paraxially), the intensity profile in the direction normal to the caustic surface is, as in any simple caustic, of the form of the (squared) Airy function profile (not to be confused with the Airy beam trajectory). Therefore, in the proposed method, control over the width of the beam intensity in the localn direction is not possible. However, the beam-width in the transverse n b coordinate and the intensity profile along σ are controllable. In what follows, we propose a simple technique for separately controlling each of the two. To that end, in our implementation, we design the aperture field amplitude in (2) to be of the general form We term w t and w l the transverse filter and longitudinal window, respectively. In the next sections we discuss the properties of these filter and window and their usage for controlling features of the beam intensity.

On-axis beam intensity profile manipulations
The longitudinal (amplitude-equalizing) window w l (σ ) controls the field intensity profile along the beam-axis. The desired intensity profile, I b (σ), can be any continuous function of σ, pending on the application. Finer adjustments of the on-axis intensity profile can be carried out iteratively, after setting the transverse filter according to section 3.4. In this procedure, we start by setting In each (n th ) step, the aperture field is propagated into the z>0 half space (by using the method described, later, in section 4). The field's intensity, denoted I (n) l (σ), is sampled along the beam-axis. We then use the mapping σ (σ) and adjust the amplitude-equalizing window via multiplication by the (real) function In (14), β (n) l ∈ [0, 1] is a feedback parameter that is chosen to facilitate convergence. This can be carried out in several ways: For example, a line search can be performed to find β (n) l that minimizes the beam-width error in the L 2 norm sense, with respect to the desired one. A simpler approach would be to maintain β (n+1) l = β (n) l , unless the error increases in the n th step, in which case the value of β (n+1) l is reduced. Examples for this procedure are presented in section 5. Following this approach, the amplitude can be designed to be roughly uniform, by setting I b (σ ) = 1 for all σ ∈ [0, σ max ] (see examples in section 5).

Off-axis beam-width manipulations
The tapering of the beam intensity along the coordinate n b is controlled by the transverse filter, w t (n , σ ) in (12) -a real function of finite support centered about n = 0. There is a trade-off between the desire for a sharp decay of the field intensity away from the beam-axis and the diffraction caused by abrupt changes in the aperture field. Therefore, the window profile, which determines the intensity profile parameters (width, decay rate, diffraction artifacts, etc.), should be selected carefully. Good candidates for the transverse filter, w t , include the Gaussian windows, as well as other gracefully decaying window functions.
In the examples presented in section 5, we use a scaled Tukey window, in order to minimize the diffraction. This window is given by where T denotes the standard definition of the Tukey window with the parameter 0<r 1, which sets the half amplitude width of T(t; r) to r/2 (see Fig. 4). In (15), the transverse window-width W t (σ ) varies along r a (σ ) (see also Fig. 2). The local source distribution width, W t , controls W b (σ) at the corresponding points along the trajectory and is adjusted iteratively, in order to obtain the desired beam width W b (σ), in the following manner: we first use the mapping σ (σ) following (3) and set W t (σ ) = W b [σ (σ)]/e a , where e a denotes the aperture efficiency, which, in the case of the Tukey window, is given by e a (r) = Then, to improve the control of the beam width, the filter-width, W t (σ ), is updated iteratively, in a manner similar to that in section 3.3: In each step, after calculating the field intensity profile in z>0 (as be described in section 4), we evaluate the resulting beam-width along the beam-axis, denoted by W (n) t (σ). Next, we assign the beam-width to the aperture curve r a (σ ) and adjust the window width W t (σ ) in (15) via multiplication by the (real) filter Here, β (n) t ∈ [0, 1] is a feedback parameter, similar to β (n) l in section 3.3.

Validation
In order to verify the proposed CB synthesis approach, the fields radiated by the aperture distribution are computed in the half-space z>0. Ideally, this should be done by applying the Kirchhoff-Huygens integral to the aperture distribution at z = 0, such that In (19), ∂ n = ∂ z denotes the normal derivative of the 3-D free-space Green's function Note that the infinite integral is truncated according to the finite aperture's support of size L x × L y . However, analytical calculation of (19) for an arbitrary u 0 is neither practical nor it is representative of a realistic implementation. In practice, it can be assumed that any realistic aperture implementation will be inherently discrete [1,38], e.g., a spatial light modulator or a patch antenna array, with uniform spacing between elemental sources. The field produced by these sources will be calculated on a regular Cartesian grid. Therefore, the infinite field integration is replaced with a discrete summation and the radiated field can, thus, be written as In (21) , ∆x and ∆y are the grid spacings, (x , y ) = (m ∆x , n ∆y ), M = L x /(∆x + 1) and N = L y /(∆y + 1) are the numbers of aperture samples in each coordinate. If the field is evaluated, at a given z, also on a regular x − y grid, with the same ∆x and ∆y, the summation can be written as a discrete convolution of the sampled u 0 and ∂ z G, such that Here, u[m, n, z], u 0 [m , n ], and ∂ z G[m − m , n − n , z] are discrete functions, consisting of the sampled u(r), u 0 (x , y ), and ∂ z G(r − r ), respectively. The computation of (22) for all points on the observer grid is accelerated by using a 2-D FFT -based algorithm, similar to that in [39]. The 2-D FFTs of u 0 [m , n ] and ∂ z G[m, n, z] are first computed, and are then multiplied in the spectral domain, and the result is transformed back to spatial domain by using a 2-D inverse FFT to obtain u[m, n, z]. Note that, for an aperture that is modeled as an array of uniformly spaced point sources and measurements on a regular grid, this fast scheme for computing the field is exact (that is, in contrast to the z-propagation of sampled plane-wave spectra of the aperture fields).
In this work, we model the aperture as an array of such point sources, weighted by the sampled values of the aperture distribution. We use a typical discretization length of ∆x = ∆y = λ/10. For more accurate modeling of a particular realistic implementation, one could replace ∂ n G(r − r ) by a function describing the field produced by an element of the array. This field representation enables the study of the proposed synthesis approach's sensitivity to implementation specifications dictated by realistic system limitations, e.g., array length, element spacing, and elemental sources' field [4,11]. Finally, it should be noted that the need for fast radiation integral computation is even more pronounced when repeated field evaluation is required. That is the case when an iterative procedure for adjusting the beam's (on-axis) intensity and (off-axis) width in sections are used, as in sections 3.3 and 3.4.

Examples
In this section, the proposed method is demonstrated for several representative examples. The approach is first validated for the simple case of a planar (beam-axis) trajectory. This case is geometrically very similar to the 2-D case, in terms of the design of aperture phase distribution. Therefore, it enables us to conveniently study the sensitivity of the design to the aperture's discretization. Then, we present results for more complex cases of helical and helical-conical CBs. These cases challenge the proposed scheme and allow us to showcase its advantages.

Planar beam: analytical example
For the first example, we chose an Airy-type planar (torsion free) CB, for which the beam-axis is defined by the parametric trajectory This beam-axis describes a planar trajectory on the x = y plane, where z b (t) = α √ t is identified as the well-known Airy beam trajectory. In this example, we synthesize a constant-width beam with W b = 15λ for α = 5. By using (3)-(6) we evaluate for r a (σ ) By using (24) in (7), we evaluate the phase along r b (σ ) as In the special case of a planar trajectory, ∂ n (σ ) = 0. By using (8), we extend the phase S(σ , n ) = S(σ , 0). In this example, the beam-axis in (23) is defined for σ ∈ [σ min , σ max ]. In order to avoid diffraction effects at the end of the beam-axis, r b (σ max ), and to obtain the desired beam-width there, we extend the beam-axis to the interval σ ∈ [σ min , 1.25σ max ] and calculate the phase and amplitude of the sources accordingly. This extension is used in the following two examples as well.
Following the procedure in section 3.4, the initial CB aperture distribution is obtained by setting W t = W b /e a = 30λ with r = 1 in (15). By using a grid spacing of λ/10, we obtain the aperture field distribution in Fig. 5. The corresponding CB-field intensity is evaluated according to section 4. Figure 6(a) depicts the field intensity over the x = y plane, demonstrating that the on-axis intensity profile follows the desired trajectory (marked by the black dashed curve). Figure 6(b) plots the projection onto the z = 0 plane, of points in the half-space z>0 of field intensity greater than half the maximum at their respective z-plane. The solid black line indicates the desired beam trajectory, whereas the dashed black lines to each side of the trajectory indicate the desired -3dB width. This point of view shows clearly the control of the beam-width in then b local direction.  This example is also a convenient test case for studying the sensitivity of the synthesized aperture distribution to constraints imposed by a realistic implementation. Specifically, the spacing between array elements that aim to mimic the continuous source distribution in the aperture may be larger than of the order of a wavelength. Such is the case, for example, of many existing spatial light modulators (SLMs) [1] [4][11] [40]. These devices for controlling the aperture phase are designed to operate at 300nm-1000nm wavelength, but often use a typical pixel edge length (also termed pixel pitch) that is, at best, currently, as small as 0.72 µm [41]. Antenna arrays are typically designed with a λ/2 spacing between elements. To evaluate the effect of aperture (under-) sampling, we examine the fields produced by aperture distributions sampled at various rates. For these fields, we calculate the standard deviations of the -3dB beam-width in the local normal axis. These are plotted in Fig. 7 as a function of the grid spacing (Note that, to maintain same grid spacing for both the aperture and the measurements, we use a fixed grid spacing of λ/20 (smaller than the numerical rule-of-thumb spacing of λ/10) for both aperture and measurement grids all cases, and for the coarser sampling cases, sparse filling of the aperture is used as needed). It can be seen that the ability to maintain the desired beam-width drops sharply for grid spacings greater than λ/2. This agrees with the well known λ/2-spacing requirement for antenna pattern design. Fig. 8 shows, for various cases of grid spacings, the points of maximal intensity for each z plane, alongside the desired beam trajectory. It can be seen that even when exceeding the λ/2-spacing limit, e. g., for ∆ = λ, the field pattern's maxima continue to follow the desired trajectory. It is only for spacings greater than 2λ that the pattern begins failing to follow the desired trajectory. In other words, in the σ direction, some under-sampling is tolerable, as the trajectory is dictated by the geometrical nature of the ray-field. That is not the case for sampling in the n direction. Aiming for the design of arbitrary trajectories, which are not necessary within a plane, the λ/2 restriction will hold for the entire aperture.

Helical beam
With the basic features of the proposed technique demonstrated, we now turn to examine its adequacy for the construction of CBs with more complex beam-axis trajectory (i.e., a trajectory with arbitrary curvature and torsion) and intensity profile. As our second example, we choose a CB that follows a right-handed helical beam-axis that is described by where R denotes the helix radius and the parameter P sets the pitch to be 2πP. The arc length of this trajectory is σ = (t − t min ) √ R 2 + P 2 . In this example, the desired beam-width profile was W b (t) = A[1 + sin 2π t−t min t max −t min /5] and the on-axis intensity profile was I b (t) = 1. By setting R = P = 20λ and following the procedure in section 3, we numerically calculate S(x , y ). Next, we calculate the amplitude A(x , y ) by setting r = 0.8 in (12), with the transverse filter w t (n , σ ) in (15) and the longitudinal window w l (σ ) in (13). Finally, we calculate the propagating fields in the z>0 half-space via (22). At every σ, the beam-width is defined as the Euclidean distance between the -3dB of the peak intensity points in the bi-normal local direction (see Fig. 9). Using the initial aperture design, we evaluate the intensity. In Fig. 10, the resulting iso-surfaces of -3dB of the peak intensity is plotted. The beam-width and the on-axis intensity profile are shown (red dot markers) in Fig. 10(c,d), respectively. Clearly the resulting intensity profile and beam-width deviate substantially for the desired ones (black, solid). In order to improve the result, we update the aperture field according to the desired on-axis profile and beam-width, following the procedures in 3.3 and 3.4. The results for the updated aperture fields are shown in Fig. 10(b) (for the -3dB iso-surface) and in Fig. 10(c,d) (red, dotted lines). Clearly, in both beam-width and intensity profiles, the updated aperture distribution meets the design requirement better.  In order to quantify the improved performance of the updated design in the on-and off-axis profiles, we use a criterion similar to the one introduced in [32], i.e., the deviation from the desired profile is assessed via an index, defined as the standard deviations of the obtained parameter minus the desired one. The computed off-axis and on-axis index are 2.87λ and 1.05, respectively for the initial CB, and 0.39λ and 0.03, for the updated CB. These numbers and Fig. 10, indicate the significant control capabilities over the desired features of the methods.

Conical-helical beam
For the last example, we chose a CB with a conical-helical beam-axis (a helix of gradually decreasing radius). This example sheds light over the aperture field structure obtained for complex r a curves. The desired beam-axis is given by for t min ≤ t ≤ t max . In (27), R is identified as the initial helix radius, P sets the pitch to be 2πP(t max − 0.8t min ) −1 , and t max sets the cone height to Pt max /(t max − 0.8t min ). In this example, we use R = 11.7πλ and P = 58.5πλ and set t min = 0.1π and t max = 4π.  The curve r a (σ ) on the aperture, corresponding to the beam-axis in (27), forms closed loops (see Fig. 11(a)). Here, we extend the approach proposed in [32] for the areas close to the intersection points. For these areas, two sets of rays need to be synthesized, each forming different sections of the caustic surface (see Fig. 11(b)). Therefore, the aperture field is a sum of ray fields of the form in (2). In the iterative process of adjusting the on-axis intensity profile or the off-axis width, we apply the procedure to each of the terms separately, according to the mapping discussed after (3). In this example, we choose W t (σ ) = 20λ and I(σ) = 1. The resulting -3dB iso-surface of the CB intensity after applying the on-axis and off-axis procedure is plotted in Fig. 12.

Conclusions
In this paper, we introduced a practical algorithm for the construction of aperture field distributions that correspond to 3D-CBs. The procedure involves the construction of the aperture field phase and amplitude, and application of an iterative algorithm to manipulate and control the intensity profile along the beam-axis as well as the beam-width profile. Numerical examples demonstrate the algorithm's capability to synthesize aperture field distributions that give rise to beams with arbitrary curvature and torsion. Along these trajectories, control of the on-axis intensity profile and beam-width in the bi-normal direction was demonstrated. The representation of the aperture and measured fields as samples on uniform Cartesian grids allowed for the convenient acceleration of the computation of exact field integrals. This was used for both validation of the approach and as a part of the iterative improvement of the design. A preliminary study shows that the sensitivity, in terms of beam-width deviations from the design specifications, to the aperture sampling rate, is similar to that acceptable in antenna array design. However, if only the trajectory formed by the peak intensity is of interest, the restrictions on the aperture sampling can be relaxed. This suggests that an experimental validation set-up should include the lowest pixel-pitch SLM available to date [41]. Future work may involve the design of optimization schemes to improve the control of the beam parameters manipulated in this work, i.e., the beam-width in the n b local direction and the intensity profile as well as for gaining control over the beam-width in the n local direction.