1. Introduction

Free-space optics, such as lenses, mirrors, and gratings, offer accurate spatial beam control and are easy to install in optical systems. Freeforms are particularly efficient and cost effective in achieving custom field transformations over a wide angular range [1]. Well-documented design algorithms are available for “zero étendue” applications and extended sources [2]. These design approaches are typically based on geometric optics and are most appropriate when diffraction or partial spatial coherence effects are negligible. Diffractive optical elements such as gratings and nanosurfaces achieve accurate far-field transformations of coherent beams but suffer from unwanted grating orders. In computer generated holograms the dominant errors arise from lateral sampling and phase discretization, whereas in geometric phase optics these limitations are absent and much higher intrinsic efficiencies are obtained [3,4]. However the overall efficiency of either type can still be low due to coupling into unused outputs (e.g. the conjugate wave of a zone plate) and length restrictions determined by the space-bandwidth product [5]. Multi-plane architectures, which apply a series of phase modifications to a propagating field, can perform arbitrary transformations with high efficiency [6,7].

Freeforms, gratings, and other free-space optics are used at a distance from their source and target. This makes them easy to integrate into larger systems, but essentially constrains the overall size to that of the largest element, typically millimeters or greater. For example, coupling between two mismatched waveguides often involves expansion to free space followed by focusing. The overall length and width of the coupler is determined by the lens dimensions, leading to “wasted” volume in which the field transitions to and from free space. In addition, the disparate size scales (several micrometers for a single-mode waveguide and several millimeters for a typical aspheric lens) leave the system vulnerable to misalignment and mechanical vibration.

In silicon photonics, beam shapers that change or taper in the direction of propagation are common [8]. A wide variety of mode converters have been designed to address the challenges of mode size and wavevector matching between silicon waveguides and lower-index media. One of the simplest structures is the inverse taper, in which the mode size grows as the waveguide width shrinks along the direction of propagation. Vertical [9] or lateral [10] arrangements of multiple tapers offer increased fiber coupling efficiencies, and longitudinal metastructures offer very high efficiency across a wide bandwidth at both TE and TM polarizations [11]. However, inverse tapers are restricted to very small lateral size scales and do not offer arbitrary control over the field profile. Consequently their use is largely limited to on-chip beam resizing at relatively low optical power.

In this paper, we propose and design volumetric GRIN elements that perform smooth and arbitrary mode conversion at lateral size scales of tens to hundreds of micrometers. The physical package is envisioned to be compact, with its lateral width matched to that of the beam. Further, since the entire transformation is completed inside the GRIN volume, it can be butt-coupled to the source and/or target, which may confer mechanical, thermal, and power handling advantages. For weak changes in refractive index, the GRIN length can be tens of millimeters, over which diffraction effects are significant. Accordingly, the refractive index profile is designed to control and compensate diffraction effects.

Fabrication techniques are available to realize arbitrary 3D refractive index profiles. Femtosecond laser writing techniques [12] such as two-photon polymerization [13] and direct writing into glass [14] are well-suited to the lateral feature size (several micrometers) and low overall refractive index contrast ($\Delta n\sim 10^{-3}$) needed for GRIN designs like the example shown in this paper. Since hours-long fabrication times are typical in scan-based femtosecond laser methods, more scalable fabrication methods may be envisioned, such as lithographic-style fabrication with a liquid photopolymer [15,16].

This paper builds upon a previous design method for 2D GRIN devices that operate on a single transverse dimension of a coherent field [17]. Whereas this older design method is suitable for fields that are reducible to one dimension, e.g. because the fields are separable or due to symmetry, it cannot address general 3D GRIN design problems in which the field may be an arbitrary 2D function. This paper presents a formal mathematical generalization to support semi-arbitrary 2D fields, within some constraints.

The aim of this paper is to describe the numerical design process for GRIN structures that transform laser beams efficiently at dimensions where diffraction is significant. In the next section, we explain the mathematics and provide details of each step of the overall method. Following this, an example GRIN design is shown which converts a Gaussian beam to a semi-arbitrary field pattern with the letters “GRIN.” Finally, limitations and advantages of the method are summarized.

2. Design method

In this section we describe the mathematical approach to designing custom GRIN mode converters. The design procedure consists of four steps:

2.1 Mapping function

Irradiance, as power per unit area, is a density function that can be manipulated by spatial redistribution of optical power in two-dimensional (2D) coordinates. New irradiance profiles can be created from a given one by moving energy from one location to another across the spatial domain. Mathematically the spatial redistribution can be described as a mapping function $f$ that converts between two sets of spatial coordinates, $f(\mathbf {x_0}) = \mathbf {x_1}$, where the original coordinates are written $\mathbf {x_0} = (x_0,y_0)$ and the mapped coordinates are written $\mathbf {x_1} = (x_1,y_1)$. Locally, the new coordinates $\mathbf {x_1}$ may be compressed or expanded relative to the original coordinates $\mathbf {x_0}$. Then an infinitesimal area $dA_0$ located at $\mathbf {x_0}$ changes size to become $dA_1$ located at $\mathbf {x_1}$ in the new coordinate plane. The factor by which this local area changes is

After writing an overall mapping function, the remainder of the design problem is to find a refractive index profile that implements this mapping. To do this we must first specify how the mapping evolves with $z$. We describe this in the next section.

2.2 Continuous deformation of the mapping

In the previous section we considered a mapping function $f$ that mapped coordinates in the input plane to those in output plane, $f(\mathbf {x_0}) = \mathbf {x_1}$. To define the field everywhere inside the GRIN volume and not just at these two planes, we seek a continuous deformation of the mapping that specifies new coordinates $\boldsymbol {\xi } = (\xi , \eta )$ as a function of axial position. Such a continuously deformed mapping is known as a homotopy in mathematics.

In this section we adopt a mathematical template for the homotopy and write the Jacobian in terms of $z$ and the derivatives of the mapping function. We first consider a single mapping and then a chain of mappings $f(\mathbf {x_0}) = f_n(\ldots f_2(f_1(\mathbf {x_0})))$ which can be used to build complicated maps from simpler ones.

2.2.1 Single mapping

A wide range of homotopies are acceptable from a mathematical point of view. The homotopy must only be continuous and satisfy the mapping function at the input and output planes. However, practically the choice of homotopy determines the field evolution throughout the GRIN medium and hence the refractive index profile that will ultimately be found. Many choices of homotopy are unphysical or undesirable.

In this work we do not seek to optimize the homotopy but rather choose the particular form

A similar cubic polynomial homotopy has been used in the design of two dimensional GRIN devices. Other forms may be more optimal in terms of optical characteristics such as reducing the overall refractive index contrast or avoiding focusing to high irradiance in the medium. But these desirable quantities are difficult to estimate a priori, and hence optimization would occur by forward comparison. For example, the cubic homotopy has been compared to one based on cosine sections for the design of 2D GRIN elements [19].

To obtain the irradiance at an arbitrary axial plane, we modify the input irradiance with the Jacobian. The Jacobian at an arbitrary axial plane is written

Note the derivatives in Eq. (9) are taken with respect to $\mathbf {x_0}$, which sets the reference plane at $z=0$ to assess the local area distortion. This choice is made without loss of generality. For instance, choosing $\mathbf {x_1}$ as the reference is equally valid since it can be found from $\mathbf {x_0}$ via the mapping function.

2.2.2 Multiple mappings

For complicated mapping functions, it will be convenient to write the overall mapping function $f$ as a chain of mapping functions $f(\mathbf {x_0}) = f_n(\ldots f_2(f_1(\mathbf {x_0})))$. This allows us to divide the overall mapping into parts and solve for sub-mappings using different methods for each individually. We now write the homotopy as a series of $n$ mappings

Each mapping will be accomplished by a cubic polynomial homotopy following the same form as Eq. (6), that is,

Coordinates for the $m^{\mathrm {th}}$ mapping and its homotopy are shown in Fig. 1. Lateral coordinates $\boldsymbol {\xi _m}$ at a plane of interest lie on trajectories $\boldsymbol {\xi _m}(\zeta )$ relating coordinates $\mathbf {x_{m-1}}$ and $\mathbf {x_m}$ at the end planes.

 

Fig. 1. Lateral coordinates as a function of normalized axial position $\zeta$.

Download Full Size | PDF

Our aim is to write the Jacobian derivatives in terms of the individual mappings referenced to the input coordinates. We proceed by differentiating Eq. (11) in both dimensions to obtain the Jacobian derivatives

The $\xi$ derivatives in Eq. (16) can now be rewritten in terms of mapped coordinates $\mathbf {x_m}$ using the chain rule. That is,

Next, the terms $\partial \xi _{i,m}/\partial x_{i,m-1}$ and $\partial \xi _{i,m}/\partial x_{j,m-1}$ in Eq. (17) can be found from the homotopy for the mapping $f_m$ because $\mathbf {x_{m-1}}$ plays the role of the initial position. In a similar pattern to Eq. (10), we write out the terms by substituting Eqs. (12) and (13) into Eq. (17) and performing differentiation to yield

2.3 Evolution of the complex field

To obtain a 3D refractive index profile, we must first specify the magnitude and phase of the complex field throughout the GRIN medium. In this section we first describe how to obtain the magnitude, then how to estimate the phase.

With the mathematical developments of the previous sections, we can write out the irradiance profile at any axial position using the Jacobian. Equation (5) applies generally to any set of mapped coordinates and may be rewritten as

The phase of the complex field can be found by wavefront construction based on the homotopy. The homotopy describes not only a distortion of lateral space but also a trajectory $\boldsymbol {\xi }(\zeta )$ for every mapped coordinate as it is brought through the length of the GRIN medium. Further, since $\boldsymbol {\xi }(\zeta )$ is continuous, the axial derivative $d\boldsymbol {\xi }/dz$ is also continuous.

To construct the wavefront, we define a collection of unit-length normal vectors $\mathbf {s}$ that are tangent to the trajectories $\boldsymbol {\xi }(\zeta )$ at a plane of interest. Each unit normal has directional components $\mathbf {s} = s_x\mathbf {\hat {x}}+ s_y\mathbf {\hat {y}} +s_z\mathbf {\hat {z}}$. Assuming that the wavefront varies slowly with $z$, then its surface is well approximated by $\mathbf {s}(\boldsymbol {\xi })$ at an axial plane. This effectively swaps movement along the trajectories $\mathbf {s}$ with movement along the $z$ axis, which is a paraxial approximation. In accordance with this approximation, the gradient of the phase $\phi$ is written

2.4 Refractive index recovery

With the magnitude and phase of the complex field now fully defined throughout the GRIN volume, the problem is well posed for refractive index recovery. To perform the calculation we adopt a model of the propagation physics and run an algorithm that finds a refractive index function consistent with the specified field evolution.

Until this point in the design problem we have been performing mathematical manipulations to specify the complex field. In principle these are separate from the physics and can be repurposed when other propagation models are used. Refractive index recovery based on the eikonal equation has been used to design GRIN optics according to a geometric optics model [19]. In the present work we use the split-step beam propagation method (BPM) to model the propagation physics and a phase retrieval method to recover the refractive index. This approach has been used previously [17] to design two-dimensional GRIN structures but can now be used to address fully 3D design problems.

We begin by describing the split-step BPM (treated in Ref. [21], for example). The model represents propagation through an inhomogeneous medium as a series of short axial steps of length $\Delta z$ and considers the scalar wave equation in terms of distinct operations that occur in each step. Local variations in refractive index are written as phase modifications $\theta (x,y)$ accrued over the discrete propagation length. Specifically, for paraxial waves traveling through a medium with refractive index function $n(\boldsymbol {\xi },z) = n_0 + \Delta n(\boldsymbol {\xi },z)$, the phase modification is simply $k \Delta z$ times the refractive index variation $\Delta n$,

We now describe phase retrieval. Phase retrieval is used in various contexts to recover the phase profile of two complex functions that are related by a Fourier transform or Fresnel transform. In the present problem we will use it to find the phase modifications $\theta (\boldsymbol {\xi },z)$ which determine the refractive index profile. The mathematical setup is illustrated in Fig. 2. Two consecutive propagation steps are shown, where each step implements the sequence of operations described in Eq. (24). Fields at the same axial position share the same magnitude and differ only in phase, e.g. $|u_B|=|{u_B}'|=|{u_B}''|$. Pairs of phase modifications between these fields are found by phase retrieval at both ends of each axial step.

 

Fig. 2. Axial locations of fields and phase modifications involved in phase retrieval.

Download Full Size | PDF

For instance, $\phi _{i,\mathrm {in}}$ and $\phi _{i,\mathrm {out}}$ are found from the complex fields $u_A$ and $u_B$ in Step 1. Since the fields $u_A$ and $u_B$ are known a priori from the field evolution, the outputs of phase retrieval (${u_A}''$ and ${u_B}'$ for Step 1) contain phase modifications that may be found by subtracting the phases of fields at the same $z$. For instance, $\phi _{i,\mathrm {out}} = \arg {u_B} - \arg {{u_B}'}$ and $\phi _{i+1,\mathrm {in}} = \arg {{u_B}''} - \arg {u_B}$. The net phase modification $\theta$ is then

Notably, the net phase modification $\theta$ is independent of the phase of the field $u_B$ in the evolution. Therefore the phase profiles in the field evolution do not directly influence the refractive index profile and only serve as an initial guess for phase retrieval. In practice we find that reasonable estimates of the phase in the field evolution greatly improve reliability and computation speed. However the approximations used in constructing the wavefront do not carry directly into the refractive index result.

In this section we have described the process of designing a 3D refractive index function to reshape a coherent field. A mapping function is continuously deformed with axial position to specify the complex field throughout the GRIN volume. Then the refractive index profile compatible with the field and the split-step BPM is recovered using phase retrieval. In the next section we implement the method and demonstrate its capabilities in an example.

3. Results

In this section we show an example GRIN design that transforms a beam with some arbitrary irradiance features. We implement a numerical algorithm that finds a mapping function between any two smooth irradiance distributions with no zeros in a closed region. We first describe the procedure and then show the numerical results. At the end of the section we numerically test the GRIN designs for accuracy and estimate their alignment tolerance.

3.1 Example GRIN design

The example GRIN design converts a Gaussian beam to a rectangular-shaped supergaussian with superimposed letters “GRIN.” The desired irradiance evolution is shown in Fig. 3. The input beam can be seen on the upper left; its irradiance has the form $I_0(x,y) = \exp {[-2(x^2+y^2)/w^2]}$, with $w = 80$ µm. The plots proceeding left to right show that the letters and supergaussian shape appear gradually and simultaneously with propagation. The output irradiance has the fully-formed letters “GRIN” summed with a rectangular beam given by $I_1(x,y) = \exp {[-2(x/w_x)^{m_x} - 2(y/w_y)^{m_y}]}$ where the beam widths are $w_x = 100$ µm, $w_y = 50$ µm and the power terms are $m_x = 10$, $m_y = 8$.

 

Fig. 3. Specified irradiance evolution for Gaussian to “GRIN” letters converter.

Download Full Size | PDF

This demonstration shows a beam transformation for an arbitrary smooth irradiance function. To form the letters “GRIN” we use a simple numerical method to obtain a sampled mapping function. However the method operates on a square spatial region and is not meant to handle irradiance zeros such as those outside the input Gaussian or the output supergaussian. To overcome this limitation, we write the overall mapping as a chain of individual mapping functions as described in Sect. 2.2.2.

The three mapping functions in the overall transformation are shown in Fig. 4. For each map, the irradiance profiles are shown in the left and middle plots, and a coordinate mesh is shown on the right. The mesh shows the initial positions $\mathbf {x_0}$ that yield a regular grid in $\mathbf {x_1}$. This represents the reference positions $\mathbf {x_0}$ that must be computed at each $\zeta$ so that $\boldsymbol {\xi }$ can be evaluated on-mesh in the homotopy. The first map (Fig. 4(a)) converts a Gaussian to a uniform square. The gridlines are concentrated in the center, signifying high irradiance, and widely spaced outside the Gaussian envelope. The second mapping, which produces the “GRIN” letters, is shown in Fig. 4(b). An image specifying the profile of the letters was created by convolving high-contrast block letters with a narrow gaussian smoothing function. The image was stretched horizontally and vertically to compensate for shrinkage in the next mapping. Alternatively this stretch could have been performed with another mapping function. The third mapping (Fig. 4(c)) converts a square-shaped uniform irradiance profile into a rectangular-shaped supergaussian. The actual input irradiance will contain letters, since it follows the second mapping. These letters would distort if they were located near the edge of the beam, but the mapping is predominantly a uniform vertical and horizontal stretch over most of the field, as is evident from the regular spacing of the gridlines. In this example the distortion of the letters is negligible.

 

Fig. 4. Mapping functions used in GRIN design. Irradiance profiles before the mapping (left column) and after the mapping (center column). Pre-distorted mesh grid that maps to a uniform mesh (right column). (a) Gaussian to uniform irradiance, (b) Uniform to “GRIN” letters with uniform background, (c) Uniform to rectangular supergaussian.

Download Full Size | PDF

The mapping that forms the “GRIN” letters (Fig. 4(b)) was found using a numerical method by Sulman [23] that solves an optimal transport problem. The goal in optimal transport is to minimize a cost function that represents the effort needed to redistribute the “mass” from a given density function in position space. The approach has been used in the design of optical freeform surfaces that convert a given irradiance profile into another. For example, freeform lens pairs have been designed to shape a coherent beam to an arbitrary irradiance pattern with a flat phase profile [24]. Similarly, ray mappings have been used to design double freeform surfaces to yield arbitrary irradiance and phase profiles [25,26] by employing Sulman’s algorithm. Though other methods have been demonstrated to solve similar optimal transport problems [27,28], we chose Sulman’s algorithm for the present work for its relative simplicity. However, Sulman’s algorithm does have some practical limitations. Firstly, it is meant to operate on a square-shaped domain. This precludes defining the support boundaries around important features in the irradiance. For example, a square domain around a circular Gaussian will unavoidably have irradiance values near zero in its corners. Secondly, the algorithm is highly sensitive to zeroes and sharp gradients in the density functions. These features tend to cause the algorithm to fail to converge. We observe that these features create large grid displacements during iteration, which can cause the newly mapped coordinate grid to become self-intersecting and therefore unusable. The tendency to create large displacements can be reduced somewhat by decreasing the iteration step size, actively detecting self-intersections, and/or adjusting the sampling. However, in practice the problems are rather persistent and the algorithm overall has a low tolerance for these features. The Gaussian and supergaussian irradiance profiles cannot be handled by Sulman’s method directly because they contain near-zero values on a square-shaped domain. But they can be joined with results from Sulman by converting them to a uniform square shaped irradiance (Figs. 4(a), 4(c)). These two mapping functions were obtained analytically. Since both mappings are separable, the 1D functions in $x$ and $y$ can be obtained by setting equal the 1D irradiance integrals,

The refractive index profile found by phase retrieval is shown in Fig. 5. The “GRIN” letter features are most noticeable because of their high contrast. They resemble the irradiance profile because in absence of strong phase gradients, light tends to be guided in regions of high refractive index. For this reason, the higher irradiance regions tend to have higher refractive index and be bounded by regions of lower refractive index. Wider gradients that perform the Gaussian-to-rectangle conversion can also be seen. Near the front of the device ($z<5$ mm), the refractive index increases with $|x|$ and decreases with $|y|$. This imparts a phase profile that encourages the beam to stretch outward in $x$ and shrink inward in $y$. Near the end of the device ($z>5$ mm), these gradients change sign to compensate the accrued phase profile and halt the transformation. The total refractive index range across the entire GRIN volume is $\Delta n = 7\times 10^{-4}$. Based on previous experiments, this index range is suitable for femtosecond laser written fabrication in Borofloat glass [29]. However the total index range can be tuned significantly by changing the length of the device. Generally, longer devices require a smaller refractive index range because they transform the beam more gradually. The length of this particular design was chosen to reach the specified index range.

 

Fig. 5. 3D refractive index profile for Gaussian to “GRIN” letters converter. For an animated version, see Visualization 1.

Download Full Size | PDF

Several physical and numerical parameters are used in the design. For the examples shown in this paper, the wavelength in vacuum is $\lambda _{\mathrm {vac}} = 520$ nm and the background refractive index is $n_0 = 1.4746$ for Borofloat glass. The lateral grid is $128\times 128$ samples and the axial step length is $\Delta z = 50$ µm. The phase retrieval algorithm is a hybrid input-output method that alternates with error-reduction every 20 iterations [30]. The refractive index $\Delta n(x,y,z)$ is assumed to be zero outside the region at each $z$ containing the highest irradiance values except the last 0.01% of the total power. 2D phase unwrapping was performed at each axial step using an available Matlab script [31,32].

3.2 Efficiency and alignment tolerance

In this section we consider the accuracy of the realized transformations and the alignment tolerance for the example GRIN design. We show the calculated output field when the input beam is displaced from center by a small distance, or tilted off-axis, or magnified in $x,y$. It is useful to understand the tolerances when testing new GRIN devices experimentally since it can be difficult to distinguish between beam distortions from misalignment versus those arising from unwanted refractive index features.

Output beams for this design are shown in Fig. 6. The top row (Fig. 6(a)) shows the calculated output with ideal alignment (lateral shift $x_s = 0$, angle tilt $\theta _x = 0$, magnification = 1). The pseudo-color plots show the 2D irradiance profiles. To the right of these are cuts of the irradiance and phase along $y=0$. The output fields were obtained by substituting perturbed fields at $z=0$ and then propagating them through the GRIN design using the split-step method. In the cut-plots, the realized field is marked by circles with dotted lines, whereas the field specification is marked by plain solid lines. Small irradiance errors are seen in the output for an ideal input beam, but no ripples, speckle, or phase abnormalities are apparent. Mode-matching efficiency, assessed by the normalized superposition

 

Fig. 6. Calculated output beam for both GRIN designs. (a) Input beam matches specification with no misalignment. (b) Input beam shifted toward $+x$ by distance $x_s$. (c) Input beam tilted toward $+x$ by angle $\theta _x$. (d) Input beam magnified in $x,y$. Cross sections along $y=0$ show ideal (solid line) and actual (circles with dotted line) irradiance and phase.

Download Full Size | PDF

Figure 6(b) shows the output field when the input beam is laterally shifted by distance $x_s = 20$ µm. In response to the lateral shift of the input beam, the output irradiance skews rightward, in the direction of the shift. The mode matching efficiency is $\eta = 0.965$. For reference, this may be compared to the matching efficiency of the ideal output beam with a shifted version of itself. Shifting by the output by $x_s$ yields $\eta = 0.91$. Perturbation values for the shift, as well as tilt and magnification below, were chosen empirically to represent a minimum recognizable distortion in the output beam. Hence, smaller perturbations produce similar-looking distortions, but to a lesser extent.

To assess tilt misalignment, the shift is reset to zero and the beam is tilted by a small angle $\theta _x = 0.15^{\circ}$. This result is shown in Fig. 6(c). The left side of the beam is lost and some ringing appears on the right. The mode matching efficiency is reduced to $\eta = 0.174$, compared to a reference value of $\eta = 0.167$ obtained from applying the tilt directly to the ideal output beam. The efficiency in both cases is low because the steep linear tilt of the phase profile ($\pi$ rad per $100$ µm) greatly diminishes the superposition.

To assess sensitivity to magnification, the tilt is reset to zero and the input beam is magnified by $\times 1.3$ in both lateral directions. The result is shown in Fig. 6(d). In this case, two irradiance peaks appear at the left and right edges of the beam, and some ringing appears in the center at spatial frequencies equal to roughly the stroke width of the “GRIN” letters. The mode matching efficiency is $\eta = 0.900$.

4. Discussion

The above example shows that a semi-arbitrary 2D field can be synthesized by a 3D GRIN device with a high degree of accuracy and efficiency. However, because we have incorporated certain properties of the input and output fields into the design process, it is worth making clear what constraints and assumptions the design approach must follow. In particular, the irradiance profiles of the input and output must be smooth and continuous, and if they contain zeroes, these must removable by manipulating the geometry. In the above example, nested mapping functions were used for this purpose. It is also possible to accommodate irradiance nulls by specifying the field evolution directly [18]. The requirements that the irradiance be smooth and contain no zeroes ensure that the Jacobian is continuous and lacks singularities. Inclusion of Sulman’s method adds further sensitivity to irradiance gradients and zeroes, but these are not much different from those imposed by the Jacobian. Anecdotally, Sulman’s method seems to impose a maximum contrast $I_{\mathrm {max}}/I_{\mathrm {min}} \approx 5$. We observed this empirically, and the same apparent limit has been reported at least once by others [33]. Also, Sulman’s method operates only on a square domain. Though it was possible to “pre-map” to such a domain in this example because of the rectangular geometry of the beam, in other cases the required pre-mapping may be more difficult.

An additional constraint on the input and output fields is that their phase profiles are assumed to be flat. This is expressed in the boundary conditions of Eqs. (7) and (13). In 2D GRIN design this was used to ensure that there were no caustics in the GRIN medium [17]. In 3D, it is a convenient simplification and should not be considered the most general form.

Let’s also consider the implications of the paraxial assumption. Paraxial waves were assumed in arriving at Eqs. (22) and (23). However, as mentioned above, Eq. (22) is used for constructing wavefronts that serve as an initial guess for phase retrieval. Upon iteration, the initial guess is lost and the converged solution will conform to the approximations of the propagation method. In the split-step BPM, we are using a full-wave scalar transfer function $H$ that does not include a paraxial approximation, but the split-step model represented by Eq. (23) does contain this approximation. Therefore, we can expect that the refractive index profile will contain errors from the split-step model. However, the errors will diminish with reduced step size $\Delta z$ to some extent [21].

The Jacobian method detailed in this paper offers several advantages over the previous GRIN design method [17]. Firstly, the Jacobian method can address general design problems where the geometry of the field cannot be simplified using independent 1D functions. The design problem with 1D fields is more straightforward because the mapping can be cast in terms of one or more 1D integrals over the irradiance, following Eq. (26). When this is the case, for instance, where the irradiance is a separable function or contains symmetry, the 2D GRIN method is suitable. In fact, we exploited separability of the Gaussian and supergaussian functions in the above example to obtain the mappings to and from a uniform square-shaped irradiance profile for Sulman (Figs. 4(a) and 4(c)). However, the Jacobian method is required for the mapping to obtain the letters (Fig. 4(b)) because those features have arbitrary shapes that cannot be simplified. Secondly, the homotopy has been reformulated as a continuous deformation of a mapping, which allows it to be computed easily at on-mesh positions in intermediate planes. In the 2D method, the homotopy was not written continuously but was uncovered discretely by dense ray sampling. This dense sampling proves to be very computationally expensive in 3D. Further, we observe that discrete sampling produces noise in the intermediate fields which interacts badly with phase retrieval. Interpolating a sampled mapping function as returned by Sulman proved to be less noisy, and optionally a denser grid can be used for Sulman’s method than for the BPM. A final difference is that, as a consequence of writing the homotopy continuously, the ray-based model from the old design method can now be seen as a sampled version of the continuous homotopy. Because the ray model is one discrete representation of the ground truth, we see that it is no longer conceptually necessary in describing the gradual implementation of the mapping(s) with $z$.

5. Conclusion

We have demonstrated a numerical design method for 3D GRIN elements that gradually change the transverse profile of a beam during propagation. The beam transformation can be prescribed by writing a mapping function to relate the input and output irradiances and then parameterizing the map as a continuous function of axial position. Once the evolution of the complex field has been specified throughout the GRIN volume, a refractive index function that implements the specification is found by phase retrieval. The resulting index profile is consistent with the paraxial scalar wave equation and diffraction effects are considered and controlled in the design. In an example GRIN design, a numerical method was utilized to obtain a mapping function for irradiance profiles with arbitrary features but no zeros or sharp gradients. Constraints to operate on a square domain were met by using multiple mapping functions to perform an additional transformation. The refractive index range of the design is suitable for fabrication by femtosecond laser writing in glass.