Gaussian Beam-Propagation Theory for Nonlinear Optics - Featuring an Exact Treatment of Orbital Angular Momentum Transfer

TL;DR

This paper introduces an exact Gaussian mode propagation formalism for nonlinear optics that accurately models orbital angular momentum transfer, enabling optimization and better understanding of multi-mode beam generation.

Contribution

It provides the first exact solution for Gaussian beam propagation in nonlinear optics with orbital angular momentum transfer, improving upon previous phenomenological models.

Findings

Excellent agreement with two experimental cases

Method enables optimization of nonlinear interactions

Applicable to various laboratory configurations

Abstract

We present a general, Gaussian spatial mode propagation formalism for describing the generation of higher order multi-spatial-mode beams generated during nonlinear interactions. Furthermore, to implement the theory, we simulate optical angular momentum transfer interactions, and show how one can optimize the interaction to reduce the undesired modes. Past theoretical treatments of this problem have often been phenomenological, at best. Here we present an exact solution for the single-pass no-cavity regime, in which the the nonlinear interaction is not overly strong. We apply our theory to two experiments, with very good agreement, and give examples of several more configurations, easily tested in the laboratory.

Equations

$${\nabla }^2\tilde{\text{E}}-\frac{1}{c^2}\frac{{\partial }^2\tilde{\text{E}}}{\partial t^2}=\frac{1}{{ϵ}_0{c}^2}\frac{{\partial }^2\tilde{\text{P}}}{\partial t^2},$$

$$(\frac{\partial }{\partial z}-\frac{i}{2k}{\nabla }_{\perp }^2)\mathcal{E}=\frac{ik}{2{ϵ}_0}\mathcal{P}.$$

$$(\frac{\partial }{\partial z}-\frac{i}{2k}{\nabla }_{\perp }^2)\mathcal{E}=\frac{ik}{2{ϵ}_0}\mathcal{P}({\mathcal{E}}_0)\equiv \wp ({\mathcal{E}}_0).$$

$$\begin{array}{rl}\mathrm{DE}:& \hat{L}\mathcal{E}=\wp ({\mathcal{E}}_0)\\ \mathrm{IV}:& {\mathcal{E}}_0,\end{array}$$

$$\begin{split}u_{\ell,p}(\vec{r})=&\dfrac{C_{\ell,p}}{w(z)}e^{-\frac{r^{2}}{w(z)^{2}}}e^{-\frac{ikr^{2}z}{2(z^{2}+z_{R}^{2})}}\big{(}\dfrac{\sqrt{2}r}{w(z)}\big{)}^{|\ell|}\\ \times&\>L_{p}^{|\ell|}\big{(}\dfrac{2r^{2}}{w(z)^{2}}\big{)}e^{i\ell\phi}e^{i(2p+|\ell|+1)\arctan(z/z_{R})},\end{split}$$

$$\int rdrd\varphi u_{\ell ,p}^{\ast }(r,\varphi ,z)u_{q,n}(r,\varphi ,z)={\delta }_{\ell q}{\delta }_{pn},$$

$$\begin{array}{r}B(r,\varphi ,z)=\sum _{\ell ,p}{c}_{\ell ,p}(w_0)u_{\ell ,p}(r,\varphi ,z,w_0),\end{array}$$

$$c_{\ell ,p}(w_0)=\int rdrd\varphi u_{\ell ,p}^{\ast }(r,\varphi ,z_0,w_0)B(r,\varphi ,z_0).$$

$$\begin{split}B(r,\phi,z)=&\int r^{\prime}dr^{\prime}d\phi^{\prime}\,B(r^{\prime},\phi^{\prime},z_{0})\\ *&\Big{(}\sum_{l,p}\,u_{l,p}^{*}(r^{\prime},\phi^{\prime},z_{0})u_{l,p}(r,\phi,z)\Big{)}.\end{split}$$

$$\sum _{\ell ,p}{u}_{\ell ,p}^{\ast }(r^{\prime },{\varphi }^{\prime },z)u_{\ell ,p}(r,\varphi ,z)=\delta (r-r^{\prime })\delta (\varphi -{\varphi }^{\prime }).$$

$$\begin{array}{rl}\mathcal{E}& =\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }K(\text{r}|{\text{r}}^{\prime }){\mathcal{E}}_0({\text{r}}^{\prime }){|}_{z=z^{\prime }}\\ & +\int d{z}^{\prime }\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }G(\text{r}|{\text{r}}^{\prime })\wp ({\text{r}}^{\prime }).\end{array}$$

$$\begin{array}{rl}K(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime })& \equiv \sum _{\ell ,p}{u}_{\ell ,p}^{\ast }(r^{\prime },{\varphi }^{\prime },z^{\prime })u_{\ell ,p}(r,\varphi ,z)\\ G(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime })& \equiv \Theta (z-z^{\prime })K(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime }),\end{array}$$

$$\wp (r,\varphi ,z)=\sum _{l,p}{c}_{l,p}(z)u_{l,p}(r,\varphi ,z),$$

$$\mathcal{E}(\text{r})={\mathcal{E}}_0(\text{r})+\sum _{l,p}{u}_{l,p}(\text{r}){\int }_{z_i}^{z_f}d{z}^{\prime }{c}_{l,p}(z^{\prime }).$$

$$\begin{array}{rl}\mathrm{DE}:& \hat{L}\mathcal{E}=0\\ \mathrm{IV}:& {\mathcal{E}}_0.\end{array}$$

$$\hat{L}K(\text{r}|{\text{r}}^{\prime })=0(z>z^{\prime })$$

$$K(\text{r}|{\text{r}}^{\prime })=\delta (r-r^{\prime })\delta (\varphi -{\varphi }^{\prime })(z=z^{\prime })$$

$$K\to 0(r\to \infty ),$$

$${\mathcal{E}}_{\mathrm{hom}}=\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }K(\text{r}|{\text{r}}^{\prime }){\mathcal{E}}_0({\text{r}}^{\prime }){|}_{z=z^{\prime }}.$$

$$\begin{split}\hat{L}\mathcal{E}_{\mathrm{hom}}=&\int r^{\prime}dr^{\prime}d\phi^{\prime}\,\big{(}\hat{L}\;K(\textbf{r}\;|\;\textbf{r}^{\prime})\big{)}\,\mathcal{E}_{0}(\textbf{r}^{\prime})\\ =&\int r^{\prime}dr^{\prime}d\phi^{\prime}\,\big{(}0\big{)}\,\mathcal{E}_{0}(\textbf{r}^{\prime})\\ =&0.\end{split}$$

$$K(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime })\equiv \sum _{\ell ,p}{u}_{\ell ,p}^{\ast }(r^{\prime },{\varphi }^{\prime },z^{\prime })u_{\ell ,p}(r,\varphi ,z).$$

$$\displaystyle\hat{L}\;K=\sum_{\ell,p}\,u_{\ell,p}^{*}(\textbf{r}\,^{\prime})\big{(}\hat{L}\,u_{\ell,p}(\textbf{r})\big{)}=0\quad(z>z^{\prime})$$

$$K(\text{r}|{\text{r}}^{\prime })=\delta (r-r^{\prime })\delta (\varphi -{\varphi }^{\prime })(z=z^{\prime })$$

$$u_{\ell ,p}\to 0⟹K\to 0(r\to \infty ).$$

$$\hat{L}G(\text{r}|{\text{r}}^{\prime })=\delta (r-r^{\prime })\delta (\varphi -{\varphi }^{\prime })\delta (z-z^{\prime })$$

$$G(\text{r}|{\text{r}}^{\prime })=0(z<z^{\prime })$$

$$G\to 0(r\to \infty ),$$

$$\mathcal{E}=\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }d{z}^{\prime }G(\text{r}|{\text{r}}^{\prime })\wp ({\text{r}}^{\prime }).$$

$$\begin{split}\hat{L}\mathcal{E}&=\int r^{\prime}dr^{\prime}d\phi^{\prime}dz^{\prime}\,\big{(}\hat{L}\,G(\textbf{r}\,|\,\textbf{r}^{\prime})\big{)}\,\wp(\textbf{r}^{\prime})\\ &=\int r^{\prime}dr^{\prime}d\phi^{\prime}dz^{\prime}\,\delta(\textbf{r}-\textbf{r}^{\prime})\,\wp(\textbf{r}^{\prime})\\ &=\wp(\textbf{r}).\end{split}$$

$$G(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime })\equiv \Theta (z-z^{\prime })K(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime }),$$

$$\begin{array}{rl}\hat{L}G& =(\frac{\partial }{\partial z}-\frac{i}{2k}{\nabla }_{\perp }^2)\Theta (z-z^{\prime })K(r,\varphi ,z|r^{\prime },{\varphi }^{\prime },z^{\prime })\\ & =K\frac{\partial }{\partial z}\Theta (z-z^{\prime })+\Theta \frac{\partial }{\partial z}K-\Theta \frac{i}{2k}{\nabla }_{\perp }^{2}K\\ & =K\delta (z-z^{\prime })+\Theta (z-z^{\prime })(\frac{\partial }{\partial z}K-\frac{i}{2k}{\nabla }_{\perp }^{2}K)\\ & =K\delta (z-z^{\prime })+\Theta (z-z^{\prime })\hat{L}K\\ & =\delta (r-r^{\prime })\delta (\varphi -{\varphi }^{\prime })\delta (z-z^{\prime }),\end{array}$$

$$\begin{array}{rl}\mathcal{E}& =\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }K(\text{r}|{\text{r}}^{\prime }){\mathcal{E}}_0({\text{r}}^{\prime }){|}_{z=z^{\prime }}\\ & +{\int }_{z_i}^{z_f}d{z}^{\prime }\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }K(\text{r}|{\text{r}}^{\prime })\wp ({\text{r}}^{\prime })\\ & ={\mathcal{E}}_0\\ & +\sum _{l,p}{u}_{l,p}(\text{r}){\int }_{z_i}^{z_f}d{z}^{\prime }\int r^{\prime }d{r}^{\prime }d{\varphi }^{\prime }{u}_{l,p}^{\ast }({\text{r}}^{\prime })\wp ({\text{r}}^{\prime }).\end{array}$$

$$\wp (r,\varphi ,z)=\sum _{l,p}{c}_{l,p}(z)u_{l,p}(r,\varphi ,z),$$

Full text

I Introduction

We are on the cusp of a new age of quantum physics and technology, where multi-spatial-mode beam propagation will play an ever more essential role. Many of the most important quantum resources are produced during nonlinear light-matter interactions. It is particularly interesting to study nonlinear effects in response to beams carrying orbital angular momentum (OAM). For example, the Laguerre-Gauss (LG) spatial modes have an azimuthal phase dependence of $\exp[i\ell\phi]$, which corresponds to OAM of $\ell\hbar$ per photon Allen et al. (1992); Yao and Padgett (2011). The intensity patterns for several LG modes, with $\ell\neq 0$, are given in Fig. 1, and a more complete description can be found Sec. II C. Conservation and storage (via slow and stopped light) of OAM has been realized in several processes including the entanglement of OAM modes in parametric down conversion (PDC) Mair et al. (2001); Caetano et al. (2002), second harmonic generation (SHG) Courtial et al. (1997), and four-wave mixing (FWM) in semiconductors Ueno et al. (2009). Unlike solid state processes, nonlinear optics in atomic vapors is highly efficient and requires low-light intensities. Transfer of OAM into atomic media Pugatch et al. (2007); Boyer et al. (2008); Cao et al. (2014); de Oliveira et al. (2015), transfer to frequency converted light Tabosa and Petrov (1999); Ding et al. (2012); Walker et al. (2012); Gariepy et al. (2014); Zhdanova et al. (2015); Akulshin et al. (2015, 2016), and amplification Borba et al. (2016) have all been observed in atomic vapors.

It has been observed that the transfer of OAM to frequency converted light in atomic media is typically accompanied by a disturbance of the radial component structure, i.e., frequency converted beams have nonzero radial index $p$, even when the pump beams have $p=0$ Walker et al. (2012); Zhang et al. (2016); Swaim et al. (2017). See Fig. 1 for examples of $p$ mode structure. These effects may be subtle, or even negligible, when analysis is limited to comparing images collected with a charge-coupled device (CCD) camera. However, these effects can be quite detrimental to quantum processes such as squeezing Zhang et al. (2016) and entanglement generation Dada et al. (2011); Romero et al. (2012); Fickler et al. (2012); Gupta et al. (2016). Nonlinear processes leading to the production of LG modes, with nonzero radial index $p$, can actually contaminate the mode structure and degrade the performance of the process. As for the $\ell$ structure, the conservation of OAM dictates the allowed $\ell$ modes of the output beam, but this is not the end of the story. Transfer to some of these modes can be quite improbable Walker et al. (2012); Akulshin et al. (2015), even though the interaction would not violate OAM conservation. Therefore, an analysis simply based on conservation can be quite misleading. With these concerns in mind, we develop a relatively simple, yet analytic semiclassical theory, which predicts the spatial mode structure of beams created during nonlinear interactions. We show that the predictions of our theory have excellent agreement with the observations we just discussed. Furthermore, the theory allows us to study and optimize the interaction to enhance performance and output mode structure. Finally, an accurate description of the mode structure is a crucial ingredient in the general second-quantization procedure which we are currently developing.

II Theoretical Background

In the following section, we present an introduction to spatial mode propagation, properties of the Laguerre-Gauss modes, and an overview our method for solving the inhomogeneous wave equation in the paraxial approximation. An analytic method is preferred, because we wish to explicitly retain the mode structure of the resulting beam, such that one can use the mode structure, for example, to perform a second-quantization procedure. Please see the appendix for full details of our Green function solution method.

II.1 Spatial Mode Propagation

The paraxial wave equation and its solutions are certainly well known, but for completeness we quickly derive a form useful for our purposes. We investigate the electric field propagating through a nonlinear medium, which for simplicity, we treat as lossless and dispersionless. Treating the light classically, we start with the inhomogeneous wave equation

[TABLE]

where the tilde indicates rapidly varying quantities. Introducing the field $\tilde{\textbf{E}}(\textbf{r},t)=\textbf{E}(\textbf{r})e^{-i\omega t}$ and polarization $\tilde{\textbf{P}}(\textbf{r},t)=\textbf{P}(\textbf{r})e^{-i\omega t}$ into Eq. (1), one finds the Helmholtz equation. Next, we let $\textbf{E}(\textbf{r})=\boldsymbol{\mathcal{E}}(\textbf{r}_{\perp},z)e^{ikz}$, $\textbf{P}(\textbf{r})=\boldsymbol{\mathcal{P}}(\textbf{r}_{\perp},z)e^{ikz}$, and along with the paraxial approximation ($|k\,\partial\mathcal{E}/\partial z|\gg|\partial^{2}\mathcal{E}/\partial z^{2}|$) we transform the Helmholtz equation into the inhomogeneous paraxial wave equation

[TABLE]

Depending on the scenario, there are several methods used for deriving the polarization $\boldsymbol{\mathcal{P}}$, and in many cases our simple phenomenological approach will lead to interesting and accurate results. Regardless of the method used, the polarization is, in general, a complicated function of the fields, that is, $\boldsymbol{\mathcal{P}}=\boldsymbol{\mathcal{P}}(\boldsymbol{\mathcal{E}})$. This makes Eq. (2) similar in structure to the nonlinear Schr$\ddot{\mathrm{o}}$dinger equation. With this similarity in mind, we make a first order Born approximation Saleh et al. (1991), and replace $\boldsymbol{\mathcal{E}}$ with the input beam $\boldsymbol{\mathcal{E}}_{0}$ instead, which establishes a much simpler differential equation of the form,

[TABLE]

The first born approximation is sufficient when the pump beam is only slightly modified and the seed or signal fields are weak compared to the driving fields. Fortunately, most nonlinear processes that involve generation of new fields start with vacuum modes. Furthermore, the nonlinear processes that we consider in this article satisfy these conditions. For example, in the polarization self-rotation squeezing scheme Matsko et al. (2002); Ries et al. (2003); Zhang et al. (2013), the seed beam is just vacuum entering into an empty port. Thus, in the Born approximation, the optical response of the nonlinear material inherent, in $\boldsymbol{\wp}(\boldsymbol{\mathcal{E}}_{0})$, behaves as a source for new mode components of the field. Regarding the vector nature of the interaction, it must must be noted that in general, the nonlinear susceptibility governing the interaction has a tensor form with, e.g., 81 terms for third-order nonlinearities. This structure offers a complication, which can lead to a set of coupled equations that are very difficult to solve. Fortunately, the majority of nonlinear processes that involve generation of new fields have symmetry, which vastly reduces the complication of the problem. This fact, paired with the Born approximation, is used to simplify the problem to uncoupled differential equations, which then can be solved using our method.

II.2 Construction of the initial-value-problem

In the first Born approximation, the right-hand-side of Eq. (3) is effectively the source of the new beam. Thus, it is helpful to restate the problem here as an initial-value problem (IVP) in the compact form,

[TABLE]

where $\hat{L}\equiv\frac{\partial}{\partial z}-\frac{i}{2k}\nabla^{2}_{\perp}$ and $\mathcal{E}_{0}$ is the input-field amplitude which is assumed to be undepleted. Although this IVP represents a first Born approximation, the following method can be used iteratively for cases when the nonlinearity and generated fields are stronger. In such cases, it may be necessary to include an attenuation factor, on the initial value, which conserves energy. Ideally, the theory does not require any free parameters for the source $\wp$, and can be used directly. In reality, however, it may be difficult to calculate the strength of the source if the model used does not capture the full complexity of the experimental system. Therefore, it maybe convenient to capture such complexity with an interaction strength factor that, if necessary, can be used as a fitting parameter. For example, a fitting parameter is quite useful when simulating the generation of squeezed light in hot atomic vapor, where the direct calculations are too complex Zhang et al. (2016).

II.3 Laguerre-Gauss Mode Properties

We elect to consider input beams with cylindrical symmetry, but we note that the following calculation can certainly be done in other coordinate systems. In cylindrical coordinates, the homogeneous paraxial wave equation gives rise to the LG family of solutions Siegman (1986):

[TABLE]

where $\ell$ is the azimuthal index, $p$ is the radial index for each mode, $C_{\ell,p}=\sqrt{2p!/\pi(|\ell|+p)!}$ is a normalization constant, $w_{0}$ is the beam waist, $w(z)=w_{0}\sqrt{1+(z/z_{R})^{2}}$ is the width function of the beam, $L_{p}^{|\ell|}$ are the generalized Laguerre polynomials, $z_{R}=\pi w_{0}^{2}/\lambda$ is the Rayleigh range, and $k=2\pi/\lambda$ is the wave number. The LG modes form a complete orthonormal set and thus can be used as a basis set to expand an arbitrary paraxial beam $B=B(r,\phi,z)$ in free space. Using the orthogonality relation

[TABLE]

we can write $B(r,\phi,z)$ as

[TABLE]

where

[TABLE]

The waist $w_{0}$ of the basis set is, in general, chosen to give the best fit and reduce the number of terms in the expansion, whereas the $c_{\ell,p}$ coefficients are independent of the position $z_{0}$. If we insert Eq. (8) in Eq. (7) and collect the LG modes we find

[TABLE]

If we impose $z=z_{0}$, then we can establish the very important completeness relation:

[TABLE]

The condition $z=z_{0}$ simply states that the LG modes are complete at equal $z$’s, i.e, when the $z$-slices chosen coincide. Furthermore, we will see that this condition reemerges as part of the mechanism which introduces the input beam in the solution of the IVP. Thus, the completeness relation is instrumental in the Green function solution method, which we will now discuss briefly.

II.4 Green Function Solution

The magic of the Green function solution method (Barton, 1989) is that once the propagator $K$, and Green function $G$ are derived, the problem is solved (for full details of the solution method see the appendix):

[TABLE]

Although we solve this problem in free space, the properties of the LG modes allow us to utilize a method developed for fixed boundary conditions. Since the LG modes are a complete orthonormal set, we can use an eigenfunction expansion. Recalling the completeness relation in Eq. (10), we define our propagator and Green function in terms of the LG modes:

[TABLE]

where $\Theta(z-z^{\prime})$ is the Heaviside step function.

At this point the problem is solved. However, we may further simplify matters by expanding the source $\wp$ in terms of the LG modes:

[TABLE]

where $c_{l,p}(z)=\int rdrd\phi\,u_{l,p}^{*}(r,\phi,z)\wp(r,\phi,z)$. To find the final form of our solution, we insert Eq. (12) and Eq. (13) into Eq. (11) and find

[TABLE]

So we see that the final solution is a superposition of the unmodified pump beam, and a collection of LG modes, which, referring to Eq. (3), will depend on the specified polarization $P$ and the spatial structure of the input beam $\mathcal{E}_{0}$. The cautious reader may be concerned that we have solved a problem concerning a possibly localized source distribution using a free space Green function method. We pause to point out that the effects of the source are totally subsumed in the spacial distribution of $\wp$. Furthermore, in most situations a Gaussian pump beam will be completely encompassed by the interaction region, effectively making the boundary in $r$ infinite.

III Simulation of Experiment

We start will several simulations which showcase our theory. We first focus on experimentally relevant simulations related to PDC and spontaneous FWM. For PDC, we take the results of the experiment a step further by suggesting ways to enhance the mode structure of the generated beam. For FWM, we address an unresolved problem pertaining to the pathways of OAM transfer in spontaneous FWM. For each simulation, we will specify the polarization governing the interaction, along with the relevant beam and material parameters. Furthermore, all beam profiles are plotted at $z=0$, and are color coded to match the visible spectrum as closely as possible.

III.1 Stimulated Down Conversion

First, we will investigate OAM transfer in frequency down-converted light. A rapidly evolving body of theoretical work is calling for entanglement in quantum systems with higher dimensions Gröblacher et al. (2006); Mirhosseini et al. (2015). Thus, as an exercise, we will study stimulated PDC as a way to understand the properties of down-converted twin beams when a cavity is not present to alter the mode structure. We model a PDC experiment, in which a $442\,\mathrm{nm}$ pump beam interacts with a $845\,\mathrm{nm}$ signal beam (in a $3\,\mathrm{mm}$ long BBO crystal) to create a $925\,\mathrm{nm}$ idler beam Caetano et al. (2002). Due to the small incident angle of the signal beam and the relatively thin crystal, effects from the non-collinear geometry are negligible (see appendix). The light-matter interaction in the crystal is governed by the polarization $\mathcal{P}(\omega_{i})^{(2)}=\epsilon_{0}\chi^{(2)}\mathcal{E}_{p}(\omega_{p})\mathcal{E}_{s}(\omega_{s})^{*}$, depicted in Fig. 2(a). In this interaction $E_{p}$ is the pump beam, $E_{s}$ is the signal beam, and the conjugation of $E_{s}$ corresponds to the creation of a photon in the input-signal mode. Therefore, we see how the choice of pump and signal modes tailor the response of the idler beam; this phenomenon has been verified experimentally for relatively simple input beams Caetano et al. (2002).

The waist of the basis set, used to represent the generated beam, is set to match the waist of the pump beam. In this simulation, we will see how the interplay of the waist and Rayleigh range, of the input beams, effects the mode structure of the generated beam.

In this simulation, we investigate how to generate a high-OAM superposition in the idler mode. First, for reference, we consider that the signal beam is prepared in the $u_{0,0}$ mode with a waist that matches the pump beam. As one can see in Fig. 2(b), the idler beam responds in $\ell=\pm 2$ superposition, with a spreading in $p$ modes, corresponding to the partial overlap of the modes in $E_{p}(\omega_{p})E_{s}(\omega_{s})^{*}$. Next, we investigate how to increase the OAM in the idler beam. We prepare a signal beam in the same superposition as the pump, but include a $\pi$ phase shift between the two modes in the superposition, i.e., the signal beam is now $u_{2,0}+e^{i\pi}u_{-2,0}$. In Fig. 2(c) we see that this combination creates a destructive interference, which suppresses response at $\ell=0$, and produces an idler beam, which responds at $\ell=\pm 4$ with a spreading in $p$ modes. More simulations using this technique can be found in the following section. Again, the spreading in $p$ modes corresponds, in part, to the poor overlap of the beams. Therefore, in Fig. 2(d), we show that one can optimize the signal-beam waist to suppress the higher order $p$ modes. The narrowing of the signal beam waist $w_{0}$ (by $\approx\,40\%$), in Fig. 2(d), corresponds to a reduction of the Rayleigh range $z_{R}\equiv\pi w_{0}^{2}/\lambda$. In effect, the beams now expand closer to the same rate, which is apparently the optimized beam overlap through the interaction region. This observation is related to the Boyd criterion (Walker et al., 2012), which states that, the nonlinear interaction is strongest when the Rayleigh ranges of the interacting beams are identical. The effect observed here extends this observation to include an improvement in the mode structure. This simulation is very useful, since one cannot determine the $p$ mode structure experimentally by simply analyzing the intensity pattern, e.g., differences in the petal structures of $\omega_{i}$, in Fig. 2(c,d), are not discernible, even though (c) has a contaminated mode structure.

III.2 Spontaneous Four-Wave Mixing

Next, we present a simulation for non-degenerate four wave mixing. In this scheme, two pump beams spontaneously create signal and conjugate beams, according to the polarization $\mathcal{P}(\omega_{c})^{(3)}=\epsilon_{0}\chi^{(3)}\mathcal{E}(\omega_{p1})\mathcal{E}(\omega_{p2})\mathcal{E}(\omega_{s})^{*}$; depicted in Fig. 3(a). In our other simulations, all of the fields which generate the new beam were carefully selected input beams. This last simulation is distinct since the signal and conjugate beams are spontaneously created during the interaction. Therefore, to predict the mode structure of the conjugate beam, one must make assumptions about the spontaneous response creating the signal beam. The question remains, how do we predict the coupling of the atoms to the vacuum LG modes at the signal and conjugate wavelengths?

To investigate this problem, we simulate an interaction, which has been experimentally realized (Walker et al., 2012). A $780\,\mathrm{nm}$ pump and a $776\,\mathrm{nm}$ pump stimulate a nonlinear interaction in Rb, which is the source of a $5230\,\mathrm{nm}$ signal beam and a $420\,\mathrm{nm}$ conjugate beam [see Fig. 3(a)]. In the experiment, the long wavelength beam remains unobserved. However unlikely, OAM transfer to this beam should not be ruled out as a possibility. Rather, the atoms can couple to any spatial light mode, which obeys OAM conservation, and thus we should account for each of these modes in our simulation. Therefore, for the sake of simulation, we input a signal beam that is a balanced incoherent superposition of all OAM modes, which obey OAM conservation, and allow the integral in Eq. (14) to determine which modes are populated. Thus, we accommodate for the possible OAM transfer to the long-wavelength beam, and only restrict our analysis to incoherent superposition at the signal wavelength. The last assumption we will impose for this simulation is the matching of the Rayleigh ranges for all beams, i.e., we invoke the Boyd criterion (Walker et al., 2012). In effect, we judiciously choose the widths of the beams to maximize the nonliniear interaction.

The results are presented in Fig. 3(b–f). From left to right, the columns correspond to the profile of the pumps, a histogram giving the mode structure of the conjugate beam, the profile of the conjugate beam, and lastly the experimental data from Ref. (Walker et al., 2012). We see excellent qualitative agreement, when the signal beam is allowed to carry OAM. In particular, our simulation appears to take account of the relative brightness and shape of the lobes in the experimental data. In each case there are lobes that are elongated and dimmer than others and, our simulation agrees with this observation. It is apparent that the richer mode structure that we predict for the conjugate beam can account for the variations in the lobe brightness, without drastically affecting the lobe structure. Therefore, we emphasize that a more careful analysis, such as the one we present here, is necessary when studying the mode structure of output beams, since a naive examination of the lobe structure can be very misleading.

IV Predictive Simulations

In the previous examples, we have shown how our theory predicts the mode structure during OAM transfer. Now we will show how it can be used as a tool to optimize the mode structure for OAM addition, subtraction, and cancellation processes. The theme of these simulations is to develop methods which enhance the mode structure of the generated beams. They also serve as predictions which can be verified experimentally, relatively easily.

IV.1 Second Harmonic Generation

First we present a simulation for second harmonic generation. An $1140\,\mathrm{nm}$ pump beam, with waist $w_{0}=0.1$ mm, interacts with a $3\,\mathrm{mm}$ long crystal to create a second harmonic beam at $570\,\mathrm{nm}$. The waist of the basis set, used to represent the generated beam, is set to match the waist of the pump beams. The light-matter interaction in the crystal is governed by the polarization $\mathcal{P}(2\omega)^{(2)}=\epsilon_{0}\chi^{(2)}\mathcal{E}(\omega_{1})\mathcal{E}(\omega_{2})$, depicted in Fig. 4(a).

We first consider that two photons are annihilated from a single pump beam in the superposition $u_{1,0}+u_{-1,0}$, and the results are given in Fig. 4(b). Each histogram depicts the probabilities $\mathscr{P}_{\ell,p}$ of modes being excited, i.e., referring to Eq. (14), $\mathscr{P}_{\ell,p}\equiv|\int dz^{\prime}c_{\ell,p}(z^{\prime})|^{2}/\sum_{\ell,p}|c_{\ell,p}|^{2}$. One can see that along with the response at $\ell=\pm 2$, there is also response at $\ell=0$, corresponding the cross terms in $(u_{1,0}+u_{-1,0})^{2}$. One can enhance the OAM transfer by creating a destructive interference to remove the response at $\ell=0$. This can be done by including a second pump beam, according to the geometry in Fig. 4(c), with a $\pi$ phase shift in the superposition, i.e., a rotated profile. The conservation of linear momentum dictates that, in order to have a response at $2\omega$ in the geometry of Fig. 4(c), a photon must be annihilated from each pump beam. As one can see in Fig. 4(d), this geometry (along with the phase shift) suppresses the response at $\ell=0$ and provides a cleaner mode structure. This technique can be quite useful and can be implemented in more complex situations, as we will show in the following examples.

IV.2 Third Harmonic Generation

Next, we present a simulation for third-harmonic generation. An $1140\,\mathrm{nm}$ pump beam, with waist $w_{0}=0.1\,\mathrm{mm}$, interacts with a $3\,\mathrm{mm}$ long crystal to create a third-harmonic beam at $380\,\mathrm{nm}$. The light-matter interaction in the crystal is governed by the polarization $\mathcal{P}(3\omega)^{(2)}=\epsilon_{0}\chi^{(2)}\mathcal{E}(\omega_{1})\mathcal{E}(\omega_{2})\mathcal{E}(\omega_{3})$, depicted in Fig. 5(a). As in the previous simulation, we first consider that three pump photons are annihilated from a single pump beam in the superposition $u_{1,0}+u_{-1,0}$, and the results are given in Fig. 5(b). One can see from the histogram that, along with the response at $\ell=\pm 3$, there is also response at $\ell=\pm 1$, corresponding the cross terms in $(u_{1,0}+u_{-1,0})^{3}$. To clean up the OAM transfer, one needs to create an interference to destroy the response at $\ell=\pm 1$. We deploy the same tactic as in the SHG simulation, but in this case include three pump beams, two of which have the phase shift $\pm 2\pi/3$ between the two modes. One can see in Fig. 5(c,d) that if one photon is annihilated from each of the three pump beams, that the response at $\ell=\pm 1$ is destroyed and a clean OAM transfer is established. Again, as we described in the previous subsection, one could realize this result by judiciously choosing the angles of the incident beams according to conservation of linear momentum, such that the generated beam is aligned along the $z\mathrm{-axis}$. So we see that this can be a very useful technique for the up-conversion of OAM, especially in the absence of a cavity, which would naturally clean the mode structure.

IV.3 Four-wave Mixing

Next, we present a simulation for degenerate four-wave mixing, in which we demonstrate the addition, subtraction, and cancellation of OAM. Two pump beams and a signal interact in a 7 cm long nonlinear cell to create a conjugate beam according to the polarization $\mathcal{P}(\omega_{c})^{(3)}=\epsilon_{0}\chi^{(3)}\mathcal{E}(\omega_{p})\mathcal{E}(\omega_{p})\mathcal{E}(\omega_{s})^{*}$, depicted in Fig. 6(a). We now investigate ways to tailor the response of the conjugate beam. In both cases the pump beams remain in the form $u_{1,0}+u_{-2,0}$ and $u_{1,0}+e^{i\pi}u_{-2,0}$, and we choose to vary the signal beam. We showed in the first three simulations that, when the pump beams are properly rotated with respect to each other, that cross terms can be canceled, i.e., response at $\ell=-1$ is suppressed. Therefore, if the seed beam were simply a $u_{0,0}$, then the conjugate beam would respond at $\ell=2$ and $\ell=-4$, as we see in Fig. 6(b). However, suppose we would prefer the field to respond at $\ell=\pm 3$, then we would choose the signal to be in a $u_{-1,0}$ mode and effectively add another unit of angular momentum to the conjugate beam, as seen in Fig. 6(c). This approach can be taken to the extreme by choosing the seed to be in a $u_{2,0}$ mode and effectively subtracting two units of angular momentum from the conjugate beam, as seen in Fig. 6(d). In this case, response at $\ell=0$ and $\ell=-6$ are the two modes allowed by OAM conservation. However we see that response at $\ell=-6$ is naturally suppressed due to the poor overlap with the pump and signal beams. Thus, we expose one of possibly many ways to completely suppress OAM transfer into unwanted modes.

V Conclusion

In conclusion, we report a general theory for calculating the Laguerre-Gauss mode structure of beams generated during nonlinear interactions. We accomplish this by making a first order Born-like approximation to the inhomogeneous paraxial wave equation. Therefore, it is akin to the weak scattering problem of the nonlinear Schr$\ddot{\mathrm{o}}$dinger equation, and the solution is consistent with the intuitive picture that the probability of a mode being excited is proportional to the overlap of the modes interacting in the nonlinearity. The theory is general in the sense that, it may be implemented for any complete orthonormal set of spatial mode functions. We use the theory to simulate orbital angular momentum transfer in several nonlinear optical processes, with an emphasis on analyzing and tailoring the resulting mode structure for optimal performance. This includes the processes of orbital angular momentum addition, subtraction, and cancellation in harmonic generation, parametric down conversion, and four-wave mixing.

Although this theory is based on a first order approximation, we show that it can model the mode structure of experimental beam data remarkably well. In doing so we are able to show that a naive analysis of the lobe structure in intensity patterns can be very misleading. One may be tempted to assume based on lobe structure that OAM transfer is limited to certain pathways. However, our analysis shows that the rich mode structure resulting from a more complete consideration of all the pathways, actually accounts for symmetries in the resulting beam pattern. On a final note, this validation of the theory is assurance that the theory is accurately predicting the classical mode structure, which can in turn be used for second-quantized treatments of these interactions.

This research was supported by the Air Force Office of Scientific Research grant FA9550-13-1-0098. In addition, RNL, ZX and JDP acknowledge additional support from the Army Research Office, National Science Foundation, and Northrop Grumman Corporation. The authors would also like to thank Sonja Franke-Arnold and her group for sharing the data appearing in Fig. 3.

Appendix A

A.1 Green Function Solution Method

In Sec. II, we introduced our spatial mode propagation equation and constructed an IVP to restate the problem in a compact form. The first step in solving our inhomogeneous IVP [Eq. (4)] is to solve the homogeneous IVP in free space, i.e., for $\wp=0$:

[TABLE]

Using a Green function method Barton (1989), we search for a propagator $K$ defined by

[TABLE]

such that when $K$ is known, the homogeneous problem is solved:

[TABLE]

The $z=z^{\prime}$ restraint is a standard condition placed by the boundary conditions which define the propagator, i.e., Eq (16).

Proof.

We check our solution Eq. (17) by operating the differential operator $\hat{L}$ and find

[TABLE]

The second line follows trivially since the propagator $K$ obeys Eq. (16a). Thus we see $\mathcal{E}_{\mathrm{hom}}$ is a solution to Eq. (15). ∎

To proceed, we recognize the property in Eq. (10) and explicitly define our propagator in terms of the LG modes:

[TABLE]

Inserting Eq. (19) into Eq. (16), we verify our choice of the propagator:

[TABLE]

Therefore, the homogeneous problem is solved, and we can turn our attention back to the inhomogeneous problem Eq. (4).

To solve the inhomogeneous problem, again in free space, we need a Green function $G$ defined by

[TABLE]

such that when $G$ is known, the inhomogeneous problem is solved:

[TABLE]

Proof.

We check our solution Eq. (22) by operating the differential operator $\hat{L}$ and find

[TABLE]

The second line follows from our definition of the Green function, Eq. (21a), and the third line follows trivially from the definition of the Dirac delta function. Thus we see $\mathcal{E}$ is a solution to Eq. (4). ∎

Building on the method for solving the homogeneous problem, we explicitly define our Green function $G$ in terms of the LG modes

[TABLE]

where $\Theta(z-z^{\prime})$ is the Heaviside step function. Next, inserting Eq. (24) into Eq. (21), we verify our choice of Green function:

[TABLE]

and thus we see that Eq. (21a) is satisfied, secondly, Eq. (21b) is satisfied since we know $G|_{z<z^{\prime}}=0$, by definition of the Heaviside function, and lastly, Eq. (21c) is satisfied since $r\rightarrow 0\Longrightarrow u_{l,p}\rightarrow 0\Longrightarrow G\rightarrow 0$. Therefore, we have derived a valid propagator and Green function and can now write the final solution by combining our two previous solutions Eq. (17) and Eq. (22). However, we would first like to simplify the notation in Eq. (22). It is likely that the source $\wp$ doesn’t contribute until some position $z^{\prime}=z_{i}$, and for $z>z^{\prime}=z_{i}$, the Heaviside function reduces to $1$, leaving only the propagator $K$. Furthermore, the source likely only contributes up to some position $z_{f}$ and thus we may safely modify the $dz$ integral. With these final considerations, we can write the complete solution to Eq. (4):

[TABLE]

Although the problem is solved, we can further employ the LG modes to simplify our calculations. Suppose we expand the nonlinear source in terms of the LG modes, i.e.,

[TABLE]

where $c_{l,p}(z)=\int rdrd\phi\,u_{l,p}^{*}(r,\phi,z)\wp(r,\phi,z)$. One should notice that, in contrast to expanding an arbitrary Gaussian beam in terms of the LG modes, the expansion coefficients for a nonlinear source distribution are not independent of the position $z$. Thus, each infinitesimal slice of the source contributes to the new field $\mathcal{E}$ in the following way. We insert Eq. (27) into Eq. (26) and invoke the orthogonality of the LG modes one last time to arrive at

[TABLE]

This result corresponds to Eq. (14) from the main text.

Lastly, we emphasize the $z$-dependence of the $c_{l,p}$ coefficients, which represent the amplitudes of the new mode structure of the beam. On one hand, these coefficients allow one to study how the new beam evolves during the interaction in the nonlinear medium. On the other, in the case that the first Born approximation does not hold, one can use this theory for successive iterations. For example, if the pump beam is known to be modified by the nonlinear material, then one can model this by a propagation equation for the pump, and use the solution as a second Born approximation for the new frequency components of the field, i.e.,

[TABLE]

where in this case $c_{l,p}$ are the result of a propagation equation for the pump.

A.2 Noncollinear Beam Geometry

We restrict our problem to only consider beams which are separated by some angle $\theta$, yet lie in the same plane and have focal points which coincide. One beam, typically the generated beam, will define the reference frame and the coordinates of the rotated beam will be transformed into this frame (see the geometry in Fig. 7).

The electric field at $P$ is described by a slice of the reference beam at $z$, and a slice of the tilted beam at $z^{\prime}$. One can easily deduce the form of $r^{\prime},\phi^{\prime},\mathrm{and}\,z^{\prime}$, in terms of $r,\phi,z,$ and $\theta$, using this geometry. However, the transformations make for cumbersome calculations. Therefore, since $\theta$ is small in our simulations, we expand around $\theta$ and keep only the first-order correction for each transformed coordinate:

[TABLE]

These corrections made no significant changes in our calculations, for which all beams interact at $\theta\leq 2$ degrees.