Open Access
Add to BibSonomy
Abstract
Graded-index optical elements are capable of shaping light precisely and in very specific ways. While classical freeform optics uses only a two-dimensional domain such as the surface of a lens, recent technological advances in laser manufacturing offer promising prospects for the realization of arbitrary three-dimensional graded-index volumes, i.e. transparent dielectric substrates with voxel-wise modified refractive index distributions. Such elements would be able to perform complex light transformations on compact scales. Here we present an algorithmic approach for computing 3D graded-index devices, which utilizes numerical beam propagation and error reduction based on gradient descent. We present solutions for millimeter-sized elements addressing important tasks in photonics: a mode sorter, a photonic lantern and a multimode intensity beam shaper. We further discuss suitable cost functions for all designs to be used in the algorithm. The 3D graded-index designs are spatially smooth and require a relatively small refractive index range in the order of 10−2, which is within the reach of direct laser writing manufacturing processes such as two-photon polymerization.
Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. Introduction
Light has become one of the most powerful and versatile tools. It is an important medium for transporting information, either between communicating partners or from an object under investigation to a detector in imaging applications. Light also plays a key role in laser-based manufacturing processes, such as welding, cutting and direct laser lithography.
However, as recently highlighted [1–3], exploiting the full physical potential of light requires solutions for its highly specific and custom shaping. For example, increasing the telecommunication bandwidth demands directing photonic signals from many single mode fibers into a single multi-core or multimode fiber with high efficiency [4,5]. Likewise, efficiently exploiting the high optical power of industrial lasers requires reshaping the often unsuitable native mode profiles [6,7].
In any case, the restructuring of light should ideally take up only little space. To this end, the processing of dielectric materials such as glasses or polymers with ultra-fast lasers has opened promising routes towards the creation of three-dimensional, miniaturized light shapers [8]. A variety of devices such as photonic lanterns [9,10], beam combiners [11,12] and mode multiplexers [13] made from glass-embedded waveguide arrangements have been successfully demonstrated and are making their way towards commercialization. Such waveguides are formed by translating a femtosecond laser focus through a glass volume along the desired guiding paths, along which the refractive index of the glass is permanently modified. Whilst undoubtedly powerful, devices using waveguides as “building blocks” are nevertheless limited in the sense that they don’t fully exploit all degrees of freedom of a volume. Ideally, one would like to have the power of arbitrarily changing the entire three-dimensional refractive index distribution at the microscale. Progress towards the fabrication of such 3D gradient index materials for optical wavelengths has been recently reported using two-photon polymerization [14–16] and 3D printed glass optics [17]. Likewise, the computational design of 3D gradient index optics in order to shape arbitrary transverse beam profiles is continuously advancing. In Ref. [18], the authors use an irradiance mapping scheme relying on geometric optics, which is found by employing algorithms used in optimal mass transport problems, on top of which they run a multiplane Gerchberg-Saxton-like algorithm in order to account for diffraction. The proposed algorithm, however, has some constraints with regards to the smoothness and continuity of input/output functions. Furthermore, it has so far only been demonstrated for single-mode shaping, that is sculpting a specific coherent output field from a single coherent input beam.
Here we propose a new algorithmic approach for the design of 3D gradient index devices. Unlike previous methods for designing gradient index or freeform optics, we refrain from using tools used in ray-optical designs. Instead, we employ an approach related to machine learning, i.e., numerical beam propagation in conjunction with suitable cost functions that are optimized using error back-propagation and gradient descent. Such inverse design approaches based on algorithmic differentiation have already been used for designing cascaded optical systems with metasurfaces [19] or photopolymers [20]. Although similar, our approach is much more computationally intensive as it involves a much larger number of input modes and planes for an accurate representation of 3D volumes. We demonstrate the effectiveness and versatility of our approach by computing a variety of highly miniaturized 3D gradient index designs that perform key tasks in integrated photonics. After a detailed explanation of our computational method and suitable cost functions in section 2, we present solutions for three application cases in section 3.
Firstly, we find a solution for a mode sorter of only 2.5 mm length, which is capable of specifically mapping 45 Gaussian input beams onto the Hermite-Gaussian (HG) modes of a multimode fiber with an average conversion efficiency of 98.6%. Secondly, we design an equally small photonic lantern, which couples light from 45 Gaussian inputs unspecifically into a multimode fiber with 99.2% efficiency. Finally, we present a highly miniaturized beam shaper of only 0.5 mm length, which turns the irradiance profile of a monochromatic $50\,\mathrm {\mu m}$-diameter multimode fiber source into a speckle-free square of $60\,\mathrm {\mu m}$ side length.
All presented designs exhibit smooth 3D refractive index modulations within a dynamic range of about $10^{-2}$ and are therefore feasible for experimental realization using 2-photon polymerization [14].
2. Inverse design algorithm
In the following, we present an algorithm that allows design of gradient-index structures for the conversion of a set $\{u_n\}_{1\leq n \leq N}$ of $N$ mutually incoherent transverse input modes to another set of output modes $\{v_n\}_{1\leq n \leq N}$, according to various cost functions depending on the application.
We follow a typical inverse design approach by first defining a forward model propagating the input modes from their definition plane $z_{in}$ to the destination plane $z_{out}$, through a structure of refractive index contrast $\Delta \textrm {RI}(x,y,z)$ that is initially unknown and will be updated iteratively. We then define different cost functions for several classes of multimode beam shaping problems and present a generic error backpropagation algorithm to compute $\Delta \textrm {RI}$ gradients. Finally, we combine these different stages into a simple gradient descent iteration scheme to solve $\Delta \textrm {RI}$ and eventually apply some constraints to it.
2.1 Forward model
We consider a refractive index distribution $\Delta \textrm {RI}(x,y,z)$ of finite extent in a rectangular coordinate system ($O,x,y,z$), $z$ being the direction of propagation of the incoming laser modes $\{u_n\}_{1\leq n \leq N}$. The wavelength of the input and output modes in vacuum is denoted $\lambda _0$. The RI volume is discretized into $P$ transverse planes located at axial positions $\{z_p\}_{1\leq p \leq P}$ separated by a small distance $\mathrm {d}z$, as illustrated in Fig. 1. The transverse resolution is denoted $(\mathrm {d}x,\mathrm {dy})$. The refractive index $n_b$ of the medium between these different planes is assumed to be the same as for the bulk material, while each plane holds a phase mask proportional to the local $\Delta \textrm {RI}$ distribution and to $\mathrm {d}z$:
for $p\in \{1,\dots,P\}$. Then, the numerical propagation of the input modes $\{u_n\}$ through the $\Delta \textrm {RI}$ distribution is achieved with the standard split-step beam propagation method (BPM), alternating propagation in a uniform medium of refractive index $n_b$ and multiplication by a phase mask defined according to Eq. (1). A symmetric split-step scheme is enforced simply by defining the input plane position $z_{in}$ at a distance $\mathrm {d}z/2$ before $z_1$ and the output plane position $z_{out}$ at $\mathrm {d}z/2$ after $z_P$.
Fig. 1. Illustration of a Gaussian to smiley gradient index mode conversion along with its discrete formulation, consisting of $P$ planes separated by a small distance $\mathrm {d}z$.
For propagating between planes we use the well-known angular spectrum (AS) method:
2.2 Cost functions and errors for different classes of problems
Depending on the beam shaping application, different cost functions can be defined, which quantify the quality of the light transform by a single scalar value $C$. Here, we consider two types of cost functions, either distances that we want to minimize to zero, or functions of merit representing a physical quantity that we want to maximize.
In order to keep the notations simple, we still name $\{u_n\}_{1\leq n\leq N}$ the elements of the N-vector $U_{out}$ representing the input modes propagated to the plane $z_{out}$. We note $\bar {U} = (\partial C/\partial u_n)_{1\leq n \leq N}$ the error vector whose elements $\{\bar {u}_n\}$ represent the variation of the cost function with respect to each propagated mode $u_n$. The $\{u_n\}$ being complex-valued, it is important to define some consistent algebraic rules for achieving such differentiation. We use the complex representation defined in [21] along with the corresponding differentiation rules. It is worth noting that other representations and derivation rules based on Wirtinger calculus can be used [22], but they should lead to the same result in the end when error computations with respect to real-valued parameters (e.g. phases) are performed. Using such rules, we see that the errors $\{\bar {u}_n\}$ have the same dimension as the input and output modes $\{u_n\}$ and $\{v_n\}$, so we call them error modes in the following. Table 1 lists a few interesting cost functions and associated error modes for different beam shaping problems, which we describe in the following.
Table 1. Examples of cost functions and associated errors for different classes of beam shaping problems.
The first cost function we refer to as "multimode matching" is a squared $l_2$ distance between the modes $\{u_n\}$ and $\{v_n\}$ with a one-to-one correspondence. We call it that way because, in the case of lossless designs, its application matches rigorously the wavefront matching method introduced in a previous article [23] in the context of mode multiplexing with planar lightwave circuits (PLC). Indeed, a lossless design implies the conservation of $\|u_n\|_2^{2}$ during $u_n$ propagation, which leads to a simpler expression for the error vector components:
We see that, up to a constant factor, the error modes correspond to the target modes $\{v_n\}$. In this case, the application of the backpropagation of errors we develop in the next subsection leads exactly to the wavefront matching method presented in [23], even if the authors obtain this result with an alternative reasoning.
If the phase offset between each mode pair $u_n$ and $v_n$ is not relevant for the multiplexing problem, a different cost function can be introduced that we name "multimode 1–1 power coupling". This time, the cost function represents the sum of overlap integrals between input and target modes and must therefore be maximized. We note that the error vector components differ from the simplified expression of multimode matching given in Eq. (3) only by a complex constant factor corresponding to the overlap integrals and by a sign which just results from the difference between minimization and maximization problems. This slight modification pre-aligns the relative constant phases between each $u_n$ and $v_n$, which can sometimes help to reduce the inverse design complexity.
We introduce a third cost function named "multimode N–N power coupling", which aims at maximizing the total power coupling between the input set of modes $\{u_n\}$ and the target modes $\{v_n\}$ without requiring a one-to-one mapping. This metric can be useful for applications such as incoherent combining of a set of modes into a multimode fiber, or for space division multiplexing (SDM) in telecommunications [4,5], as we will illustrate in the numerical results section. Compared to the multimode 1–1 power coupling metric, this one gives much more flexibility to the optimization algorithm since it allows for mapping each input mode onto an arbitrary linear combination of target modes $\{v_n\}$. Usually, this additional degree of freedom results in smoother and more symmetrical designs.
Finally, we introduce a cost function named "multimode intensity shaping", which allows to shape the total intensity distribution $I_S$ of mutually incoherent input modes $\{u_n\}$ into a desired target intensity shape $I_T$. For this particular beam shaping application, it is not necessary to define a set of target modes $\{v_n\}$, as only $I_T$ is relevant. Moreover, the $\{u_n\}$ being mutually incoherent, their time-averaged total intensity profile is simply given by:
In the following, we describe the backpropagation algorithm that allows to compute the gradients $\partial C/\partial \Delta \textrm {RI}$ for any cost function $C$ in a unified fashion, starting from its associated and case specific error vector $\bar {U}$.
2.3 Backpropagation of errors
Algorithm S2 of Supplement 1 describes the backpropagation of the error modes vector $\bar {U}$ in order to compute the refractive index gradients associated to each plane $\nabla _\textrm {RI}[p] = \partial C/\partial \Delta \textrm {RI}[p]$ in a reverse fashion. Here, the term "backpropagation" refers to the traditional definition found in machine learning, where the chain rule is applied in order to compute the partial derivatives of the cost function with respect to some parameters of interest backwards. Following the complex representation of errors defined in [21], it appears that this procedure also corresponds exactly to the physical backpropagation of the error modes through the optical system.
2.4 Optimization algorithm
The optimization algorithm relies on gradient descent and can be summarized in the following steps:
- 1. Initialize $\Delta \textrm {RI}$.
- 2. Compute the propagated modes $U_{out}$ by calling the function propagate_forward (Algorithm S1) on the input modes $U_{in}$.
- 3. Eval $C$ and $\bar {U}$, associated to a particular beam shaping problem, using the propagated modes $U_{out}$.
- 4. Compute the refractive index gradients $\nabla _\textrm {RI}$ by calling backpropagate_gradient (Algorithm S2) on the error modes $\bar {U}$.
- 5. Update $\Delta \textrm {RI}$ with a gradient step along $\nabla _\textrm {RI}$.
- 7. If $C$ satisfies a stopping criterion then return $\Delta \textrm {RI}$, otherwise go to step 2.
Concerning step 1, we often initialize $\Delta \textrm {RI}$ to zero even if it is worth mentioning that some better initialization schemes could be used [18]. A clever initialization can help to achieve a faster convergence in case the propagated modes through a constant refractive index medium have a very low spatial overlap with the target modes.
Concerning step 5, one needs to introduce a constant step size, small enough to avoid divergence, or use an adaptive learning rate method. Of course the direction of the gradient step has to be chosen consistently, depending on whether the cost function has to be minimized or maximized.
Concerning step 6, it is worth mentioning that enforcing constraints right after performing gradient steps finds a mathematical justification in the proximal algorithms framework [24]. Moreover, a notion of Plug-and-Play Priors (PPP) has been introduced in order to use more sophisticated regularizers not necessarily corresponding to an optimization objective [27]. Regarding our gradient design, a typical constraint that we wish to apply is the restriction of $\Delta \textrm {RI}$ to a fixed interval, in order to satisfy manufacturing constraints for instance.
3. Numerical results
In this section, we present some numerical results that illustrate the versatility of the optimization algorithm in addressing different beam shaping problems involving fiber optics. We aim at designing fully integrated devices of only millimeter length, which do not require any collimation or focusing of the input and output beams. The refractive index contrast required to realize some of the proposed devices may be difficult to achieve with current glass or polymer gradient manufacturing processes, but the main aim of this work is to propose a proof of concept and give an idea of the objectives to be reached so that these types of devices can be industrially manufactured.
3.1 HG mode sorter
The first system we design is HG mode sorter, inspired by recent work [28], where the authors perform a conversion of 210 modes from a set of Gaussian input beams arranged in a triangle to a set of co-propagating HG modes. Mode sorters are a proposed solution for increasing the telecommunication bandwidth by exploiting the concept of selective mode-division multiplexing, i.e., coupling many signals into a few-mode or multimode fiber [8] with control of the particular excitation of a given output mode.
In [28], the HG conversion is performed by only 7 phase masks separated by free space propagation and is implemented experimentally by a multi-pass cavity formed by a spatial light modulator (SLM) and a flat mirror. The algorithm used in Ref. [28] is equivalent to a coordinate descent formulation of the wavefront matching method presented in [23]. It is important to underline that this complex transformation is only made possible by a clever one-to-one mapping between the cartesian coordinates of a Gaussian mode in the triangle arrangement and the order $(m,n)$ of the corresponding HG output mode. With so few phase patterns, any other arrangement would lead to poor convergence of the mode transformation.
In order to limit the computational resources to something reasonable, we restrict the mode conversion to 45 modes, which still includes interesting applications for telecommunications [29]. We set the working wavelength $\lambda _0 = 1.55\,\mathrm {\mu m}$ and use Gaussian modes of waist $\omega _{in} = 5.2\,\mathrm {\mu m}$ as inputs, typical of single mode fibers (SMF-28). Figure 2 shows the triangle arrangement of the 45 input modes, with a separation parameter $\Delta = 4\,\omega _{in}$. For the outputs, we use the HG representation of standard graded-index multimode fibers’ eigenmodes, with a maximum refractive index contrast $dn_\textrm {max} = 15\times 10^{-3}$ between the core and the cladding as in [29]. With a mode solver [30], we compute the eigenmodes for such a fiber profile and estimate the value $\omega _{out} = 7.7\,\mathrm {\mu m}$ for the HG output modes’ waist. To each of the 45 input modes, we assign one of the output fiber HG eigenmodes according to the mapping described in [28]. In short, if the Cartesian position of an input mode in the triangle arrangement (in a $45^{\circ }$ rotated frame compared to Fig. 2) is expressed as $(m\Delta, n\Delta )$, then the associated output mode is $\textrm {HG}_{m,n}$.
In simulation, the device is represented by a $512\times 512\times 250$ array, the first two dimensions corresponding to the transverse planes and the third one to the plane index. The bulk refractive index used for the propagation between planes is set to $n_b = 1.444$ corresponding to fused silica. The transverse resolution is set to $\mathrm {d}x = \mathrm {d}y = 0.65\,\mathrm {\mu m}$ and the separation between planes is $\mathrm {d}z = 10\,\mathrm {\mu m}$, leading to a total device length of 2.5 mm. At first glance, the axial resolution may appear a bit low for BPM accuracy, but we verify afterwards that it is feasible, by linearly interpolating the final simulated RI profile with a better resolution $\mathrm {d}z = 0.5\,\mathrm {\mu m}$ and by propagating the input modes through it. In the following, all the results presented concerning the performance of the different simulated devices are obtained after interpolation with the finest resolution. Limiting the number of planes in the simulation is crucial because we need to store all the input modes at each plane during the forward pass (Algorithm S1) in order to be able to compute a gradient afterwards (Algorithm S2). Even with single-precision complex numbers, this already requires 23.6 GB of memory allocation.
In order to achieve this 45-mode conversion, we run 1500 iterations of the optimization algorithm using the multimode 1-1 power coupling cost function defined in Table 1. After each iteration, we force the values of $\Delta \textrm {RI}$ to stay within a range of $12\times 10^{-3}$, which is slightly smaller than the RI contrast $dn_\textrm {max}$ required to guide the output HG modes. We observed that enforcing stricter constraints on $\Delta \textrm {RI}$ causes losses in the mode conversion due to the device’s inability to strongly guide the small modes involved. In the end of the iterations, the cost function reaches the value $C \simeq 44.37$, corresponding to an average conversion efficiency of 98.6% per mode. As in [28], we introduce the $N\times N$ transfer matrix $T$ of the device, whose elements $T_{ij} = \int {v_i^{*}u_j}$ are the overlaps between the output and the propagated input modes. We also introduce the crosstalk matrix $X = T^{\dagger } T$, whose elements $X_{ij}$ can be shown to be equal to
where $\textrm {proj}_{\{v_n\}}$ denotes the projection operation on the finite and orthogonal set of modes $\{v_n\}$. Physically, the crosstalk matrix $X$ represents the consecutive multiplexing and demultiplexing operations of the input modes $\{u_n\}$ by two identical devices having the same transfer matrix $T$. The diagonal elements of $X$ represent the system losses, and the off-diagonal terms represent the crosstalks between modes. Usually, the off-diagonal terms are normalized by dividing them by the corresponding diagonal value, in order to obtain crosstalk values independent of the power transmitted in each channel. For this particular HG mode sorter, we find a largest normalized off-diagonal term $X_\textrm {max} = -46.16$ dB, corresponding to an extremely small crosstalk between modes. The intensity and phase of each converted mode are shown in Fig. 3, where we recognize almost perfectly shaped HG modes.Fig. 3. Intensity and phase patterns of the 45 created HG modes. The average conversion efficiency is 98.6%.
Fig. 4 compares the transfer matrix $T$ and the crosstalk matrix $X$ by displaying their off-diagonal normalized square values, and their losses corresponding to the diagonal terms expressed in dB (physically the diagonal terms represent the power coupling efficiencies but we name them losses when they are expressed in dB as they also correspond to attenuation coefficients). We see that for a mode selective device like this HG mode sorter, the matrix $X$ is almost diagonal with a maximum off-diagonal term $T_\textrm {max} = -28.8$ dB. We also notice that the mode losses of the matrix $X$ are about 2 times higher than the ones of matrix $T$, which is consistent with the fact that $X$ accounts for both multiplexing and demultiplexing operations with two identical systems.
Fig. 4. Normalized off-diagonal elements of the transfer matrix $T$ and the crosstalk matrix $X$, along with the losses corresponding to the diagonal elements of the same matrices expressed in dB. The modes are ordered according to the lines of Fig. 3 from left to right.
Figure 5 illustrates horizontal and vertical sections through the center of the mode converter, indicating a very smooth evolution of the refractive index along the propagation direction. The input facet $z=z_{in}$ defines multiple waveguides for each of the 45 Gaussian input modes while the end facet $z=z_{out}$ defines a single larger waveguide holding the co-propagating HG modes. In between, the waveguides are merged by the algorithm in order to achieve the desired mode to mode mapping. We provide an additional animation of the joint evolution along $z$ of the RI profile and the total intensity of all modes in Visualization 1.
Fig. 5. Evolution of the simulated transverse RI profile sections through the center of the mode sorter device along the propagation direction, for the multimode HG conversion.
3.2 Photonic lanterns
Photonic lanterns are optical components allowing for a low-loss conversion between a set of fundamental modes belonging to independent single-mode waveguides and a set of high-order modes belonging to a common multimode waveguide [10].
The functional difference between a photonic lantern and a mode sorter as discussed previously is that the lantern is a non-selective device that does not aim for a pre-defined mode-to-mode mapping. The only goal is to couple as much power into the set of modes supported by the output fiber, while preserving mode orthogonality. Such objects are usually fabricated by tapering single-mode waveguides together or by femtosecond laser writing and they find applications in astrophotonics [31] and SDM [32] for telecommunications.
The empirical adiabatic merging of different fiber cores gives rise to super-modes by evanescent coupling, but it is challenging to ensure that this process is lossless and that the final modal content exactly matches a unitary superposition of modes belonging to a given few-mode or multimode fiber. Assuming that we can control the RI distribution as we wish, we will show that the multimode N-N power coupling metric proposed in Table 1 is perfectly adapted to overcome this difficulty. If the target waveguide was a graded-index multimode fiber with HG eigenmodes, then the HG mode sorter presented previously would already be an appropriate solution. However, it is not necessarily possible to find a Cartesian representation for the eigenmodes of an arbitrary fiber profile, one example are the linearly polarized (LP) modes of a step-index fiber.
Here, we still define the working wavelength $\lambda _0 = 1.55\,\mathrm {\mu m}$ and the bulk refractive index $n_b = 1.444$. The device is represented by a $350\times 350 \times 250$ array, with the same transverse and axial resolution as before, leading again to a total device length of 2.5 mm. Concerning the input modes, we still use the fundamental Gaussians of single-mode fibers with $\omega _{in} = \mathrm {5.2 \mu m}$, but this time with a concentric ring arrangement which has been shown to lead to the best coupling efficiency for a photonic lantern design with tapered waveguides [33]. Figure 6 represents the concentric ring arrangement for 45 input modes, with a separation parameter $\Delta = 4\,\omega _{in}$ as before.
Concerning the output modes, we use the first 45 Laguerre-Gaussian (LG) modes with $\omega _{out} = 7.7 \mu\textrm{m}$, using their real-valued representation with circular and radial nodal lines. For the multimode N-N power coupling optimization, the exact base representation is not relevant since it allows for a linear rearrangement of the modes, but we find it nicer for visualizing the intensity profiles later on.
After only 300 iterations of the optimization algorithm, we reach a total power coupling efficiency of 99.2%. For this particular non-selective device, Fig. 7 shows that the converted modes look very different from the LG modes defining the output basis. Indeed, the transfer matrix $T$ displayed in Fig. 8 is far from being diagonal and seems instead to distribute the energy almost uniformly over the LG output basis. However, the crosstalk matrix $X$ displayed in Fig. 8 is nearly diagonal with very low losses and a maximum crosstalk $X_\textrm {max} = -46.5$ dB.
Fig. 8. The first image represents the raw transfer matrix $T$. The vertical indices denote the LG modes defining the output basis as displayed in Fig. 9, and the horizontal indices denote the converted modes as displayed in Fig. 7. The modes are numbered line by line from left to right. The second image represents the normalized off-diagonal elements of the crosstalk matrix $X$, similarly as in Fig. 4. For $X$, the vertical and horizontal indices both denotes the converted modes projected on the output basis. The third image represents the losses associated with the crosstalk matrix $X$.
Usually, these non-selective devices are used in conjunction with a transmission link which is also non-selective, typically a graded index multimode fiber, and can introduce arbitrary unitary mode mixing between the multiplexing and demultiplexing stages. It is then common to use performance metrics that are independent of this arbitrary unitary mixing, such as insertion losses (IL) and mode dependent losses (MDL). As in [28], we define the IL and MDL of the device using the singular values $\{\sigma _i\}_{1\leq i \leq N}$ resulting from the singular value decomposition (SVD) of the transfer matrix $T = U\Sigma V^{\dagger }$:
For this device, we obtain $\textrm {IL} = 0.035$ dB and $\textrm {MDL} = 0.097$ dB. Since $\Sigma \simeq \mathbb {I}$, the transmission matrix $T$, representing the projection of the converted modes on the LG basis, is almost unitary and is close to $T \simeq UV^{\dagger }$. Thus, the inverse of $T$ is also almost unitary and is close to $T^{-1} \simeq VU^{\dagger }$. We can then verify that it is possible to revover the original LG output mode basis simply by rotating the converted modes displayed in Fig. 7 with the unitary matrix $VU^{\dagger }$ resulting from the SVD decomposition of the transfer matrix $T$. The result is shown in Fig. 9.
Fig. 9. Intensity and phase profiles of the converted modes rotated by the unitary matrix $VU^{\dagger }$ resulting from the SVD.
Finally, Fig. 10 shows that the computed design is very smooth and symmetrical, resembling a tapering of several single-mode waveguides, thus truly deserving the name photonic lantern. Our results prove that our method is able to engineer integrated photonic lantern devices in a very precise fashion, with perfectly controlled input and output modes. To our knowledge, this inverse design approach is new in this context. As for the HG mode sorter, an animation of the joint evolution along $z$ of the RI profile and the total intensity of all modes in the photonic lantern is available in Visualization 2.
Fig. 10. Transverse RI profiles through the center of the photonic lantern design. The total power coupling efficiency is 99.2%.
3.3 Multimode intensity shaping
Finally, we would like to treat the application of shaping the intensity profile of a light source emitting many mutually incoherent modes. This class of beam shaping has importance in various fields, such as digital holography or laser material processing, where the native intensity profile of the light source is often not ideal for the task at hand. The associated cost function is called “multimode intensity shaping” in Table 1. In a recent article, we have used the same optimization metric to find two cascaded phase patterns for shaping the intensity profile of a light emitting diode (LED) [34]. Here, we aim at realizing the same kind of intensity shaping with a volumetric gradient index design, starting from the eigenmodes of a step-index fiber with a diameter of $50\,\mathrm {\mu m}$ and $\textrm {NA} = 0.1$.
In simulation, we use $\lambda _0 = 0.640\,\mathrm {\mu m}$ as the working wavelength and $n_b = 1.457$ as the refractive index of fused silica at this wavelength. With such parameters, a mode solver finds 158 eigenmodes for this particular step-index fiber. We further assume that each input mode contributes with the same optical power. The device is represented by a $256\times 256\times 50$ array with a transverse resolution $\mathrm {d}x = \mathrm {d}y = 0.5\,\mathrm {\mu m}$ and a separation between planes $\mathrm {d}z = 10\,\mathrm {\mu m}$. This corresponds to a total device length of only 0.5 mm. As previously, we validate our sampling choice by linearly interpolating the final RI profile along the z-direction at a finer sampling interval of $\mathrm {d}z = 0.5\,\mathrm {\mu m}$, followed by a simulated readout of the gradient index design.
As desired output, we define a square target intensity profile $I_T$ with a side length of $60\,\mathrm {\mu m}$. The square is modeled as the Cartesian product of two 1D supergaussian profiles of order 20 to provide smooth edges.
For this simulation, we restrict the $\Delta \textrm {RI}$ range to $5\times 10^{-3}$ and we also enforce a smooth design, which is an example of a constraint that could facilitate the manufacturing of the device, by Fourier filtering the 3D RI profile with a sharp low-pass filter at each iteration. This Fourier filtering corresponds to a projection of the RI distribution on a set of band-limited functions, which makes it a theoretically valid proximal operator [24]. Here, we enforce the smallest details in the RI profile to be around $4\,\mathrm {\mu m}$. After 1500 iterations, even under these constraints, we reach an almost perfect intensity conversion as illustrated in Fig. 11.
Fig. 11. Theoretical target intensity profile and converted multimode intensity profile obtained by simulation.
Figure 12 shows the RI profile and the total intensity evolution of the beam along the propagation direction. An animation representing the same evolution in 3D is available in Visualization 3. We observe that the obtained RI profile is very smooth, and its limited $\Delta \textrm {RI}$ range makes it already manufacturable with commercial 3D printers with polymers [35], even if it may require some slice stitching in order to reach the full 0.5 mm length. We also see that the beam is guided in the computed RI structure while its intensity is continuously shaped into the desired profile. The fact that the total length is 5 times shorter for this multimode shaping device than for the previous mode sorter or photonic lantern shows that multimode intensity shaping requires in general much less degrees of freedom than individual or collective mode matching. The same conclusion holds for discrete designs such as presented in Refs. [34] and [28].
4. Conclusion
We have demonstrated a powerful and versatile computational method to design millimeter-range integrated devices performing several important multimode light transformations. Different cost functions were introduced for solving particular mode multiplexing or intensity shaping problems, but the optimization approach is stated in such a general way that any other user-defined cost function would fall within the scope presented here. The presented designs can already be manufactured with the available two-photon polymerization technologies and we can envision that their realization in a glass material will also be accessible in the coming years. Our work could be extended further to include broad-spectrum or multicolor sources, or to account for light polarization, which would cover many more applications in optics with fully integrated designs.
Funding
Austrian Science Fund (I3984-N36).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. H. Rubinsztein-Dunlop, A. Forbes, M. V. Berry, M. R. Dennis, D. L. Andrews, M. Mansuripur, C. Denz, C. Alpmann, P. Banzer, T. Bauer, E. Karimi, L. Marrucci, M. Padgett, M. Ritsch-Marte, N. M. Litchinitser, N. P. Bigelow, C. Rosales-Guzmán, A. Belmonte, J. P. Torres, T. W. Neely, M. Baker, R. Gordon, A. B. Stilgoe, J. Romero, A. G. White, R. Fickler, A. E. Willner, G. Xie, B. McMorran, and A. M. Weiner, “Roadmap on structured light,” J. Opt. 19(1), 013001 (2016). [CrossRef]
2. A. Forbes, M. de Oliveira, and M. R. Dennis, “Structured light,” Nat. Photonics 15(4), 253–262 (2021). [CrossRef]
3. M. Piccardo, V. Ginis, A. Forbes, S. Mahler, A. A. Friesem, N. Davidson, H. Ren, A. H. Dorrah, F. Capasso, F. T. Dullo, B. S. Ahluwalia, A. Ambrosio, S. Gigan, N. Treps, M. Hiekkamäki, R. Fickler, M. Kues, D. Moss, R. Morandotti, J. Riemensberger, T. J. Kippenberg, J. Faist, G. Scalari, N. Picqué, T. W. Hänsch, G. Cerullo, C. Manzoni, L. A. Lugiato, M. Brambilla, L. Columbo, A. Gatti, F. Prati, A. Shiri, A. F. Abouraddy, A. Alù, E. Galiffi, J. B. Pendry, and P. A. Huidobro, “Roadmap on multimode light shaping,” arXiv preprint arXiv:2104.03550 (2021).
4. G. Li, N. Bai, N. Zhao, and C. Xia, “Space-division multiplexing: the next frontier in optical communication,” Adv. Opt. Photonics 6(4), 413–487 (2014). [CrossRef]
5. B. J. Puttnam, G. Rademacher, and R. S. Luís, “Space-division multiplexing for optical fiber communications,” Optica 8(9), 1186–1203 (2021). [CrossRef]
6. S. Li, Y. Wang, Z. Lu, L. Ding, P. Du, Y. Chen, Z. Zheng, D. Ba, Y. Dong, H. Yuan, Z. Bai, Z. Liu, and C. Cui, “High-quality near-field beam achieved in a high-power laser based on SLM adaptive beam-shaping system,” Opt. Express 23(2), 681–689 (2015). [CrossRef]
7. L. Ackermann, C. Roider, and M. Schmidt, “Uniform and efficient beam shaping for high-energy lasers,” Opt. Express 29(12), 17997–18009 (2021). [CrossRef]
8. S. Gross and M. Withford, “Ultrafast-laser-inscribed 3d integrated photonics: challenges and emerging applications,” Nanophotonics 4(3), 332–352 (2015). [CrossRef]
9. R. R. Thomson, T. A. Birks, S. G. Leon-Saval, A. K. Kar, and J. Bland-Hawthorn, “Ultrafast laser inscription of an integrated photonic lantern,” Opt. Express 19(6), 5698–5705 (2011). [CrossRef]
10. T. A. Birks, I. Gris-Sánchez, S. Yerolatsitis, S. Leon-Saval, and R. R. Thomson, “The photonic lantern,” Adv. Opt. Photonics 7(2), 107–167 (2015). [CrossRef]
11. A. Ródenas, G. Martin, B. Arezki, N. Psaila, G. Jose, A. Jha, L. Labadie, P. Kern, A. Kar, and R. Thomson, “Three-dimensional mid-infrared photonic circuits in chalcogenide glass,” Opt. Lett. 37(3), 392–394 (2012). [CrossRef]
12. R. Diener, J. Tepper, L. Labadie, T. Pertsch, S. Nolte, and S. Minardi, “Towards 3d-photonic, multi-telescope beam combiners for mid-infrared astrointerferometry,” Opt. Express 25(16), 19262–19274 (2017). [CrossRef]
13. N. Riesen, S. Gross, J. D. Love, and M. J. Withford, “Femtosecond direct-written integrated mode couplers,” Opt. Express 22(24), 29855–29861 (2014). [CrossRef]
14. 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 Reviews 9(6), 706–712 (2015). [CrossRef]
15. M. Schmid, D. Ludescher, and H. Giessen, “Optical properties of photoresists for femtosecond 3d printing: refractive index, extinction, luminescence-dose dependence, aging, heat treatment and comparison between 1-photon and 2-photon exposure,” Opt. Mater. Express 9(12), 4564–4577 (2019). [CrossRef]
16. C. R. Ocier, C. A. Richards, D. A. Bacon-Brown, Q. Ding, R. Kumar, T. J. Garcia, J. van de Groep, J.-H. Song, A. J. Cyphersmith, A. Rhode, A. N. Perry, A. J. Littlefield, J. Zhu, D. Xie, H. Gao, J. F. Messinger, M. L. Brongersma, K. C. Toussaint, L. L. Goddard, and P. V. Braun, “Direct laser writing of volumetric gradient index lenses and waveguides,” Light. Sci. Appl. 9(1), 196 (2020). [CrossRef]
17. R. Dylla-Spears, T. D. Yee, K. Sasan, D. T. Nguyen, N. A. Dudukovic, J. M. Ortega, M. A. Johnson, O. D. Herrera, F. J. Ryerson, and L. L. Wong, “3d printed gradient index glass optics,” Sci. Adv. 6(47), eabc7429 (2020). [CrossRef]
18. W. M. Kunkel and J. R. Leger, “Numerical design of three-dimensional gradient refractive index structures for beam shaping,” Opt. Express 28(21), 32061–32076 (2020). [CrossRef]
19. A. S. Backer, “Computational inverse design for cascaded systems of metasurface optics,” Opt. Express 27(21), 30308–30331 (2019). [CrossRef]
20. N. U. Dinc, J. Lim, E. Kakkava, C. Moser, and D. Psaltis, “Computer generated optical volume elements by additive manufacturing,” Nanophotonics 9(13), 4173–4181 (2020). [CrossRef]
21. A. S. Jurling and J. R. Fienup, “Applications of algorithmic differentiation to phase retrieval algorithms,” J. Opt. Soc. Am. 31(7), 1348–1359 (2014). [CrossRef]
22. P. Chakravarthula, Y. Peng, J. Kollin, H. Fuchs, and F. Heide, “Wirtinger holography for near-eye displays,” ACM Trans. Graph. 38(6), 1–13 (2019). [CrossRef]
23. Y. Sakamaki, T. Saida, T. Hashimoto, and H. Takahashi, “New optical waveguide design based on wavefront matching method,” J. Lightwave Technol. 25(11), 3511–3518 (2007). [CrossRef]
24. N. Parikh and S. Boyd, “Proximal algorithms,” FNT in Optimization 1(3), 127–239 (2014). [CrossRef]
25. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. on Image Process. 18(11), 2419–2434 (2009). [CrossRef]
26. S. Lefkimmiatis, J. P. Ward, and M. Unser, “Hessian Schatten-norm regularization for linear inverse problems,” IEEE Trans. on Image Process. 22(5), 1873–1888 (2013). [CrossRef]
27. U. S. Kamilov, H. Mansour, and B. Wohlberg, “A plug-and-play priors approach for solving nonlinear imaging inverse problems,” IEEE Signal Process. Lett. 24(12), 1872–1876 (2017). [CrossRef]
28. N. K. Fontaine, R. Ryf, H. Chen, D. T. Neilson, K. Kim, and J. Carpenter, “Laguerre-gaussian mode sorter,” Nat. Commun. 10(1), 1865 (2019). [CrossRef]
29. N. K. Fontaine, R. Ryf, H. Chen, S. Wittek, J. Li, J. C. Alvarado, J. E. Antonio-Lopez, M. Cappuzzo, R. Kopf, A. Tate, H. Safar, C. A. Bolle, D. Neilson, E. C. Burrows, K. W. Kim, P. Sillard, F. Achten, M. Bigot, A. A. Correa, R. A. Correa, J. Du, Z. He, and J. Carpenter, “Packaged 45-mode multiplexers for a $50\mathrm {-\mu m}$ graded index fiber,” in 2018 European Conference on Optical Communication (ECOC), (IEEE, 2018), pp. 1–3.
30. A. B. Fallahkhair, K. S. Li, and T. E. Murphy, “Vector finite difference modesolver for anisotropic dielectric waveguides,” J. Lightwave Technol. 26(11), 1423–1431 (2008). [CrossRef]
31. B. Norris and J. Bland-Hawthorn, “Astrophotonics: The rise of integrated photonics in astronomy,” Opt. Photonics News 30(5), 26–33 (2019). [CrossRef]
32. S. G. Leon-Saval, N. K. Fontaine, J. R. Salazar-Gil, B. Ercan, R. Ryf, and J. Bland-Hawthorn, “Mode-selective photonic lanterns for space-division multiplexing,” Opt. Express 22(1), 1036–1044 (2014). [CrossRef]
33. N. K. Fontaine, R. Ryf, J. Bland-Hawthorn, and S. G. Leon-Saval, “Geometric requirements for photonic lanterns in space division multiplexing,” Opt. Express 20(24), 27123–27132 (2012). [CrossRef]
34. N. Barré and A. Jesacher, “Holographic beam shaping of partially coherent light,” Opt. Lett. 47(2), 425–428 (2022). [CrossRef]
35. X. Porte, N. U. Dinc, J. Moughames, G. Panusa, C. Juliano, M. Kadic, C. Moser, D. Brunner, and D. Psaltis, “Direct (3+1)d laser writing of graded-index optical elements,” Optica 8(10), 1281–1287 (2021). [CrossRef]
