## Abstract

We introduce a Lüneburg lens design where Kerr nonlinearity is used to compensate for the focal point shift caused by diffraction of a Gaussian source. A computationally efficient iterative method introduced in [Opt. Lett. **35**, 4148 (2010)] is used to provide ray diagrams in the nonlinear case and verify the focal shift compensation. We study the joint dependence of focal shift on waist size and intensity of Gaussian source, and show how to compensate spherical aberration caused by the nonlinearity by a small perturbation of the Lüneburg profile. Our results are specific to Lüneburg lens but our approach is applicable to more general cases of nonlinear nonperiodic metamaterials.

© 2011 Optical Society of America

## 1. Introduction

The gradient index (GRIN) element proposed by R. K. Lüneburg and bearing his name [2] is spherically symmetric and able to focus a plane wave incoming from one side of the sphere into a perfect geometric image point at the opposite edge of the sphere. It is widely used in applications such as antenna arrays [3], reflectors [4], optical cloaking [5], etc. However, if the input is a Gaussian beam, the focal point (i.e. waist of Gaussian beam) will be shifted beyond the edge of the sphere. Thus the initial desirable property of Lüneburg lens is lost with Gaussian beam illumination.

It has been shown that nonlinearity is useful in a number of applications, such as nanophotonic devices [6–8] and imaging [9]. In this paper, we propose that nonlinear Kerr effect [10–13] can be used to compensate the focal shift of Lüneburg lens due to diffraction. In order to do so, optical wave propagation in nonlinear inhomogeneous media needs to be studied. Previously, a lot of work has been done on this topic, both theoretically and experimentally, such as nonlinear light propagation in photonic lattices using optical induction [14,15], nonlinear pulse propagation in optical fibers [16], and beam propagation in nonlinear periodic potential with Wannier functions [17]. Recently, we introduced an iterative ray tracing method for nonlinear beam propagation [1] based on Hamiltonian ray tracing [18–20] and the Wigner distribution function (WDF) [21–23], which is applicable to both bulk media and metamaterials. This method is adopted here to validate the idea presented in this paper. Nonlinear beam propagation is studied in a subwavelength aperiodic dielectric Lüneburg lens structure composed of silicon rod arrays embedded in air [24], with consideration of refractive index changes of silicon rods due to Kerr effect [25]. Our results show that the focal shift compensation is successfully achieved as a result of focusing induced by Kerr effect. In addition, nonlinear Lüneburg lens is also investigated with conventional finite-difference time-domain (FDTD) method and a good agreement between the results from both methods is found. Besides, a modified Lüneburg lens is presented to minimize the spherical abberation caused by the introduction of nonlinearity.

An additional benefit of introducing nonlinearity into design of nanophotonic devices is that it provides more design flexibility. As an example, finite conjugate imaging, i.e. the object is at finite distance rather than infinity, is demonstrated by tuning nonlinearity in the aperiodic subwavelength nanostructured Lüneburg lens.

More importantly, though a lot of work has been done on metamaterials [26, 27], nonlinear metamaterials have not been investigated thoroughly yet. Thus we seek to fill this gap by taking nonlinear subwavelength aperiodic Lüneburg lens as an example.

In Section 2, we describe the simulation setup and the design of the subwavelength aperiodic Lüneburg lens. In Section 3, the nonlinear iterative method used for investigating nonlinear Lüneburg lens is explained. In Section 4, we show the simulation results demonstrating the compensation of focal shift by nonlinearity. In addition, focal shift under different types of sources and optical intensities is further discussed. Finally in Section 5, a modified nonlinear Lüneburg lens is designed to minimize the spherical aberration caused by Kerr effect.

## 2. Setup

Lüneburg lens is a spherically symmetric lens with index distribution
$n\left(r\right)={n}_{0}\sqrt{2-{\left(r/R\right)}^{2}}$, where *n*_{0} is the ambient index outside the lens, *R* is the radius of the lens, and *r* is the distance to the lens center. To investigate the focal shift due to diffraction and focal shift compensation by the Kerr effect, a Lüneburg lens illuminated by a Gaussian beam is studied (Fig. 1). The radius of the lens is *R*, and distances between the beam waist of input/output and the left/right edge of the lens are *z* and *z*′, respectively.

The exact Lüneburg lens profile is not achievable by using bulk media. Therefore, in this paper, we study a two-dimensional (2D) subwavelength aperiodic nanostructured Lüneburg lens [24] (Fig. 2). This structure consists of silicon rods (*n* = 3.46) centered at each square lattice with lattice constant *a*_{0} = *λ*/8, where *λ* = 1550 nm is the free space wavelength. The rod radii profile is
$a\left(r\right)={a}_{1}\sqrt{2-{\left(r/R\right)}^{2}}$, where *a*_{1} = 0.367*a*_{0}, *a*_{2} = −0.101*a*_{0} and *R* = 30*a*_{0}. The ambient medium is air (*n* = 1). Since *λ* is much larger than the lattice constant, the nanostructure could be treated as an effective medium. In this case, effective refractive indices are locally-modulated by controlling the radii of the rods, mimicking the Lüneburg lens index profile. Effective refractive index distribution of this lens is
${n}_{\mathit{eff}}={n}_{\mathit{eff}1}\times \sqrt{2-{\left(r/R\right)}^{2}}+{n}_{\mathit{eff}2}$, where *n _{eff}*

_{1}= 1.9239 and

*n*

_{eff}_{2}= −0.0678. Throughout this paper, our discussions are based on this 2D subwavelength aperiodic nanostructure example.

## 3. Method

To validate the focal shift compensation, both the nonlinear iterative method [1] and FDTD are used. More emphasis is put on the former for the reason that it is not only computationally efficient but also physically more intuitive, providing a ray diagram for nonlinear beam propagation. In this section, nonlinear iterative method is briefly described for convenience. For more details, the reader may consult [1].

Figure 3 presents the block diagram of the iterative method. Each iteration consists of 3 steps to generate a new estimate of the intensity profile (i.e. refractive index distribution) based on the previous iteration. Final beam propagation result is obtained when intensity profile estimates converge in the end.

In the first step, initial condition (position and momentum) of all the rays are defined. To take the wave effects into account, multiple rays with different momenta (directions) are emanating from each point of the source plane, as suggested by Huygens’ principle. Besides, each ray is assigned with a generalized radiance, defined by the value of WDF for given position and momentum. Generalized radiance is conserved along the ray trajectory [23]. In the second step, Hamiltonian ray tracing method is applied to generate the paths for all the rays. This method is computationally much more efficient than the FDTD method. Moreover, it works well for subwavelength nanostructure case under the ‘local periodicity’ assumption [19], while traditional ray tracing solution to this structure is unknown. For example, for the Lüneburg nanostructure investigated in this paper, we can apply the Hamiltonian method because radii of neighboring rods are almost the same thus it satisfies the ‘local periodicity’ assumption. Under this assumption, Hamiltonian ray tracing is implemented by solving the following two coupled differential equations (Eq. (1)):

*ω*is frequency, and

**x**and

**k**are position and momentum, respectively, along the path and

*σ*parameterizes the ray trajectories. Here

*ω*(

**x**,

**k**) is a function determined by the dispersion diagram of the unit cell. In the third step, a new estimate of the intensity distribution is generated from the ray trajectories as the projections of the WDF along the momentum direction. Besides, a new estimate of refractive index distribution is also calculated based on Kerr effect for use in the next iteration.

Each iteration simulates a new intensity distribution estimation. After a certain number of iterations, the intensity profile converges and this gives the final beam propagation result in the Kerr effect medium.

## 4. Simulation results

Focal shift and its compensation in the nonlinear Lüneburg lens are investigated with the iterative method. The results are are compared with FDTD results performed by the simulation software MIT Electromagnetic Equation Propagation (MEEP) [28].

As a simple example, a Gaussian beam with beam width 36*a*_{0} and peak intensity *I*_{0} is illuminating the subwavelength Lüneburg lens from the left. Without the Kerr effect, i.e. in the linear Lüneburg lens case, a focal point shift to the right can be clearly seen as shown in Figs. 4(a) and 4(b). The focal shift is due to diffraction of the input Gaussian source as compared with the ideal plane wave input case. Therefore, with a Gaussian beam input, the property that Lüneburg lens produces a geometrical focal point exactly at the opposite edge of the lens is invalid.

However, when taking Kerr effect into consideration, i.e. in the nonlinear Lüneburg case, the focal shift could be compensated and focal point is shifted back to the edge of the lens, as illustrated in Figs. 4(c) and 4(d).

In the above simulations, nonlinear Kerr coefficient for silicon is 2.7 × 10^{−14} cm^{2}/W [10]. For the iterative method simulation, grid size used is 0.03*a*_{0} × 0.03*a*_{0}. After six iterations, the beam propagation profile converges. For the FDTD simulation, grid size used is 0.083*a*_{0} × 0.083*a*_{0}.

In Fig. 4, the Kerr nonlinearity can compensate the diffraction because it increases the effective refractive index according to the intensity distribution. Effective refractive index is determined by the radius and refractive index of the dielectric rod in it. In the nonlinear Lüneburg lens case, the refractive index of the rod increases according to the local optical intensity. As a result, the effective refractive index distribution of the nonlinear case is higher than the linear case. Figure 4(e) illustrates the effective refractive index difference between the nonlinear and linear Lüneburg lens. Two peaks are observed in this figure. The peak at the right edge of the lens is due to the high optical intensity near the focal point, resulting higher refractive indices of these rods. The other peak at the center is due to the relative larger radii of these rods. Assuming similar intensity, effective indices increase more for inner part of the lens than the outer part, since the rods near the center occupy a larger fraction of space in a unit cell. Therefore, effective refractive index gradient is increased in the nonlinear case.

Figure 4 should be compared with Gutman’s “modified Lüneburg lens” [29] where a spherically symmetric perturbation to the refractive index is employed to similarly shift the focus to the left (i.e. inside the lens). The nonlinear refractive index change in our case is of the same order of magnitude with Gutman’s for the range of focal shifts considered; however, in our case the change is obviously not spherically symmetric, due to the high intensity in the focal region.

From the discussions above, the rays in the nonlinear Lüneburg lens experience more effective refractive index gradient and bend more as compared to the linear lens case. Thus, the Kerr effect is able to move the focal point towards the center of the lens, compensating the focal shift caused by diffraction. The nonlinearity may also alter the beam profile, but the change is not significant in this case, as shown in Fig. 4(f).

Furthermore, the relationship between the focal points, sources and optical intensities is investigated, as shown in Fig. 5. Three type of sources (plane wave, point source and Gaussian source) with various optical intensities are discussed. In the simulation setup, the beam waist location is kept fixed at *z* = 2*R*. For plane wave illumination, which is an extreme case of Gaussian source with infinite size of waist, the focal point is at the opposite edge when the intensity *I* approaches 0. Increasing the intensity gradually moves the focal point inside the lens, due to the resulting higher Kerr nonlinearity. For the point source case, which is another extreme case of Gaussian source where the waist approaches zero, focal point without nonlinearity is at *z*′ = 0.5*R*. This matches the prediction from geometrical imaging condition for Lüneburg lens under paraxial approximation, i.e. *zz*′ = *R*^{2}. Similar to the plane wave case, the focal point also moves towards the left when the intensity is higher. However, it will not reach the edge within the intensity range discussed here since diffraction is too large to be compensated. Between the two extremes is Gaussian incidence case. We show an example of a Gaussian beam with waist of 9*a*_{0}. The nonlinear effect exactly cancels the diffraction when *I* = 2.8*I*_{0}, i.e. the output beam focuses at *z*′ = 0. In general, the smaller the beam waist is, the higher the intensity is required for the cancellation between diffraction and nonlinearity.

## 5. Modified Lüneburg lens

As shown in the previous section, nonlinearity can be utilized to cancel diffraction and a proper designed nonlinear Lüneburg lens can focus Gaussian beam input to the opposite edge of the lens. However, spherical aberration is introduced in this treatment. For example, a Gaussian beam with waist 125*a*_{0} and peak intensity *I*_{0} is incident on the nonlinear Lüneburg lens. The beam waist is chosen so that it is much larger than the diameter of the lens. The spherical aberration is shown by the ray tracing results in Fig. 6(a). The aberration is resulting from non-uniform nonlinearity-induced refractive index gradient. As a result, rays near the center experience more nonlinearity-induced bending than the off-center ones, thus resulting in in spherical aberration.

Here, a modified nonlinear Lüneburg lens is designed to minimize the aberration. The modified structure is chosen through an optimization method. In this structure a polynomial is used to describe the radii distribution of the rods. The particular values of the coefficients in the polynomial vary with the input condition. In the example discussed above, the modified radii distribution of the rods across the diameter of the lens is shown in Fig. 6(c), along with comparison to the original profile. The radii decrease for the outer part of the lens while the inner part remains almost the same. In this way, the effective index gradient of the outer part of the lens is increased so the rays in this region bend more, resulting in reduced aberration, as shown in Fig. 6(b). In addition, spherical symmetry is preserved, ensuring the same focusing property for all directions. Furthermore, comparison of the focal point profile between the original and modified nonlinear Lüneburg lens is illustrated in Fig. 6(d). It can be seen that quality of focus is improved significantly for the modified lens.

## 6. Conclusion

In conclusion, we present a nonlinear subwavelength aperiodic Lüneburg lens which can actively control its focal point by tuning the nonlinearity. A detailed discussion on the focal shift and its relationship to different source with different optical intensity is provided. In addition, a modified Lüneburg lens is also designed to reduce the aberration caused by the Kerr effect. We fill the gap of investigating nonlinear metamaterials by taking nonlinear Lüneburg lens as an example.

## Acknowledgments

We are grateful to Baile Zhang and Se Baek Oh for useful discussions, and to one anonymous reviewer for constructive criticism. Funding support is from the following sponsors: Singapore’s National Research Foundation through the Singapore-MIT Alliance for Research and Technology (SMART) Centre and Air Force Office of Scientific Research MURI program on Nanomembranes under contract No. FA9550-08-1-0379.

## References and links

**1. **H. Gao, L. Tian, B. Zhang, and G. Barbastathis, “Iterative nonlinear beam propagation using Hamiltonian ray tracing and Wigner distribution function,” Opt. Lett. **35**, 4148–4150 (2010). [CrossRef]

**2. **R. K. Lüneburg, Mathematical Theory of Optics (Brown U.P., Providence, 1944).

**3. **H. Mosallaei and Y. Rahmat-Samii, “Nonuniform Lüneburg and two-shell lens antennas: radiation characteristics and design optimization,” IEEE Trans. Antenn. Propag. **49**, 60–69 (2001). [CrossRef]

**4. **C. S. Liang, D. A. Streater, J.-M. Jin, E. Dunn, and T. Rozendal, “A quantitative study of Lüneburg-lens reflectors,” IEEE Antennas Propag. Mag. **47**, 30–42 (2005). [CrossRef]

**5. **N. A. Mortensen, O. Sigmund, and O. Breinbjerg, “Prospects for poor-man’s cloaking with low-contrast all-dielectric optical elements,” J. Eur. Opt. Soc. Rapid Publ. **4**, 09008 (2009). [CrossRef]

**6. **M. Takahashi and H. Goto, *Progress in Nonlinear Optics Research* (Nova Science Publishers, 2008).

**7. **M. Soljačić, C. Luo, J. Joannopoulos, and S. Fan, “Nonlinear photonic crystal microdevices for optical integration,” Opt. Lett. **28**, 637–639 (2003). [CrossRef]

**8. **J. Bravo-Abad, S. Fan, S. G. Johnson, J. D. Joannopoulos, and M. Soljačić, “Modeling nonlinear optical phenomena in nanophotonics,” J. Lightwave Technol. **25**, 2539–2546 (2007). [CrossRef]

**9. **D. V. Dylov and J. W. Fleischer, “Nonlinear self-filtering of noisy images via dynamical stochastic resonance,” Nat. Photonics **4**, 323–328 (2010). [CrossRef]

**10. **R. Boyd, Nonlinear optics , (3rd ed.) (Academic, 2008).

**11. **R. Y. Chiao, E. Garmire, and C. H. Townes, “Self-Trapping of Optical Beams,” Phys. Rev. Lett. **13**, 479 (1964). [CrossRef]

**12. **P. L. Kelley, “Self-Focusing of Optical Beams,” Phys. Rev. Lett. **15**, 1005 (1965). [CrossRef]

**13. **Y. Kivshar and G. Agrawal, Optical solitons , (Academic, 2003).

**14. **J. W. Fleischer, M. Segev, N. K. Efremidis, and D. N. Christodoulides, “Observation of two-dimensional discrete solitons in optically induced nonlinear photonic lattices,” Nature **422**, 147–150 (2003). [CrossRef]

**15. **D. N. Christodoulides, F. Lederer, and Y. Silberberg, “Discretizing light behavior in linear and nonlinear waveguide lattices,” Nature **424**, 817–823 (2003). [CrossRef]

**16. **D. Anderson, “Variational approach to nonlinear pulse propagation in optical fibers,” Phys. Rev. A **27**, 3135–3145 (1983). [CrossRef]

**17. **G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, “Wannier functions analysis of the nonlinear Schrödinger equation with a periodic potential,” Phys. Rev. E **66**, 046608 (2002). [CrossRef]

**18. **K. B. Wolf, *Geometric optics on phase space*, (Springer, 2004).

**19. **Y. Jiao, S. Fan, and D. A. B. Miller, “Designing for beam propagation in periodic and nonperiodic photonic nanostructures: Extended Hamiltonian method,” Phys. Rev. E **70**, 036612 (2004). [CrossRef]

**20. **P. S. J. Russel and T. A. Birks, “Hamiltonian optics of nonuniform photonic crystals,” J. Lightwave Technol. **17**, 1982–1988 (1999). [CrossRef]

**21. **A. Walther, “Radiometry and coherence,” J. Opt. Soc. Am. **58**, 1256–1259 (1968). [CrossRef]

**22. **E. Wolf, “Coherence and radiometry,” J. Opt. Soc. Am. **68**, 6–17 (1978). [CrossRef]

**23. **M. Bastiaans, “Transport equations for the Wigner distribution function,” Opt. Acta **26**, 1265–1272 (1979). [CrossRef]

**24. **S. Takahashi, C. Chang, S. Y. Yang, and G. Barbastathis, “Design and fabrication of dielectric nanostructured Luneburg lens in optical frequencies,” in *Optical MEMS and Nanophotonics*, (IEEE Photonics Society, 2010), Paper Th1-1, pp. 177–178.

**25. **H. Gao, S. Takahashi, L. Tian, and G. Barbastathis, “Nonlinear Kerr effect aperiodic Lüneburg lens,” in *Optical MEMS and Nanophotonics*, (IEEE Photonics Society, 2010), Paper Th1-2, pp. 179–180. [CrossRef]

**26. **D. Schurig, J. Mock, B. Justice, S. A. Cummer, J. Pendry, A. Starr, and D. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science **314**, 977 (2006). [CrossRef]

**27. **J. Valentine, S. Zhang, T. Zentgraf, E. Ulin-Avila, D. A. Genov, G. Bartal, and X. Zhang, “Three-dimensional optical metamaterial with a negative refractive index,” Nature **455**, 376–379 (2008). [CrossRef]

**28. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. **181**, 687–702 (2010). [CrossRef]

**29. **A. Gutman, “Modified Lüneburg lens,” J. Appl. Phys. **25**, 855–859 (1954). [CrossRef]