## Abstract

Gradient-index (GRIN) media offer advantages over thin optical elements for beam shaping of strongly diffracting fields. A numerical GRIN design method is presented, where diffraction effects are considered in solving for the refractive index profile. The index profile is found by specifying a desired beam transformation throughout the GRIN volume and solving a series of phase retrieval problems. A Gaussian to flat-top beam shaper and a coherent beam combiner are shown as examples. Reduced beam distortion is demonstrated in comparison to a purely geometric design method.

© 2016 Optical Society of America

## 1. Introduction

In the geometric approach to beam shaping with phase-only elements, a given input field is converted into a desired one using two optical elements spaced by a finite axial distance. The initial field is modeled as a collection of rays, and the first phase element tilts each ray appropriately so that propagation to the second element results in the desired intensity profile. The ray directions at the second plane are then corrected by the second optical phase element, which produces the desired phase profile on the beam. This arrangement has been implemented with several types of optical elements, including aspheric lenses [1], selectively aberrated spherical lens systems [2], and gradient-index (GRIN) lenses [3].

Although this arrangement is appealing due to its simplicity, it cannot achieve perfect beam shaping in general due to diffraction. Diffractive effects in this system are conveniently described by the parameter [4,5]

where*r*and

_{i}*r*are the half-widths of the beam at the input and output plane,

_{o}*λ*is the wavelength in vacuum,

*n*

_{0}is the background refractive index,

*L*is the propagation length, and

*C*is a constant that depends on the particular beam geometries. If the phase elements are designed using a geometric model, the system will provide a reasonably good beam transformation only if

*β*is large. For the conversion of a Gaussian beam to a flat-top, the required value of

*β*is typically

*β*> 32 [5], where

*C*= 2(2π)

^{1/2}.

Beam shapers that account for diffraction incorporate the wave nature of light in their design. Again, the light is modified at two distinct planes by two phase elements. The first imparts a phase profile that produces the desired irradiance upon propagation to the second plane [6]. As before, the second element is used to correct the leftover phase distortion. Perfect conversion is not possible in general, however, and performance degrades as *β* becomes small [4]. This leaves a need for a beam shaping method that remains accurate and efficient at all size scales.

Recently, we introduced a different paradigm of beam shaping whereby a continuous GRIN medium was employed to slowly convert a laser beam from one irradiance distribution to another [7]. Instead of imparting the necessary phase in a single plane, the inhomogeneous medium reshaped the beam gradually during propagation. The design methodology utilized a ray-based approach, limiting the technique to large beam sizes that propagate over short distances. Nevertheless, this technique was shown to produce beam shaping with theoretical efficiencies of 100%. The design is suited to several new fabrication techniques for arbitrary GRIN structures, including lithographic direct-write techniques [8,9], photopolymer lithography methods [10,11], and 3-D printing of silica powder beds [12].

In the current work we extend the design methodology of the previous GRIN beam shaper to include the effects of diffraction that are encountered with smaller input beams, tighter spot focus, and longer propagation lengths. Phase retrieval methods are used to calculate a modified inhomogeneous index profile that compensates for these diffractive effects. The advantages of 100% conversion efficiency and simplified phase correction are maintained, making the optical elements ideal for applications in the areas of semiconductor laser arrays and integrated waveguides.

## 2. Phase retrieval design method

To create a design method compatible with diffraction, we begin by choosing a suitable model for field propagation in inhomogeneous media. We use the split-step beam propagation method (BPM) [13], where the medium is divided into a series of homogeneous propagation steps of length ∆*z* with an infinitesimally thin phase element between each step. 2-D notation is used for simplicity; extensions to 3-D are straightforward. Defining *x* as the transverse position and *z _{i}* as the axial position of the

*i*

^{th}phase element, the phase

*φ*(

_{i}*x*,

*z*) acquired on the

_{i}*i*

^{th}propagation step due to the index profile

*n*(

_{i}*x*,

*z*) is applied by multiplying the field by the term exp[

_{i}*jφ*(

_{i}*x*,

*z*) ]. According to scalar paraxial BPM theory (e.g [14].) the phase is given by

_{i}*n*(

_{i}*x*,

*z*) isIn these equations,

_{i}*k*

_{0}= 2

*πn*

_{0}/

*λ*, and

*n*

_{0}can be taken as the on-axis refractive index value.

By virtue of Eqs. (2) and (3), the index profiles *n _{i}* (

*x*,

*z*) at each axial position are known if the phase profiles

_{i}*φ*(

_{i}*x*,

*z*) are known. We can use this fact to create an index profile for the desired GRIN element. The phases

_{i}*φ*(

_{i}*x*,

*z*) are not known initially, but they can be found via phase retrieval if field magnitude information is provided at the input and output of each propagation step. We will see shortly that magnitude profiles can be found by providing a desired field evolution throughout the GRIN region. This field evolution will then define the refractive index profile found by the method described presently.

_{i}To pose each propagation step as a phase retrieval problem, we divide each step into a three-part sequence composed of a homogeneous propagation length sandwiched by two phase elements *φ _{A}* (

*x*) and

*φ*(

_{B}*x*) to be determined. A single sequence is illustrated in Fig. 1(a). Assume for the moment that the complex fields just outside both phase elements,

*u*(

_{A}*x*) and

*u*(

_{B}*x*), are known in both magnitude and phase. Then the magnitudes of the fields inside the phase elements,

*u*(

_{A}'*x*) and

*u*(

_{B}'*x*), are known, but the phases are not. Dropping the understood

*x*-dependence, this can be written: |

*u*| = |

_{A}*u*| and |

_{A}'*u*| = |

_{B}*u*|, but since

_{B}'*φ*and

_{A}*φ*are unknown, ∠

_{B}*u*and ∠

_{A}'*u*are also unknown. Given the known magnitudes, however, it is straightforward to find the unknown phases using retrieval methods like the Gerchberg-Saxton algorithm [15] or input-output methods [16].

_{B}'Once the phase profiles of *u _{A}'* and

*u*are known, the phase elements

_{B}'*φ*and

_{A}*φ*are given simply by the phase difference of the fields on either side of the elements. Recalling our assumption that

_{B}*u*and

_{A}*u*are fully known, we solve the transmission expressions

_{B}*u*=

_{A}'*u*exp(

_{A}*jφ*) and

_{A}*u*=

_{B}*u*exp(

_{B}'*jφ*) for

_{B}*φ*and

_{A}*φ*. The phase functions of the neighboring propagation sequence are found in the same way, with new input and output fields. After every sequence is solved, abutting pairs of phase elements can be summed directly since their axial separation is zero. Referring to Fig. 1(b), composite phase elements are formed by

_{B}*φ*+

_{n−1,B}*φ*and

_{n,A}*φ*+

_{n,B}*φ*. The index profiles

_{n + 1,A}*n*(

_{i}*x*,

*z*) at each step are then obtained from the summed phase profiles using Eqs. (2) and (3).

_{i}## 3. Specifying the desired field evolution

The above method is useful in finding an index profile if the intermediate fields *u _{A}* and

*u*are known in advance at every propagation step. In practice the complete evolution of the field is not known and must be set by the designer. The complete process of obtaining the desired index profile then consists of two parts:

_{B}- I. Specify a field evolution that describes the gradual transformation of the field through the GRIN volume.
- II. Compute the index profile that will realize the specified field evolution.

The algorithm described in Sect. 2 addresses Part II of the design process, that is, it computes an index profile given a desired field evolution. Part I is independently accomplished by adopting a ray-model heuristic.

In this approach to Part I, the field is expressed by a set of rays with cubic polynomial trajectories through the GRIN volume [7]. Boundary conditions to determine the polynomial coefficients are given by the position and slope of rays at the input and output planes. These are found by converting the desired fields from a complex scalar representation to a ray representation. By assuming that none of the rays intersect, there is just one possible mapping between the input and output rays that follows this form. Finally, with ray trajectories filling the GRIN volume, the field can be found at each axial position.

Using a ray model to specify the field evolution does not diminish the accuracy of the GRIN design method for diffracting beams. Although wave effects normally cause a propagating field to deviate from a ray description, in this situation the error is compensated by the index profile. The phase-retrieval method ensures that the realized evolution in the GRIN element will match specification, and thus the ray-based field evolution will result in phase elements that cancel the wave effects normally responsible for the disagreement.

A few numerical details must be considered if the ray mapping is used for the field evolution. In particular, one must convert the field at each axial step from a ray-based representation to a complex-valued scalar representation. Both the magnitude and phase are needed for a complete description of the optical field.

The magnitude conversion is performed as described in [7], where the local field irradiance is found from the ray density. The phase conversion is performed as shown in Fig. 2, where the phase profile is inferred from the ray slopes. At a plane of interest (*z* = *z _{i}*) where the phase profile is desired, the ray slopes

*θ*

_{1},

*θ*

_{2}, … are known from the cubic ray paths. These are incident from the left in Fig. 2(a). The wavefront is constructed from the planes orthogonal to each ray, and phase values are then obtained from the

*z*height of the constructed wavefront as shown in Fig. 2(b). Since the ray slope is given by

*dz*/

*dx*= tan

*θ*, and for a paraxial wave

*dφ*=

*k*

_{0}

*dz*, the interpolated phase can be written as

*x*. Once the magnitude and phase profiles have been found from the ray mapping, the field is defined at every step position. This completes the specification of the field evolution; the index profile is then found by phase retrieval.

## 4. Example applications of the GRIN design method

To demonstrate the method and to assess its accuracy, we consider two example beam-shaping problems. The first is a Gaussian to flat-top mode conversion. We solve for the index profile using phase retrieval as described in Sect. 2 and then compare the solution to one obtained from a purely geometric design method. The second example is a GRIN element that performs coherent addition of multiple beams. The GRIN medium performs complete aperture-filling by preventing interference between neighboring beams.

#### 4.1 Gaussian to flat-top beam shaper

In our first example, we apply the algorithm to the design of an index profile that converts a Gaussian beam to a flat-top. The complex fields at the input and the output have the general form

where*a*is the half-width of the beam and

*m*is the supergaussian order. A Gaussian at the input is obtained by setting

*m*= 2, and a supergaussian with

*m*= 10 is used at the output. The phase profiles of both fields are flat. We choose

*a*= 50 μm for both fields, a GRIN length of 5 mm, a wavelength

*λ*= 1 µm, and background index

*n*

_{0}= 1.5. The step size ∆

*z*is 5 μm.

With these specifications, we create a cubic ray mapping and then run the phase retrieval algorithm to find the index profile that performs the desired conversion. An example field evolution is shown in Fig. 3, where we have used different beam sizes at input and output to emphasize the ray bending. The local irradiance is given by the ray spacing; the input Gaussian at *z* = 0 is represented by rays that are placed farther apart upon increasing distance from the axis, and the ray spacing is largely uniform for the output beam at *z* = 5 mm. These ray trajectories imply the irradiance evolution shown in the figure (red), where it is clear that the beam changes continuously from Gaussian to flat-top over the full length of the converter. The phase profiles at each plane (not shown) are recovered from the ray slopes. The mapped field is then used as an input to the phase retrieval algorithm.

By constructing the desired field evolution and then running the method of Sect. 2, we obtain the index profile shown in Fig. 4(a). We can understand on a basic level how the index profile performs the desired conversion by considering the concavity of *n*(*x*) at a fixed *z*. Near the input plane (*z* = 0), the index profile acts as a negative lens near the axis and as a positive lens farther out. The combined effect is to pull beam energy outward from the center of the Gaussian and inward from the tails. This continues as the beam propagates along *z* until about the middle of the converter, at which point the transverse index gradient changes direction and the maxima on either side of the axis are replaced by minima. This slows the beam expansion and compensates the transverse phase profile that has built up during propagation in the first half of the GRIN element. By the very end of the GRIN element, the beam is a collimated supergaussian and the index has developed positive lobes at the edges of the beam. This resembles the step-index structure that is used in waveguides for flat-top modes [17]. We verify the solution by propagating the input field through the GRIN element using the split-step method. The propagated field is compared to the output specification in Fig. 4(b), and good agreement is clear.

We illustrate the advantage of the phase-retrieval design method by comparing this beam shaper to one that was designed using a geometric method [7]. We use the same field evolution, but the resulting index profile will neglect diffraction since a geometric-optics propagation model is assumed. The index profile is shown in Fig. 4(c). Its general shape is similar to the diffraction-compensated index profile of Fig. 4(a), but differences are apparent particularly at large *z*. Most notably, the direction of the index gradient is positive with increasing distance from the axis. This anti-waveguiding behavior clearly does not compensate diffraction and will cause the beam to distort with further propagation.

The shortcomings of the geometric index profile become further evident when the input Gaussian beam is propagated through it. The output intensity for the geometric GRIN element is compared to specification in Fig. 4(d). The lack of guiding index structure causes the beam to spread while it passes through the GRIN element, resulting in an incomplete mode conversion. Both the geometric and the phase retrieval methods produce an index profile that causes the field to follow the specified evolution according to the respective propagation model. However, the geometric model does not accurately describe the wave-like character of the beam, so the actual field does not obey the specified evolution if it is strongly diffracting. As a quick gauge, this example has a Fresnel number of 0.75, so diffraction effects are expected to be noticeable. Given this, it is not surprising that the index profile of the phase retrieval method differs from that of the geometric method.

#### 4.2 Beam addition for three mutually coherent Gaussian beams

As a second example, we consider a GRIN element designed for the coherent combination of multiple beams. GRIN media are advantageous for this since typical beam combiners in free-space utilize lossy, alignment-sensitive elements like gratings, spatial filters, and beam splitters whereas GRIN media are low-loss and low-complexity. In integrated optics, the step-index Y-branch waveguide is common. This structure is efficient but necessarily single-mode. Its small size puts it at risk of optical damage or nonlinear effects at high power, whereas a GRIN element can be operated at larger size scales.

In this demonstration we consider a linear array of three identical Gaussian beams which we wish to convert into a single Gaussian. This is accomplished using two sequential GRIN stages. The first stage performs an aperture-filling operation, and the second stage reshapes and resizes the beam. The design of such a GRIN element is shown in Fig. 5. Each of the input beams has a spot radius of 50 µm and the beam centers are separated by 125 µm. At this scale the beams will be both strongly diffracting and strongly interacting. The beams are assumed to be mutually coherent and in-phase. As with the mode converter, the desired field evolution is described by a collection of rays with cubic polynomial trajectories. Two separate sets of polynomials are used for the two stages of the GRIN element. The rays are shown with the desired intensity evolution in Fig. 5(a). Cross-sections of the desired intensity at the input, middle, and output planes are shown in Fig. 5(b). A supergaussian with *m* = 10 is used at the midplane to specify the aperture filling.

The index profile is again found using the phase retrieval approach, and the three Gaussians are propagated through the GRIN element afterwards to check the performance. The index solution is shown in Fig. 5(c). The two stages of the index profile each have a distinct character: the first stage pulls the three Gaussians together and then corrects the accrued phase, while the second focuses the beam and follows with recollimation. The intensity after propagation is shown in Fig. 5(d), which again shows good agreement. The *M*^{2} value and Strehl ratio are both found to be within two decimal places of unity.

Although beam resizing is commonly needed in application, the aperture-filling operation is generally more important in coherent beam addition [18]. When performing aperture-filling with thin diffractive optics, there is usually a tradeoff between fill factor and beam interaction. As the beams are made to overlap more to increase fill factor, their mismatched wavefronts cause undesired interference. This is one area in which a GRIN element provides a distinct advantage: guiding effects provided by the continuous medium ensure that the wavefronts match one another so there is no interference and the beams behave as a single entity throughout.

## 5. Discussion

Whether fabrication is feasible for the above examples can be checked by calculating the maximum contrast and gradient of the refractive index. Both numbers depend on the width of the GRIN structure since the gradient becomes steep and monotonic at large |*x*|. A practical width must be chosen by the designer. We set the width of the GRIN channel to vary with *z* such that 0.02% of the beam energy lies outside it. This results in virtually no distortion of the propagated beam. For the mode converter of Fig. 4(a), this results in a maximum index contrast of Δ*n* ~1 × 10^{−3} and a gradient of *dn*/*dx* ~0.4 mm^{−1}. The beam-combining structure of Fig. 5(c) requires greater index contrast Δ*n* ~2.5 × 10^{−2}, and a similar gradient *dn*/*dx* ~0.2 mm^{−1}. Based on these values, fabrication of the first device should be possible with e.g. either lithographic photopolymer methods [10] or direct-write [8], while the second device is compatible only with the former method due to the need for high index contrast in the second half of the GRIN profile.

The above numerical values apply to the particular examples shown here, but in general the contrast and gradient of the refractive index are sensitive to the size scale. In particular, the necessary index contrast shrinks as the length of the GRIN element increases. This is expected since the intensity remapping will occur more gradually, and the necessary phase to achieve the transformation is applied over a greater length. In principle, then, the length of the GRIN element can be chosen to suit the given fabrication capabilities. This design flexibility emphasizes the advantage of the phase retrieval design method over the geometric method. As the device length increases, diffraction becomes more prominent and the geometric method becomes inaccurate, while the phase retrieval method continues to yield the desired beam shape.

The examples have shown that the phase retrieval method performs well in some cases that are challenging in a geometric approach, but the method is not perfectly general. To start, the convenient cubic ray technique employed in Part I of the algorithm can put limitations on the phase profiles of the input and output fields. This is because ray intersections must be avoided inside the GRIN element. Intersections imply infinite local irradiance accompanied by sudden jumps in phase, both of which are unphysical. Whether ray intersections actually occur for a particular pairing of fields depends on the phase and intensity profiles. However, it is guaranteed that no intersections occur if both phase profiles are flat. In our examples we obey the flat-phase condition, but in fact the use of non-flat phase profiles is acceptable in many cases.

Part II of the algorithm (finding the index from a given field evolution) applies quite generally, but we have made paraxial approximations in our implementation. The division of the inhomogeneous medium into discrete phase elements continues to be valid at smaller size scales and broader angular bandwidths. The BPM has been cast at many levels of approximation [19], including one that uses a fully vectorial Helmholtz equation. In this work we have used the paraxial BPM to simplify our demonstrations, but better accuracy can be achieved under different conditions by switching to a more rigorous BPM.

Fortunately, changes made in one part of the approach have little bearing on the other. If a better strategy for defining the field evolution (Part I) becomes available, it can be substituted without necessitating changes to the phase retrieval method (Part II). Likewise, if a different propagation method is needed to improve accuracy at smaller size scales, the field evolution need not change unless it contains overlapping assumptions. Equation (4), for instance, uses a slowly-varying envelope approximation. If a more precise propagation method is to be used, then Eq. (4) should be considered invalid and the wavefront should be constructed with a more general technique, as in [7].

## 6. Conclusion

We have demonstrated a GRIN design method for coherent beam shapers where diffraction is important. The GRIN volume is represented by a series of axially spaced phase elements which are found by phase retrieval given an arbitrarily-chosen field evolution. To produce the field evolution we have used a particularly simple ray mapping scheme. The method was demonstrated by designing a GRIN beam shaper that performs a Gaussian to flat-top conversion. The resulting element was shown to perform better than one designed using a purely geometric approach. Coherent addition of three Gaussian beams was also demonstrated in a separate GRIN element.

GRIN beam shapers enjoy a significant advantage over thin-element systems. By compensating beam diffraction in the index profile, the GRIN element behaves simultaneously as beam shaper and waveguide. The guiding property allows the phase and intensity changes to occur gradually, resulting in near 100% conversion efficiency at size scales that are inaccessible to other beam shaping methods.

## Acknowledgments

The authors acknowledge financial support from Defense Advanced Research Projects Agency (DARPA) (HQ0034-14-D-0001) and the United States Department of Defense Air Force (FA9550-14-1-0382). We gratefully acknowledge Di Lin for providing his geometric GRIN design code.

## References and links

**1. **P. W. Rhodes and D. L. Shealy, “Refractive optical systems for irradiance redistribution of collimated radiation: their design and analysis,” Appl. Opt. **19**(20), 3545–3553 (1980). [CrossRef] [PubMed]

**2. **B. R. Frieden, “Lossless conversion of a plane laser wave to a plane wave of uniform irradiance,” Appl. Opt. **4**(11), 1400–1403 (1965). [CrossRef]

**3. **C. Wang and D. L. Shealy, “Design of gradient-index lens systems for laser beam reshaping,” Appl. Opt. **32**(25), 4763–4769 (1993). [CrossRef] [PubMed]

**4. **L. A. Romero and F. M. Dickey, “Lossless laser beam shaping,” J. Opt. Soc. Am. A **13**(4), 751–760 (1996). [CrossRef]

**5. **F. M. Dickey and S. C. Holswade, “Beam shaping: a review,” in *Laser Beam Shaping Applications*, F. M. Dickey, S. C. Holswade, and D. L. Shealy, eds., 269–303 (CRC, 2006).

**6. **W.-H. Lee, “Method for converting a Gaussian laser beam into a uniform beam,” Opt. Commun. **36**(6), 469–471 (1981). [CrossRef]

**7. **D. Lin and J. R. Leger, “Numerical gradient-index design for coherent mode conversion,” Adv. Opt. Technol. **1**, 195–202 (2012).

**8. **R. R. Gattass and E. Mazur, “Femtosecond laser micromachining in transparent materials,” Nat. Photonics **2**(4), 219–225 (2008). [CrossRef]

**9. **A. Žukauskas, I. Matulaitienė, D. Paipulas, G. Niaura, M. Malinauskas, and R. Gadonas, “Tuning the refractive index in 3D direct laser writing lithography: towards GRIN microoptics,” Laser Photonics Rev. **9**(6), 706–712 (2015). [CrossRef]

**10. **A. C. Urness, E. D. Moore, K. K. Kamysiak, M. C. Cole, and R. R. McLeod, “Liquid deposition photolithography for submicrometer resolution three-dimensional index structuring with large throughput,” Light Sci. Appl. **2**(3), e56 (2013). [CrossRef]

**11. **A. C. Urness, K. Anderson, C. Ye, W. L. Wilson, and R. R. McLeod, “Arbitrary GRIN component fabrication in optically driven diffusive photopolymers,” Opt. Express **23**(1), 264–273 (2015). [CrossRef] [PubMed]

**12. **H. R. Wang, M. J. Cima, B. D. Kernan, and E. M. Sachs, “Alumina-doped silica gradient-index (GRIN) lenses by slurry-based three-dimensional printing (S-3DP),” J. Non-Cryst. Solids **349**, 360–367 (2004). [CrossRef]

**13. **J. A. Fleck Jr, J. R. Morris, and M. D. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. (Berl.) **10**(2), 129–160 (1976). [CrossRef]

**14. **T.-C. Poon and T. Kim, *Engineering Optics with MATLAB* (World Scientific, 2006).

**15. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttg.) **35**(2), 237–246 (1972).

**16. **J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. **3**(1), 27–29 (1978). [CrossRef] [PubMed]

**17. **A. K. Ghatak, I. C. Goyal, and R. Jindal, “Design of a waveguide refractive index profile to obtain a flat modal field,” Proc. SPIE **3666**, 40–44 (1999). [CrossRef]

**18. **J. R. Leger, “External methods of phase locking and coherent beam addition of diode lasers,” in *Surface Emitting Semiconductor Lasers and Arrays*, G. A. Evans and J. M. Hammer, eds. (Academic, 1993).

**19. **R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical techniques for modeling guided-wave photonic devices,” IEEE J. Sel. Top. Quantum Electron. **6**(1), 150–162 (2000). [CrossRef]