1. Introduction

With recent advancements in quantum information and communication [1, 2] and the increasing demand on the information capacity of optical communication systems, higher-order spatial modes of laser light, such as Hermite-Gaussian (HG) modes or Laguerre-Gaussian modes are being actively considered for encoding and sending information [3–8]. The sorting and demultiplexing [9, 10] of these modes is therefore instrumental in high capacity optical communication. The ability to analyze and decompose an optical beam into a specific modal basis is also essential for optical system analysis [11–14]. Currently, there exist a variety of approaches that can handle specific pre-defined modal sets. For example, q-plates can unwrap the azimuthal phase of optical beams carrying orbital angular momentum to facilitate the sorting of these specific modes [15]. Modal content of a beam comprising of incoherent mixtures of HG modes can be evaluated using the intensity profile of the beam [16]. Mode sorting of spatial modes with well-defined parity can be realized using Mach-Zehnder and Sagnac interferometers [17–19]. Each of the previous modalities handles a specific set of modes and it would be desirable to have reconfigurable optical hardware that can sort arbitray modes on the fly — spatial, temporal, spectral and polarization. This is not only critical for studying simple optical beams, but it is essential to interrogate optical beams carrying so-called classical entanglement [20–23].

In the following, we describe an approach to arbitrary spatial mode analysis based on two-path interferometry. Recently spatial mode sorting based on controlled strong scattering was also reported in [24]. As described in [25, 26], by generalizing the delay line in a conventional Mach-Zehnder interferometer to an arbitrary unitary transformation, it becomes possible to unlock the full mode-analyzing ability of the Mach-Zehnder. Specifically, the eigenfunctions of the generalized delay line are the basis in which beam modal analysis is possible. A Fourier transform of the output interferogram reveals which eigenmodes of the generalized delay line are active in the input optical beam.

Our proof-of-principle mode-analyzing demonstration is based on a Michelson fractional Fourier transform (fFT) delay line. The fFT has a long and rich history in optics [27, 28]. It emerges naturally in considering scalar field propagation between two parallel planes [29] and also describes the propagation of fields in a quadratic GRIN medium [30]. It can be nicely visualized as a rotation of the field’s Wigner function in phase space [31] where the angle of rotation is related to the order of the fFT. Important for the following studies is the fact that the HG modes are the fFT eigenfunctions and that this can be used for analyzing the HG modal content of an optical beam [32]. When replacing the conventional temporal delay line in a Michelson interferometer, the interferometer is able to decompose an input optical beam in terms of its constituent HG modes. In the following we first describe our optical implementation of both 1D and 2D fFT. Next, we provide a proof-of-principle experiment where we demonstrate a mode-analyzing Michelson based on a fFT generalized delay line. We use this to measure the modal content of a 1D optical slit and a 2D superposition of HG modes by extracting the modal weights of the constituting HG modes. In contrast to [26] the present work utilizes a two-path Michelson interferometer and demonstrates complete two-dimensional spatial mode analysis. In this work, we report the extraction of the amplitudes of the modal weights only, their complex phase are not extracted. However, the complex phases can also be retrieved using phase retrieval techniques as we have recently shown in [33].

2. Concept

2.1. Interferometric spatial mode analysis

A Michelson interferometer is an amplitude dividing interferometer in which the two beams travel different paths (a two-path geometry) and one of the beams is delayed with respect to the other before they interfere. When a quasi-monochromatic beam with mean frequency ω and amplitude E₀ illuminates this interferometer, the part of the beam that experiences the delay (∆t) in temporal domain acquires a phase delay proportional to (−iω∆t) in the spectral domain. When this phase-delayed beam interferes with the beam from the reference arm, the resultant interferogram varies sinusoidally as |E₀(1 + e^iω∆t)|² = 2|E₀|²[1 + cos(ω∆t)] when ∆t is varied. Harmonic analysis (in the spectral domain) of this interferogram reveals a delta function at the spectral frequency ω. When an optical pulse, which can be expressed as a linear superposition of time harmonics ∫ c(ω)e^iωt dω with weights c(ω), illuminates this interferometer, harmonic analysis reveals the spectral content of the pulse |c(ω)|². This works since the time harmonics $\left\{e^{i{\omega }_{n}t}\right\}$ are the eigenfunctions of the temporal delay with eigenvalues $e^{i{\omega }_n\mathrm{\Delta }t}$. Recasting the interferogram delay as a unitary transformation provides a bridge to utilize interferometers to analyze optical beams in arbitrary basis sets.

In this current work, we will experimentally demonstrate how with a bucket detector the fractional Fourier transform, implemented in the delay arm of the Michelson interferometer, helps reveal the transverse spatial modal structure of an optical beam in the Hilbert space spanned by HG basis vectors {HG_m,n(x, y)}. Note that {HG_m,n(x, y)} are eigenfunctions of the linear canonical transform- fFT with eigenvalues $e^{i(m{\alpha }_x+n{\alpha }_y)}$ [29], where α_x,y is the parameter of the transform along x (y) dimension. HG modes of order (m, n ) are defined as $H{G}_{m,n}^{w_0}(x,y)~H_m(x)H_n(y)\exp (-\pi \frac{x^2+y^2}{w_0^2})$ where H_m(x) [H_n(y)] is the Hermite polynomial in x [y] of order m[n]. These 2D spatial modes are orthonormal and form a complete set with the Gaussian mode being the first member of this family (n =m =0). Each function HG_m,n(x, y) is a basis vector of countably infinite dimensional Hilbert space where (m, n ) ∈ (0, 1, 2, 3, …∞). Any spatial distribution U (x, y) can be written in this Hilbert space of modes as $U(x,y)={\displaystyle {\sum }_{m,n}^{\infty }{c}_{mn}H{G}_{m,n}(x,y)}$, where c_mn are the expansion coefficients or modal weights and can be complex numbers. Since HG_m,n(x, y) is not a complex function, modal weights can be computed from $c_{mn}={\displaystyle {\int }_{-\infty }^{\infty }U(x,y)H{G}_{m,n}(x,y)\mathrm{d}x\mathrm{d}y}\forall m,n=0,1,2\dots$

Figure 1 describes the schematic of a Michelson interferometer with a fFT delay line, for parameters (α_x, α_y), designed for spatial modal analysis in the HG modal basis. As shown in Fig. 1, an input beam U(x, y), which in the Hilbert space of HG modes is comprised of 40% of HG_0,2(x, y) (|c₀₂|²), 20% of HG_1,3(x, y) (|c₁₃|²) and 40% of HG_2,1(x, y) (|c₂₁|²) (approx.), illuminates the interferometer. Half of the amplitude of this beam reflects at the beam splitter onto the fFT delay line and the other half transmits to the mirror (shown as ‘m’ in the figure) and reflects back as the reference beam. When the entire beam undergoes the fractional Fourier transform with parameters (α_x, α_y), each of its modes (in the Hilbert space of modes) acquires a phase proportional to its eigenvalue: $e^{i2{\alpha }_x},e^{i({\alpha }_x+3{\alpha }_y)}$ and $e^{i(2{\alpha }_x+{\alpha }_y)}$ respectively. When this spatially transformed beam interferes with the reference beam on a bucket detector, which integrates the intensity over all space, only those constituting modes interfere with each other which have the same (mn) index (orthogonality). The resultant interferogram goes as $P({\alpha }_x,{\alpha }_y)=\left|c_{02}{|}^2\right|(1+e^{i2{\alpha }_y}){|}^2+|c_{13}{|}^2(1+e^{i({\alpha }_x+3{\alpha }_y)}){|}^2+|c_{21}{|}^2(1+e^{i(2{\alpha }_x+{\alpha }_y)}){|}^2$ which is shown as the interferogram in Fig. 1. Since this has weighted sum of three 2D sinusoids: 2(|c₀₂|² + |c₁₃|² + |c₂₁|² + |c₀₂|² cos(2α_x)+ |c₁₃|² cos(α_x +3α_y)+ |c₂₁|² cos(2α_x +α_y)), harmonic analysis of this interferogram reveals 3 delta functions at corresponding |c_mn|² and is shown as the modal weights in the same figure.

 

Fig. 1 Schematic of a Michelson interferometer designed for spatial mode analysis in the Hermite-Gaussian (HG) basis. The delay line in the interferometer is replaced with a fractional Fourier transform (fFT delay line) of parameters α_x, α_y. An exampleinput field which comprises of three HG modes is shown (intensity of this field is shown as ‘Input’). The intensities of the constituent HG modes in the input, and their weights, are shown adjacent to the input in the schematic (see text for details). The phase delays imparted to each of these HG modes by the fFT, which depend on the index (m, n) of the mode and the delay parameter α_x, α_y, are shown above the fFT delay line. The corresponding numerically simulated interferogram, recorded by the bucket detector, is shown as ‘Interferograms’. Finally, the modal weights revealed by Fourier analysis are shown next to the interferogram. Note that the fFT delay line in the interferometer is implemented using a set of three generalized lenses. 1(a) shows the schematic of this fFT set-up in transmission for simplicity.

Download Full Size | PDF

Therefore, the fractional Fourier transform is a spatial transformation which delays the optical beam in the Hilbert space spanned by discrete modes by an amount which is linear in the modal index and the order of the transform. This delay, when implemented in an interferometer, facilitates the extraction of modal weights from the interferogram by Fourier analysis. The general form of the interferogram is given as:

2.2. Generalized delay line: optical implementation of the fractional Fourier transform

The 2D fFT for a distribution f (x, y) with parameters α_x and α_y is defined as:

Our optical implementation of the fFT is based on the set-up described in [34]. Figure 1(a) shows the schematic of this fFT set-up in transmission for simplicity. This set-up implements the full fFT (including the phase of the fFT) of the input field for the transformation parameters α ∈ [π/2, 3π/2]. Implementing the fFT along with its phase requires three generalized lenses (L₁, L₂, L₃) separated by a fixed distance z where the first (L₁) and the last lens (L₃) are identical [34]. Each of these generalized lens is an assembly of two cylindrical lenses crossed at angle of π/2. The optical power of each of the generalized lenses lenses varies with the fFT parameter α_x,y and the separation between the lenses z. The power is given as:

These three generalized lenses needed to obtain the fFT of the input are implemented electronically in our set-up using SLMs. We use reflective phase only Hamamatsu SLMs having 792×800 LCD pixel displays with each pixel measuring 20 µm on each side which respond to horizontal polarization only. The phase ϕ(x, y) corresponding to the transmittance function of a generalized lens with power p_x and p_y [as defined in Eq. (3)] is written on the SLM by addressing its each pixel (x,y) using LabView. The family of HG modes that are eigenfunctions to the fFT implemented by this set-up have a waist w_o given by $\sqrt{2\lambda z}$ where z is the separation between the generalized lenses.

3. Methods

3.1. Experimental setup

The layout of our experimental set-up is given in Fig. 2. The gray box shows the Michelson interferometer where the reflected arm of the interferometer implements the fFT of the input and the transmitted arm acts as the reference. The “mode generation input” channel prepares and sends a horizontally polarized input beam either as a HG mode or a back-illuminated slit with a wavelength of 795 nm. The HG mode with appropriate waist is prepared using a hologram written on a spatial light modulator (SLM), not shown in the schematic, as described in [35]. Both the slit and the HG modes are imaged on the input plane (shown as a black dashed line) of the interferometer.

 

Fig. 2 Schematic of the optical set-up: the gray box shows the Michelson interferometer. The fFT is implemented using SLMs in the upper arm. The stabilization channel sends a second optical beam cross-polarized with respect to the signal input in order to generate an error signal which is used to path stabilize the interferometer. Experimentally generated test modes are also shown.

Download Full Size | PDF

This field distribution at the input plane of the interferometer is imaged on the first generalized lens, L₁ (implemented on SLM1) using a 4f lens system with unit magnification (f(L1)=f(L2)= 200 mm) through the reflecting port of a plate beam splitter (BS2). The second generalized lens L₂ (implemented on SLM2) is placed at a distance z = 600 mm away from the first generalized lens at normal incidence and the light is reflected back onto SLM1 (the power of the first and the third generalized lenses is same so SLM1 implements both these lenses). This field distribution on SLM1 gives the fFT of the input field. In the reference arm, the image relayed from first the 4f system (L1, L2) is relayed to the reference mirror on a piezo using another 4f system with 0.5 magnification (f(L3)=200 mm, f(L4)=100 mm). The fFT of the input interferes with the reference field on BS2 where the total field is transmitted and imaged using a third 4f lens system (f(L5)=200 mm, f(L6)=100 mm) with magnification of 0.5 on a bucket detector.

Note that the set-up described in [34] and in Fig. 1(a), only outputs the intensity corresponding to the fFT of the input. In order to obtain full fFT with phase [36], an additional global phase $e^{i({\alpha }_x+{\alpha }_y-\pi )/2}$ is added to SLM2 in our set-up. To record an interferogram P(α_x, α_y), all of the optical elements are fixed and only the phases on the two SLMs are changed. The phase displayed on each SLM corresponds to the phase transmittance t_A(x, y) of the generalized lens displayed. This transmittance depends on the power of the lens $p_{x,y}^{(L)}$ which depends upon the fFT parameter (α_x, α_y) and is given in Eq. (3):

The path length of each arm of the interferometer was 600 mm. To ensure that the two arms stay in phase during the acquisition time, we stabilize the interferometer using control loop feedback mechanism. The error signal is derived from the interference of a second optical beam collinear to the input signal and prepared in vertical polarization and Gaussian profile. A polarizing beam splitter (BS3) is placed before the bucket detector to sort the signal beam to detector D1 and the error signal to detector D2. If the two arms dephase due to slight change in path length, the intensity at D2 varies and the PID adjusts the piezo on the reference mirror to change the reference arm path length accordingly.

3.2. Numerical modelling

To compute the interferograms for various inputs into our mode analyzer, we numerically modeled the propagation of the optical beam through our set-up in MATLAB. The fractional Fourier transform of the input distribution was computed by modeling the phase delays imparted by each of the three generalized lenses [Eq. (3)] and each of the two free space propagation intervals. We used the Fresnel transfer function propagator (Fourier transform of the impulse response of free space propagation under Fresnel’s approximation), H (f_x, f_y), to model propagation over each free space interval from the input plane U₁(x, y) to the output plane [U_out(x, y)] [37]:

4. Results and discussion

4.1. 1D mode analysis

Figure 3(a), 3(d) and 3(g) show the input intensities corresponding to HG_1,0(x), HG_2,0(x) and HG_10,0(x) respectively. These modes are eigenfunctions of the fFT and preserve their spatial profile as the fFT parameter α is swept from 0 to 2π and acquire a phase corresponding to their eigenvalues of e^iα, e^i2α and e^i10α. Therefore, as explained in Sec 2.1 and as shown from theoretical modeling in the upper panels of Figs. 3(b), 3(e) and 3(h), the interferograms are pure sinusoids varying as cos(α), cos(2α) and cos(10α). As shown in lower panels in Figs. 3(b), 3(e) and 3(h) experimentally generated interferograms also exhibit sinusoidal variations of one full cycle for HG_1,0(x), two cycles for HG_2,0(x) and ten cycles for HG_10,0(x). Modal weights, |c_n|², were computed by fitting these experimentally generated interferograms with a sum of cosines function. Figs. 3(c), 3(f) and 3(i) show the modal weights for first ten HG_n,0(x) modes where the ordinate shows the square of modal weights, |c_n|². For HG_1,0(x) input, we obtained a maximum contribution of |c₁|²=0.8 from n=1 mode and for the HG_2,0(x) input, |c₂|²=0.85 was obtained from n=2 mode and for HG_10,0(x) maximum contribution from n = 10 mode was obtained with |c₁₀|²=0.91.

 

Fig. 3 Results for 1D mode analysis: First column shows the intensities of the various input distributions as recorded on the CCD (located in place of the bucket detector): 3(a) HG_1,0(x), 3(d) HG_2,0(x), (g) HG_10,0(x) and 3(j) 1D slit, rect(x/a), with width a= 2.3 mm. Numerically modeled interferograms are shown in red circles in upper panels and experimentally generated interferograms are shown in blue circles in lower panels of 3(b), 3(e), 3(h) and 3(k) for each of the inputs respectively. The red lines in the experimentally generated interferograms in 3(b), 3(e) and 3(h) show the fits from which the modal weights are extracted. Extracted modal weights are shown in 3(c), 3(f), 3(i) and 3(l) (red bars show theoretical values of |c_th|² for slit input). The error in the measured modal weights for the slit input are shown in blue bars in lower panel of 3(l).

Download Full Size | PDF

Next, we used the mode analyzer to measure the modal deconstruction of a rect(x/a) function, i.e. 1D slit, which is not an eigenfunction of the fFT. The number of HG_n,0(x) modes needed to completely describe the rect(x/a) function depends on the width of the slit and the waist of the HG_n,0(x) modes used as the basis vectors. In this experiment, we used a slit with width 2.3 mm. The upper panel of Fig. 3(h) shows the numerically modeled interferogram for this slit and the lower panel shows the experimentally generated interferogram which had a similar functional dependence on the fFT parameter α as expected from the numerical modeling. Figure 3(i) shows both the numerically modeled modalweights, |c_th|², as obtained form $c_{th}={\displaystyle \int \text{rect}(x/a)H{G}_n^{w_0}(x)dx}$ in red bars and the modal weights extracted from the experimentally generated interferogram, |c_n|², in yellow-blue bars. The lower panel of Fig. 3(i) shows the difference between these modal weights defined as Δ|c|²=|c_th|² − |c_n|².

4.2. 2D mode analysis

Figure 4 shows the results for 2D optical field inputs to the mode analyzer. The input intensity distribution is shown corresponding to two HG_m,n(x, y) modes: Fig. 4(a) HG_2,2(x, y) and 4(d) HG_4,1(x, y). A 2D fFT of this input was implemented in the delay line by sweeping the fFT paramter α_x for each α_y. P(α_x, α_y) behaves like a 2D sinusoid function of the form cos(2x + 2y) and cos(4x + y) for these two HG_2,2(x, y) and HG_4,1(x, y) input states respectively. This was also achieved in the numerical modeling of the interferograms and is shown in the left column of Figs. 4(b) and 4(e) where the maximum/minimum intensity value of P(α_x, α_y) is color coded as yellow and black respectively. The experimentally generated interferograms also showed similar 2D sinusoidal behavior and are shown in the second column of Figs. 4(b) and 4(e) corresponding to HG_2,2(x, y) and HG_4,1(x, y) respectively. The extracted modal weights, |c_mn|² for these two input states, corresponding to the 25 HG_m,n(x, y) modes with m, n ∈ (0, 4), are shown in Figs. 4(c) and 4(f). The modal weights in Figs. 4(c) and 4(f) show strongest contributions from m=2, n=2 and m=4, n=1 modes as expected.

 

Fig. 4 Results for 2D mode analysis: First column shows the intensities of the various input distributions: 4(a) HG_2,2(x, y), 4(d) HG_4,1(x) and 4(i) superposition state $L{G}_{0,1}^w(x,y)=H{G}_{1,0}(x,y)+iH{G}_{1,0}(x,y)$, with w=1.6 mm. Second and third columns show the numerically modeled interferograms and experimentally generated interferograms in 4(b), 4(e) and 4(h) for each of the inputs respectively. Modal weights extracted from data are shown in 4(d), 4(f) and 4(i).

Download Full Size | PDF

Next, we measured the modal weights of a 2D superposition state of multiple HG modes. We used $L{G}_{0,1}^w(x,y)$ mode as the input state with waist w=1.6w₀. For w= w₀ (eigen waist for the HG modes), this state is a superposition of first two HG modes HG_1,0(x, y) and HG_0,1(x, y) only. However, as the waist is increased, higher order HG modes also contribute to the field. The intensity of the input mode as recorded on a CCD is shown in Fig. 4(g). The numerically modeled interferogram and the experimentally measured interferogram are shown in Fig. 4(h). Modal weights extracted from the experimentally generated interferogram are shown in Fig. 4(i). The strongest contributions from the m =1, n =0 and m =0, n =1 modes can be seen with |c₁₀|² =0.32 and |c₀₁|² =0.37 as compared to theoretically expected values of 0.32 each.

4.3. Discussion

Note that both 1D and 2D interferograms have an overall shift, for example in Fig. 4 (e), the interfered intensity should go as cos(2α) but it is measured to be shifted by ∝ π/4. This phase shift results from stabilizing the path difference between the two arms of the interferometer at a constant value. Also, note that both 1D and 2D interferograms appear distorted for certain fFT parameters around values 0 and 2π. We believe this is a result of very short focal lengths of the generalized lenses required to implement those orders. Since the SLMs are pixelated, spatial phase transformations that can be encoded on the SLM are limited by the Nyquist sampling limit which depends on the size of the SLM pixel, number of the pixels and the wavelength [26]. Regarding signal-to-noise ratio, we operate the system such that constructive interference saturates the CCD. Using the characterized detector noise, we estimate a maximum SNR of ≈ 30.

We have demonstrated an optical mode analyzer capable of measuring the spatial modal content in Hermite-Gaussian basis of an optical beam with a bucket detector. Our experimental results support and agree with the theoretical modeling of the interferograms and the modal weights. Specifically, we used this set-up to experimentally extract modal weights of 1) 1D superposition state, a slit and 2) 2D superposition of two first order HG modes. To the best of our knowledge, this is a first demonstration of a two-path Michelson interferometer capable of doing such a spatial modal analysis. As mentioned earlier, our mode analyzer only outputs the amplitudes |c_mn| and not the complex phase φ_mn of the modal weights $|c_{mn}{|}^2=|c_{mn}{e}^{i{\phi }_{mn}}{|}^2$. This limitation, however, can be overcome by using phase retrieval techniques. We have reported in [33] on retrieving phase using an array detector, by minimizing the error function constructed from measured interferogram and an estimated interferogram with c_n phases, using a nonlinear gradient-based optimization to obtain the phase of c_n from the experimentally measured c_n amplitudes. Additionally, we have also successfully created an algorithm to retrieve phase from a bucket detector [38].

5. Conclusions and future work

We have built a Michelson interferometer with a non-conventional delay line- fractional Fourier transform. This set-up enables analysis of the transverse structure of an optical beam in the Hilbert space spanned by Hermite-Gaussian modes. We used this mode analyzer to measure the modal wights of 1D superpositions and 2D superpositions of Hermite-Gaussian modes. In the future, we will report on phase retrieval algorithms to recover the phase information of the modal weights engaged by a given optical beam.