Self-focusing of multiple interacting Laguerre-Gauss beams in Kerr media

TL;DR

This paper derives the self-focusing critical power for multiple interacting Laguerre-Gauss beams in Kerr media, revealing complex phenomena and conditions where interactions alter focusing behavior compared to single beams.

Contribution

It generalizes the critical power concept to superpositions of Laguerre-Gauss beams, accounting for interactions and revealing new self-focusing phenomena.

Findings

Interactions can increase or decrease the critical power needed for self-focusing.

High-orbital angular momentum modes can focus with less power when interacting with low-angular momentum beams.

Stationary equilibrium configurations may not exist in superpositions, unlike single-beam cases.

Abstract

Using a variational approach, we obtain the self-focusing critical power for a single and for any number of interacting Laguerre-Gauss beams propagating in a Kerr nonlinear optical medium. As is known, the critical power for freely propagating higher-order modes is always greater than that of the fundamental Gaussian mode. Here, we generalize that result for an arbitrary incoherent superposition of Laguerre-Gauss beams, adding interactions between them. This leads to a vast and rich spectrum of self-focusing phenomena, which is absent in the single-beam case. Specifically, we find that interactions between different modes may increase or decrease the required critical power relative to the sum of individual powers. In particular, high-orbital angular momentum modes can be focused with less power in the presence of low-orbital angular momentum beams than when propagating alone. The decrease in required critical power can be made arbitrarily large by choosing the appropriate combinations of modes. Besides, in the presence of interactions, an equilibrium configuration of stationary spot-size for all modes in a superposition may not even exist, a fundamental difference to the single-beam case in which a critical power for self-focusing always exists.

Equations

$$\left(2i{k}_0\frac{\partial }{\partial z}+{\nabla }_{\perp }^2+2{\kappa }^2\text{absolutevalue}{a}^2\right)a=0,$$

$$\frac{1}{2}a({\mathbf{r}}_{\perp },z,t)\exp -i{k}_0(ct-z)+c.c.=\frac{e}{m{c}^2}A({\mathbf{r}}_{\perp },z,t),$$

$$\mathcal{L}=i{k}_0\left(a\frac{\partial a^{\ast }}{\partial z}-a^{\ast }\frac{\partial a}{\partial z}\right)+{\nabla }_{\perp }{a}^{\ast }\cdot {\nabla }_{\perp }a-{\kappa }^2{a}^2{a^{\ast }}^2,$$

$$\frac{\partial }{\partial z}\left(\frac{\partial \mathcal{L}}{\partial (\partial a^{\ast }/\partial z)}\right)+{\nabla }_{\perp }\cdot \left(\frac{\partial \mathcal{L}}{\partial ({\nabla }_{\perp }{a}^{\ast })}\right)-\frac{\partial \mathcal{L}}{\partial a^{\ast }}=0.$$

$$\frac{d}{dz}\left(\frac{\partial L}{\partial \dot{{\beta }_i}}\right)-\frac{\partial L}{\partial {\beta }_i}=0,$$

$$\begin{array}{rl}a=& A{C}_{p\ell }{\left(\frac{\sqrt{2}r}{W}\right)}^{\text{absolutevalue}\ell }{L}_p^{\text{absolutevalue}\ell }\left(\frac{2r^2}{W^2}\right)\exp -\frac{r^2}{W^2}\times \\ \times & \exp i\left(\ell \phi +k_0\frac{r^2}{2R}-\psi \right),\end{array}$$

$$\begin{split}L=&\frac{\pi}{2}A^{2}C_{p\ell}^{2}W^{2}\Biggl{[}-2I_{100}k_{0}\dot{\psi}+\frac{1}{2}I_{101}k_{0}^{2}W^{2}\left(\frac{1}{R^{2}}-\frac{\dot{R}}{R^{2}}\right)+\\ +&\frac{4\ell^{2}}{W^{2}}I_{10-1}-\frac{4\absolutevalue{\ell}}{W^{2}}(I_{100}+2I_{110})+\\ +&\frac{8}{W^{2}}\left(\frac{1}{4}I_{101}+I_{111}+I_{121}\right)-\kappa^{2}I_{200}A^{2}C_{p\ell}^{2}\Biggr{]},\end{split}$$

$$I_{mns}(p,\ell )={\int }_0^{\infty }{e}^{-mx}{x}^{m\text{absolutevalue}\ell +s}{\left[L_p^{\text{absolutevalue}\ell }(x)\right]}^{2m-n}{\left[L_{p-1}^{\text{absolutevalue}\ell +1}(x)\right]}^{n}dx.$$

$$\begin{array}{rl}& I_{100}=\frac{(p+\text{absolutevalue}\ell )!}{p!},I_{101}=\frac{(p+\text{absolutevalue}\ell )!}{p!}(2p+\text{absolutevalue}\ell +1),\\ & I_{10-1}=\frac{1}{\text{absolutevalue}\ell }\frac{(p+\text{absolutevalue}\ell )!}{p!},I_{110}=0,\\ & I_{121}=-I_{111}=\frac{(p+\text{absolutevalue}\ell )!}{(p-1)!},\\ & I_{200}=\frac{1}{2^{4p+2\text{absolutevalue}\ell +1}}{\left[\frac{(p+\text{absolutevalue}\ell )!}{p!}\right]}^2{S}_{p\ell },\\ & S_{p\ell }=\sum _{n=0}^p\frac{(2n)!{\left[(2p-2n)!\right]}^2(2\text{absolutevalue}\ell +2n)!}{(n!{)}^2{\left[(p-n)!\right]}^4{\left[(\text{absolutevalue}\ell +n)!\right]}^2}.\end{array}$$

$$\frac{d}{dz}\left(\frac{\pi }{2}{W}^2{A}^2{C}_{p\ell }^2{I}_{100}\right)=0⟹\frac{dP}{dz}=0.$$

$$\ddot{W}+\frac{4}{k_0^2{W}^3}\left[\frac{{\kappa }^{2}P}{4}\frac{2I_{200}{C}_{p\ell }^2}{I_{101}}-1\right]=0.$$

$$\begin{array}{rl}{P}_c& =\frac{4}{{\kappa }^2}\frac{I_{101}}{2I_{200}{C}_{p\ell }^2}=P_G\frac{1}{2I_{200}}{\left[\frac{(p+\text{absolutevalue}\ell )!}{p!}\right]}^2(2p+\text{absolutevalue}\ell +1)\\ & =P_G{4}^{2p+\text{absolutevalue}\ell }(2p+\text{absolutevalue}\ell +1)S_{p\ell }^{-1},\end{array}$$

$$P_c=P_G{4}^{\text{absolutevalue}\ell }\frac{\text{absolutevalue}\ell !(\text{absolutevalue}\ell +1)!}{(2\text{absolutevalue}\ell )!}.$$

$${\mathcal{L}}_{\mathrm{int}}=-{\kappa }^2\sum _{i,j\ne i}{a}_i^{\ast }{a}_i{a}_j^{\ast }{a}_j=-{\kappa }^2\sum _{i,j\ne i}\text{absolutevalue}{a_i}^2\text{absolutevalue}{a_j}^2.$$

$$\begin{split}\frac{\partial L_{\mathrm{int}}}{\partial W_{i}}&\bigg{\rvert}_{W_{i}=W}=\frac{4\kappa^{2}}{W_{i}^{3}}\left(\frac{\pi}{2}P_{i}C_{p_{i}\ell_{i}}^{2}\right)\\ &\times\sum_{j\neq i}P_{j}C_{p_{j}\ell_{j}}^{2}(\absolutevalue{\ell_{i}}H^{ij}_{00}-H^{ij}_{01}-2H^{ij}_{11})\,,\end{split}$$

$$\begin{array}{r}{H}_{ns}^{ij}(p_i,{\ell }_i,p_j,{\ell }_j)={\int }_0^{\infty }{e}^{-2x}{x}^{\text{absolutevalue}{\ell }_i+\text{absolutevalue}{\ell }_j+s}{\left[L_{p_i}^{\text{absolutevalue}{\ell }_i}(x)\right]}^{2-n}\\ \times {\left[L_{p_i-1}^{\text{absolutevalue}{\ell }_i+1}(x)\right]}^n{\left[L_{p_j}^{\text{absolutevalue}{\ell }_j}(x)\right]}^{2}dx.\end{array}$$

$$\begin{split}\ddot{W_{i}}&+\frac{4}{k_{0}^{2}W_{i}^{3}}\Bigg{[}\frac{\kappa^{2}P_{i}}{4}\frac{2I^{i}_{200}C_{p_{i}\ell_{i}}^{2}}{I^{i}_{101}}\\ &+\kappa^{2}\sum_{j\neq i}P_{j}C_{p_{j}\ell_{j}}^{2}\frac{\absolutevalue{\ell_{i}}H^{ij}_{00}-H^{ij}_{01}-2H^{ij}_{11}}{I_{101}^{i}}-1\Bigg{]}=0\,.\end{split}$$

$$\sum _j\left[{\delta }_{ij}+(1-{\delta }_{ij})2\frac{C_{p_j{\ell }_j}^2}{C_{p_i{\ell }_i}^2}\frac{\text{absolutevalue}{\ell }_i{H}_{00}^{ij}-H_{01}^{ij}-2H_{11}^{ij}}{I_{200}^i}\right]P_j^c=P_i^0,$$

$$\ddot{W_i}+\frac{4}{k_0^2{W}_i^3}\frac{{\kappa }^2}{4}\frac{2I_{200}^i{C}_{p_i{\ell }_i}^2}{I_{101}^i}\sum _j{M}_{ij}\delta P_j=0.$$

$$\begin{array}{rl}a=& A{C}_{mn}{H}_m\left(\frac{\sqrt{2}x}{W}\right)H_n\left(\frac{\sqrt{2}y}{W}\right)\exp -\frac{x^2+y^2}{W^2}\times \\ \times & \exp i\left(k_0\frac{x^2+y^2}{2R}-\psi \right),\end{array}$$

$$\begin{split}L=&\frac{1}{2}A^{2}C_{mn}^{2}W^{2}\Biggl{[}-2J^{m}_{100}J^{n}_{100}k_{0}\dot{\psi}\\ +&\frac{1}{2}(J^{m}_{102}J^{n}_{100}+J^{m}_{100}J^{n}_{102})k_{0}^{2}W^{2}\left(\frac{1}{R^{2}}-\frac{\dot{R}}{R^{2}}\right)+\\ +&\frac{8}{W^{2}}(m^{2}J^{m}_{120}J^{n}_{100}+n^{2}J^{m}_{100}J^{n}_{120}-mJ^{m}_{111}J^{n}_{100}-nJ^{m}_{100}J^{n}_{111})\\ +&\frac{2}{W^{2}}(J^{m}_{102}J^{n}_{100}+J^{m}_{100}J^{n}_{102})+-\kappa^{2}J^{m}_{200}J^{n}_{200}A^{2}C_{mn}^{2}\Biggr{]},\end{split}$$

$$J_{\alpha \beta \gamma }^q={\int }_{-\infty }^{+\infty }{e}^{-\alpha {\xi }^2}{\xi }^{\gamma }{\left[H_q(\xi )\right]}^{2\alpha -\beta }{\left[H_{q-1}(\xi )\right]}^{\beta }d\xi .$$

$$\begin{array}{r}{J}_{100}^q=\sqrt{\pi }{2}^{q}q!,J_{102}^q=\left(q+\frac{1}{2}\right)J_{100}^q,\\ J_{120}^q=\frac{1}{2q}{J}_{100}^q,J_{111}^q=\frac{1}{2}{J}_{100}^1.\end{array}$$

$$P_c=P_G\frac{\pi }{2J_{200}^m{J}_{200}^n}(2^{m+n}m!n!{)}^2(m+n+1).$$

$$\begin{split}\ddot{W}+\frac{4}{k_{0}W^{3}}\Biggl{[}&\frac{\kappa^{2}P_{i}}{8}\frac{2C_{n_{i}m_{i}}^{2}J_{200}^{m_{i}}J_{200}^{n_{i}}}{J_{102}^{m_{i}}J_{100}^{n_{i}}+J_{100}^{m_{i}}J_{102}^{n_{i}}}+\frac{k_{p}^{2}}{8}\sum_{j\neq i}P_{j}C_{m_{j}n_{j}}^{2}\\ &\times\frac{2m_{i}K_{11}^{m_{i}m_{j}}K_{00}^{n_{i}n_{j}}+2n_{i}K_{00}^{m_{i}m_{j}}K_{11}^{n_{i}n_{j}}-K_{02}^{m_{i}m_{j}}K_{00}^{n_{i}n_{j}}-K_{00}^{m_{i}m_{j}}K_{02}^{n_{i}n_{j}}}{J_{102}^{m_{i}}J_{100}^{n_{i}}+J_{100}^{m_{i}}J_{102}^{n_{i}}}-1\Biggr{]}=0\,,\end{split}$$

$$K_{\beta \gamma }^{qr}={\int }_{-\infty }^{+\infty }{e}^{-2{\xi }^2}{\xi }^{\gamma }{\left[H_q(\xi )\right]}^{2-\beta }{\left[H_{q-1}(\xi )\right]}^{\beta }{\left[H_r(\xi )\right]}^{2}d\xi .$$

$$\sum _j\left[{\delta }_{ij}+(1-{\delta }_{ij})2\frac{C_{m_j{n}_j}^2}{C_{m_i{n}_i}^2}\times \frac{2m_i{K}_{11}^{m_i{m}_j}{K}_{00}^{n_i{n}_j}+2n_i{K}_{00}^{m_i{m}_j}{K}_{11}^{n_i{n}_j}-K_{02}^{m_i{m}_j}{K}_{00}^{n_i{n}_j}-K_{00}^{m_i{m}_j}{K}_{02}^{n_i{n}_j}}{J_{200}^{m_i}{J}_{200}^{n_i}}\right]P_j=P_i^0,$$

$$\mathrm{Im}\left(a\frac{\partial a^{\ast }}{\partial z}\right)=\text{absolutevalue}{a}^2\left(\dot{\psi }+k_0\frac{r^2}{2R^2}\dot{R}\right).$$

$$\begin{split}\bm{\nabla}_{\perp}a^{*}&\cdot\bm{\nabla}_{\perp}a=\frac{\absolutevalue{a}^{2}}{r^{2}}\Biggl{[}\frac{k_{0}r^{4}}{R^{2}}+\ell^{2}+\\ &+\left(\absolutevalue{\ell}-\frac{2r^{2}}{W^{2}}\left(1+2\frac{L_{p-1}^{\absolutevalue{\ell}+1}(\frac{2r^{2}}{W^{2}})}{L_{p}^{\absolutevalue{\ell}}(\frac{2r^{2}}{W^{2}})}\right)\right)^{2}\Biggr{]}\,.\end{split}$$

$${\int }_0^{\infty }{e}^{-x}{x}^{\text{absolutevalue}\ell }{L}_p^{\text{absolutevalue}\ell }(x)L_q^{\text{absolutevalue}\ell }(x)dx=\frac{(p+\text{absolutevalue}\ell )!}{p!}{\delta }_{pq},$$

Full text

I Introduction

Lasers carrying orbital angular momentum (OAM) Allen et al. (1992) have attracted much attention recently due to a vast and interesting set of possible applications, including in plasma-based acceleration Vieira and Mendonça (2014); Vieira et al. (2016a), optical tweezers Padgett and Bowman (2011); Padgett (2014), quantum computation Mair et al. (2001); Molina-Terriza et al. (2007), super-resolution microscopy Vieira et al. (2016b), optical communications Wang et al. (2012); Bozinovic et al. (2013); Zhou et al. (2019), imaging Jack et al. (2009), and astrophysics Tamburini et al. (2011). Since many applications of OAM depend on the propagation of the beams in nonlinear media it is important to establish the self-interaction of OAM beams as well as the interactions between different OAM modes in those media.

One way to treat the laser-medium interaction is to consider the envelope evolution of the vector potential of the laser beam, which is described by some nonlinear partial differential equation, for instance, the nonlinear Schrödinger equation (NLSE) for the paraxial propagation of an ultra-intense short-pulse laser. The optical medium is described by a set of nonlinearities, which can be instantaneous (local along the beam), notably quadratic (Kerr type), or noninstantaneous, breaking translational invariance along the beam. An example of a medium with both types of nonlinearity is a relativistic underdense plasma, where an instantaneous quadratic nonlinearity arises from the relativistic transverse quiver motion of the electrons (relativistic mass correction) and a noninstantaneous one results from the coupling of the laser to the plasma waves. Neglecting noninstantaneous nonlinearities effectively reduces the problem to $(2+1)$-dimensional propagation (meaning that the dynamics along the beam decouple from the dynamics in a plane transverse to it; also, the coordinate along the beam is then essentially equivalent to time, and hence only one of them needs to be considered in addition to the two transverse coordinates).

Even with the above approximation, there are in general no exact solutions for the envelope of the beam. However, we might not even be interested in this exact solution, but instead in its dependence on a certain set of macroscopic parameters which have a clear physical interpretation or can be easily analyzed and controlled in an experiment. Examples of such parameters are the spot size of the beam, its centroid, or its phase. To this end, several methods have been developed and widely applied to the study of Gaussian beams, namely, the source-dependent expansion method Sprangle et al. (1987a, b); Esarey et al. (1994, 1997), the moment method Lam et al. (1975, 1975), fully numerical methods Sun et al. (1987), and the variational method Anderson and Bonnedal (1979); Duda et al. (1999); Duda and Mori (2000); Ren et al. (2001a, 2000, b, 2002). We shall employ the last in this Paper. It has been extensively used for Gaussian beams, in media both with and without noninstantaneous nonlinearities, leading to self-focusing, self-phase modulation, spot size self-modulation, and centroid hosing of a single beam Anderson and Bonnedal (1979); Duda et al. (1999); Duda and Mori (2000); Ren et al. (2001a), as well as the interaction of two beams Ren et al. (2000, 2001b, 2002); Dong et al. (2002); Wu et al. (2004), which leads to mutual attraction and spiralling, braiding, and merging. The interaction of any number of Gaussian beams in simple configurations has also been considered Ren et al. (2002).

The Gaussian beam is only the fundamental mode in the expansion of an arbitrary beam, and explicit expressions for the critical power for self-focusing of arbitrary higher-order transverse modes, free or interacting, are also of interest. In particular, we consider the decomposition of an arbitrary beam with circular symmetry into Laguerre-Gauss (LG) modes, which are especially relevant since they carry OAM. LG modes are characterized by two integers (see Section III below): the radial index $p\geq 0$ and the azimuthal index $\ell$, which is related to the vortex structure of the beam and directly quantifies its OAM.

The special case of $p=0$ consists of a radial profile with a single ring, and its critical power has been addressed in Refs. Kruglov et al. (1992); Vuong et al. (2006). Assuming the beam maintains an LG profile with OAM $\ell$ throughout the focusing process (what is called a self-similar collapse or aberrationless approximation) Kruglov et al. in Ref. Kruglov et al. (1992) directly integrated the NLSE and extracted the critical power from the conditions of existence of periodic solutions. Such a procedure is not easily generalized to profiles with higher $p$ (i.e. more radial nodes). By employing the variational method, Chen and Wang were able to obtain the critical power for propagation in a cubic-quintic medium Chen and Wang (2010). Although not explicitly written out, the critical power for an LG mode with arbitrary $p$ propagating in a Kerr medium can be obtained from the results of Ref. Chen and Wang (2010). The assumption of self-similar collapse is also built in the variational method (for a discussion see, e.g., Ref. Desaix et al. (1991)) and hence used throughout this Paper. In Refs. Fibich and Ilan (2000); Fibich and Gavish (2008a, b); Fibich (2015), it was shown that when the collapse is not self-similar the analytic prediction for the critical power given in Ref. Kruglov et al. (1992) becomes an upper bound for the critical power. For non-self-similar collapse, only for certain input vortex profiles can an analytic estimate be given for the critical power. Unfortunately, this is not the case for the LG modes.

Besides the analytical works already mentioned, the propagation and stability of LG modes and other types of vortex beams (e.g., Airy vortex beams Chen et al. (2010, 2013, 2015, 2016)) have also been extensively studied numerically, in instantaneous Kerr Kruglov and Vlasov (1985); Kruglov et al. (1988); Soljačić et al. (1998); Soljačić and Segev (2000, 2001); Grow et al. (2007); Ishaaya et al. (2008), nonlocal Buccoliero et al. (2007); Zhong and Yi (2007); Shen et al. (2011), and saturable Firth and Skryabin (1997); Skryabin and Firth (1998); Desyatnikov and Kivshar (2001); Bigelow et al. (2002); Desyatnikov and Kivshar (2002); Bouchard et al. (2016) media; see also Ref. Desyatnikov et al. (2005) for a review. In particular, Refs. Desyatnikov and Kivshar (2001); Bigelow et al. (2002); Bouchard et al. (2016) considered incoherent superpositions of two vortices.

In this Paper, we use the variational method to analytically study an arbitrary superposition of LG modes of any order, neglecting interference effects between them (this condition can be met, for example, when the phase of each mode varies arbitrarily). This analysis leads to a rich phenomenology of focusing phenomena, with the following two main novel findings. First, we show that an equilibrium configuration where all modes evolve with stationary spot size may not always exist, a fundamental difference regarding single-beam propagation where matched spot size evolution can always be attained. Second, we show that in some cases the total power required for self-focusing is lower for a set of interacting beams than for those beams propagating alone, which is the case of a high OAM mode being guided by low OAM modes. Since it is possible to fully sort LG modes Fu et al. (2018), the results of this work should be important for the guided propagation of intense pulses, with particular implications in compact laser-plasma accelerators and optical communications.

The Paper is organized as follows. In Section II we briefly review the variational method. As a steppingstone to interacting LG modes, we apply it to obtain the explicit formula for the self-focusing critical power of a single LG beam with arbitrary $p$ (Section III), relating it to known results in the literature. Section IV contains our central result, the critical power for an incoherent superposition of LG beams, Eq. (15), and a discussion of the rich associated phenomenology. We draw our conclusions and discuss possible applications and extensions to this work in Section V. For completeness, and as an illustration of the power of the variational method, in Appendix A we derive the critical power for arbitrary Hermite-Gauss modes (in a superposition or not), for which no known expressions exist in the literature, to the best of our knowledge. Appendices B and C present algebraic details and useful properties of special functions, respectively.

II Brief review of the variational method

Considering the instantaneous response of an optical medium with Kerr nonlinearity, the paraxial evolution of the envelope of a laser is described by the nonlinear Schrödinger equation

[TABLE]

where $a$ is the linearly polarized, normalized envelope of the vector potential $A$,

[TABLE]

$e$, $m$, $c$ are the electron charge, electron mass, and speed of light, respectively, $k_{0}$ is the laser wave number, $\kappa^{2}$ is a medium-dependent constant proportional to the nonlinear part of the refractive index (e.g., $\kappa^{2}\equiv k_{p}^{2}/8$ in a plasma, with $k_{p}^{2}$ the plasma wave number), $z$ is the coordinate along the beam, $\bm{r}_{\perp}$ are the coordinates transverse to $z$, and $\nabla^{2}_{\perp}$ is the Laplacian in the transverse plane. The paraxial wave equation is derived in speed of light frame variables $\tau\equiv z$, $\psi\equiv ct-z$. Yet, a single longitudinal variable $z$ suffices here, since we are assuming translational invariance along the beam. Some media, such as plasmas, have a minimum propagation (cutoff) frequency $\omega_{p}<\omega_{0}=k_{0}c$ (in the plasma, $\omega_{p}$ is the electron plasma frequency), and in such cases an extra term $-k_{p}^{2}a$ exists inside the bracket of Eq. (1). In this type of media, the envelope description (and hence this work) is valid only if the frequency of the beam is considerably larger than the cutoff frequency. However, the extra term in the NLSE just contributes with an overall constant in our Lagrangians below, and hence we simply drop it throughout.

Equation (1) can be obtained by minimizing the action $S=\int\mathcal{L}\,dzd\bm{r}_{\perp}$, where the appropriate Lagrangian density is

[TABLE]

using the Euler-Lagrange equations

[TABLE]

Instead of using the action in terms of the Lagrangian density $\mathcal{L}$ to solve the problem exactly, which is in general not possible, we seek an approximate solution. To that end, we can make an ansatz for the functional form of the envelope, introducing a set of parameters $\beta_{i}(z)$ depending only on $z$ which fully characterize the envelope $a=a(\beta_{i})$ and define a reduced Lagrangian $L(\beta_{i})=\int\mathcal{L}\,d\bm{r}_{\perp}$. The set of parameters act as generalized coordinates, with respect to which the action $S=\int L\,dz$ can be varied. The resulting reduced Euler-Lagrange equations,

[TABLE]

(with $\dot{\beta_{i}}=d\beta_{i}/dz)$, determine the evolution of the parameters, thus fully characterizing the evolution of the envelope $a$. By integrating out the transverse coordinates from the action we are thus substituting an infinite number of (transverse) degrees of freedom in the field $a$ by a finite number of mechanical coordinates $\beta_{i}$, reducing the problem to a system of coupled ordinary differential equations. In principle, the solution can be made arbitrarily exact by considering more parameters $\beta_{i}$, at the cost of computational complexity. The degree of approximation of the variational method depends, therefore, on the choice of trial function.

III Self-focusing of a single Laguerre-Gauss mode

We start by considering the $(p,\ell)$ Laguerre-Gauss mode as our trial function,

[TABLE]

where the parameters $\beta_{i}$ are the amplitude $A$, the spot size $W$, the radius of curvature $R$ and the phase $\psi$. The transverse coordinates are polar $(r,\varphi)$, with the beam centered at the origin, $p$ is the radial mode number (it gives the number of radial nodes), $\ell$ is the azimuthal mode number or OAM of the vortex, and $L_{p}^{\absolutevalue{\ell}}$ is an associated Laguerre polynomial. The normalization constant $C_{p\ell}=\sqrt{p!/(p+\absolutevalue{\ell})!}$ is such that the power $P$ of each mode is $P=(2/\pi)\int\absolutevalue{a}^{2}d\bm{r}_{\perp}=A^{2}W^{2}$ as is usually done for Gaussian beams. Other parameters could be considered, e.g., the centroid of the beam and some momenta transverse to the propagation axis (i.e. components $\bm{k}_{\perp}\perp\bm{k}_{0}$), but they are not crucial in what follows, since we consider only different beams centered at the origin as is appropriate for the modes resulting from a decomposition of an arbitrary beam. However, these parameters have to be considered if one wishes to study the interaction of several beams at different transverse locations.

Inserting the trial function into Eq. (2) and integrating over all $r$ and $\varphi$ we obtain the reduced Lagrangian (additional details are presented in Appendix B)

[TABLE]

where $I_{mns}$ is the following integral, which depends only on the beam mode:

[TABLE]

The evaluation of the relevant integrals (using the properties in Appendix C) yields

[TABLE]

We now apply the Euler-Lagrange equations to the Lagrangian of Eq. (4). Variation with respect to the phase $\psi$ gives rise to power conservation, since

[TABLE]

Since the variable $A$ arises in the Lagrangian only through the combination $A^{2}W^{2}$, we can replace it by the constant $P$, therefore avoiding variations with respect to $A$. Variation with respect to the radius of curvature $R$ relates $R$ to $W$ through $R=W/\dot{W}$. Both this auxiliary condition and power conservation are also present in the Gaussian case. Using these intermediate results, as well as Eqs. (6), the variation with respect to the spot size gives the equation for spot size dynamics:

[TABLE]

A stationary spot size is obtained when the term inside brackets vanishes, i.e. for $P=P_{c}$ where $P_{c}$ is the critical power given by

[TABLE]

where $P_{G}=4/\kappa^{2}$ is the critical power for a Gaussian beam, which is corrected by a factor characteristic of each LG mode. By setting $p=\ell=0$ in Eq. (9) we recover $P_{c}=P_{G}$ as we should. For the particular case where the laser intensity profile consists of a single ring ($p=0$), we find $I_{200}=4^{-\absolutevalue{\ell}}(2\absolutevalue{\ell})!/2$ and consequently recover the known result Kruglov et al. (1992):

[TABLE]

The expression for arbitrary $p$, Eq. (9), agrees with the result of Ref. Chen and Wang (2010), once model-specific constants are related and the quintic medium constant is set to zero.

The evolution of the critical power with OAM for various values of $p$ is presented in Fig. 1. It shows that the critical power rises monotonically with the OAM of the beams. The critical power grows with the radial number $p$. As a result, the mode which is most easily focused is the fundamental Gaussian mode. The fact that the critical power rises with $\ell$ can be understood heuristically using a simple physical picture, based on the centrifugal force felt by the photons of the OAM beam. The OAM beams have helical wavefronts, whereby photons at radius $r$ undergo azimuthal motion in the transverse plane with projected velocity $v/c=k_{\perp}/k_{0}$ where $k_{\perp}\sim-i\bm{\nabla}_{\perp}\sim\ell/r$ is transverse momentum. In their frame of reference, the photons are thus subjected to an outward centrifugal force $\absolutevalue{F}\sim k_{\perp}^{2}/r$. For a ring-shaped beam ($p=0$) the profile is peaked at $(r/W)^{2}=\absolutevalue{\ell}/2$ and we have $F\sim\sqrt{\absolutevalue{\ell}}$. Hence, the higher the OAM mode, the more the photons are pushed out, and the more difficult it is to focus them inwards, requiring a higher focusing power.

For systems with rectangular instead of circular symmetry, the expansion of a beam is best done in terms of Hermite-Gauss (HG) modes. Proceeding in the same way, one obtains the critical power for HG beams; see Appendix A.

IV Self-focusing of incoherent interacting Laguerre-Gauss modes

A general treatment of interacting higher-order modes must include the interference between them and is not readily treated by the variational method. We will consider only the incoherent case, where we neglect all interference between beams. This approximation is reasonable if the beams have (slightly) different frequencies or, for a sufficient number of different modes, random initial phases, whereby in both cases the interference would be averaged out during the propagation. If we consider only two beams, the incoherent approximation becomes exact if the beams have orthogonal polarizations. We decompose the beam envelope $a=\sum_{i}a_{i}$ into a linear combination of LG modes $a_{i}$, given by Eq. (3) with an index $i$ in the parameters $A$, $W$, $R$, and $\psi$ and mode numbers $(p_{i},\ell_{i})$. The amplitudes $A_{i}$ of each mode act as linear coefficients in the expansion. Then, the incoherent approximation leads to $\absolutevalue{a}^{2}\approx\sum_{i}\absolutevalue{a_{i}}^{2}$, which can be substituted in Eq. (1) for each mode envelope $a_{i}$. Hence, the interaction of incoherent beams is accounted by an extra term $2\kappa^{2}\sum_{j\neq i}\absolutevalue{a_{j}}^{2}a_{i}$ in the NLSE or equivalently by an extra pairwise amplitude-amplitude interaction term in the Lagrangian density Ren et al. (2002),

[TABLE]

We can again understand the form of the interaction in terms of a physical picture. A relativistic beam will increase locally the index of refraction $n$ of the medium according to $n=n_{0}+\Delta n$, with $n_{0}$ the linear index of refraction and $\Delta n\propto\kappa^{2}\absolutevalue{a}^{2}$. LG modes propagate in the $z$-direction, have a gradient of $n$ in the $r$ direction, and hence are attracted towards radii of higher $n$. When two beams interact, the peaks of the first attract the second and vice versa, while the nodes of intensity do not attract at all, in accordance with Eq. (11). In the case of an OAM beam, this effective attractive force has to compete with the centrifugal force discussed previously. This physical picture proves valuable in analyzing the results at the end of this section.

Since the new interaction term does not depend on the phases of the modes, only the equation for the spot sizes is modified. Setting $W_{i}=W$ for all modes after differentiating, the variation with respect to $W_{i}$ results in an extra term

[TABLE]

with the $H_{ns}$ integrals defined analogously to before:

[TABLE]

Adding Eq. (12) to Eq. (8), leads to a corrected equation of motion for the spot size of each mode given by

[TABLE]

Equating the term inside square brackets to zero yields the critical power for each beam, while the set of powers which leads to all modes evolving with constant $W_{i}$ is the solution to the following system of linear equations:

[TABLE]

where $P_{i}^{0}$ is the free critical power given by Eq. (9) and $\delta_{ij}$ is the Kronecker delta.

Equation (15) admits solutions where the waists of all beams propagate with constant spot size $W_{i}$. This kind of solution represents an equilibrium configuration where all spot sizes are stationary. In addition to these self-focusing solutions, Eq. (15) also predicts the existence of mode combinations where no stationary solution can be found. This is in stark contrast with the single-mode analysis that characterizes the usual self-focusing theory, valid for a single beam.

We start with the cases where there is no equilibrium. For some combinations of interacting modes, a number of $P^{c}_{j}$ come out negative, and since power is non-negative one concludes that no physical solutions exist and hence the equilibrium cannot be attained. In terms of the above physical picture, one interprets this as an impossibility to exactly balance all the attractions of peaks and the centrifugal forces on the beams, and to attain the equilibrium one would require repulsive forces (expressed as negative power), to stabilize these configurations. Examples of this situation are the interaction of three beams with $p=0$ and $\ell=0$, $\pm 1$ with the $(1,0)$-mode, in which all four required powers are negative, and the interaction of five beams with $p=0$ and $\ell=0$, $\pm 1$, $\pm 2$ in which only the $\ell=\pm 2$ would require negative power. Let us be more specific. For any of these cases, if one inserts any positive (that is, physical) set of values for powers $P_{i}$ into Eq. (14), then one obtains at least one $W_{i}>0$, i.e. at least one beam defocuses. By varying $P_{i}$ smoothly, all $\ddot{W}_{i}$ will also vary smoothly. But we have seen that one can never go through the equilibrium configuration $\ddot{W}_{i}=0$ for all $i$ (because the powers required for that laid outside the set of physical solutions). So, we cannot reach $\ddot{W}_{i}<0$ (all beams focus), because any smooth variation would have to go through $\ddot{W}_{i}=0$. We thus conclude that for a set of interacting higher-order modes, it may not be possible to reach self-focusing propagation simultaneously for all beams. This result is fundamentally different from the single-beam case, where self-focusing is always possible, given high enough power.

In the cases where the equilibrium configuration does indeed exist, it is not enough to raise slightly the power of one beam to focus it, since each power now depends on the values of all others, and changing one could break the equilibrium of the others. One has to define a perturbation of the power for all beams and reinsert it into Eq. (15). We denote the term inside square brackets in Eq. (15) by $M_{ij}$ and perturb $P_{j}=P^{c}_{j}+\delta\!P_{j}$. Then the $i$-th beam will focus if $\sum_{j}M_{ij}\,\delta\!P_{j}$ is positive and will defocus if it is negative (the limiting case $\sum_{j}M_{ij}\,\delta\!P_{j}=0$ corresponds to the equilibrium configuration). To see this, we note that, using Eq. (15) written as $\sum_{ij}M_{ij}(P_{j}^{c}+\delta\!P_{j})=P^{0}_{j}$, Eq. (14) can be rewritten as

[TABLE]

Since the prefactor of $\sum_{j}M_{ij}\,\delta\!P_{j}$ is always positive, focusing (resp. defocusing) of the $i$-th beam, i.e. $\ddot{W}_{i}<0$ (resp. $\ddot{W}_{i}>0$), occurs if $\sum_{j}M_{ij}\,\delta\!P_{j}>0$ (resp. $<0$). If we can find a vector $\delta\!P$ such that all of the entries of $M\cdot\delta\!P$ are positive, then it is possible to focus all of the beams simultaneously. This procedure is illustrated in several examples to follow.

A remarkable result is that the sum of the critical powers for all the beams can be less when they interact than when they do not. Also, for free modes it was seen that the Gaussian beam always has the lowest critical power. However, in the presence of interactions, this no longer holds. An example illustrative of both these results is the interaction of four modes with $p=0$ and $\ell=0$, $1$, $2$, $4$, whose profiles are depicted in Fig. 2, so that the physical picture can be visualized, together with their free and interacting critical powers. In this case, any perturbation around the equilibrium in which the powers of the $\ell=0$, $1$ beams are slightly higher than the critical power will lead to focusing of all four beams. Thus, it is possible to focus the interacting beams with only around two-thirds of the power required for focusing the four beams if they were propagating individually through the medium. Note also that, in this superposition, the Gaussian beam no longer has the lowest critical power, having indeed the highest.

Even more interesting is the possibility of using low OAM modes to help focus a high OAM beam. For instance, the mode with $p=0$ and $\ell=10$ has a (free) critical power $P_{c}\approx 62.43P_{G}$. However, if we consider the simultaneous propagation of this mode together with the three modes with $p=0$ and $\ell=0,5,8$, the total power required to focus all four beams at once is only $P_{c}\approx 37.87P_{G}$. Hence, it is energetically favorable to focus this set of beams instead of the $\ell=10$ beam individually, i.e. there is power saved by injecting further modes in the medium. In principle, by choosing modes with higher $p$ and $\ell$ and combinations with a higher number of modes, the decrease in critical power can be made arbitrarily large. However, the final required power would still be very high (although much less than the free power) since the critical power rises very fast with increasing mode. Nonetheless, it could be that, in some cases, the interactions bring the critical power from above the current technological capabilities to below them; the above, complicated, combinations may, then, prove to be important.

In some other cases, all of the critical powers increase. As an example, the interaction of the $(1,0)$, $(0,2)$, and $(0,3)$ modes, has a fivefold increase in critical power required as is depicted in Fig. 3. It is again qualitatively explained by the physical picture. As before, the increase of critical power can be made arbitrarily large. These cases are also of considerable importance when one wants the beams not to focus, and the calculations of critical power must include the interactions to ensure that the power is in fact below the threshold.

Given an arbitrary beam envelope with circular symmetry, one can decompose it into LG modes, and knowing the spot size of the beam, the coefficients of the expansion squared give the power of each mode, which can then be compared to the solutions of Eq. (15) to check if the beam as a whole will focus or not, as long as the modes can be considered incoherent superpositions. This analysis results from a situation where one solves Eq. (15) for all beams involved in the interaction. However, it is also possible to tune and fix the power of one (or more) of the modes to increase or decrease $P_{c}$ for the remaining ones. In this case, only a subsystem of Eq. (15) is to be solved.

V Conclusions

In summary, we have obtained explicit expressions for the critical self-focusing power of all basis modes with circular symmetry in a Kerr medium. Here, the main conclusion is that the fundamental Gaussian mode has the lowest critical power, which increases monotonically with both $p$ and $\ell$. In addition, we have shown that, as was already known for Gaussian beams Ren et al. (2000, 2001b, 2002), LG modes can attract each other, altering the power thresholds for simultaneously propagating beams and for different modes of a decomposition of an arbitrary beam. Although the interaction is always attractive, the existence of multiple peaks and troughs in the intensity profile of LG beams allows for combinations of modes where the interaction leads to focusing or defocusing of some or all beams. Together with the natural diffraction of the laser beams and the enhanced defocusing of LG modes due to its OAM, this allows for rich possibilities of focusing phenomena. In particular, in major contrast with single-beam propagation, propagation for all beams with stationary spot size may be impossible to attain.

Our results could prove useful in a variety of situations. First, if we want to perform experiments with multiple modes where self-focusing is important, one can perform the experiments simultaneously, saving power in the focusing process, since it is now possible to fully sort LG modes Fu et al. (2018). Second, even if one is interested in the propagation of a single beam (in particular, a high OAM one), it may be energetically favorable to use other beams (particularly low OAM ones) to help focus it. Third, when one wants the beams not to focus (for instance, to prevent optical damage to the nonlinear medium), secondary beams could be used to raise the critical power of the main beam, avoiding self-focusing.

The quantitative expressions derived in this work should be most reliable at the onset of self-focusing and for moderate powers (not too high compared with $P_{c}$). On the one hand, this ensures that higher-order corrections in the NLSE may be neglected and also that the paraxial approximation holds. On the other hand, it has been pointed out Chen and Wang (2010) that the results of the variational method are most accurate for moderate powers. This shortcoming at intense powers could, in principle, be overcome by choosing a trial function with more variational parameters, and therefore does not rule out the variational method itself as an accurate approach to self-focusing. Notwithstanding the (controllable) loss of quantitative accuracy in certain situations, the qualitative picture introduced in this work remains valuable.

Furthermore, the set of variational parameters considered in this Paper allows for the study of self-focusing only. The inclusion of centroid positions and transverse momenta would allow other possibilities, namely, hosing of a single beam or interaction of several beams at different transverse positions, as has been done for the Gaussian case. Furthermore, this work considered only instantaneous quadratic nonlinearities, neglecting spatiotemporal evolution along the beam (e.g., the coupling of a laser to plasma waves). The Lagrangian density for this more general case is known Duda and Mori (2000); Ren et al. (2001b, 2002), but a study of this phenomenon for higher-order modes is still lacking. Also of interest would be to extend the variational method to include interference between modes. In that case, power would no longer be constant and could be exchanged between the beams, and dynamical equations for the amplitude and the phase would need to be taken into account. The phenomenon of filamentation could possibly be treated this way.

Finally, other trial functions (say, Hermite-Gauss, Bessel, or Airy beams) could be employed to $(i)$ obtain the respective critical powers and $(ii)$ investigate some new possible phenomenology arising due to interactions. We followed path $(i)$ for Hermite-Gauss beams, and the expressions thus obtained are given in Appendix A.

Acknowledgements.

We acknowledge fruitful discussions with Luis Oliveira e Silva. We thank the anonymous reviewer whose comments have greatly helped improve this work. We acknowledge the EU Accelerator Research and Innovation for European Science and Society (EU ARIES) Grant Agreement No. 730871 (H2020-INFRAIA-2016-1). J. V. acknowledges the support of FCT (Portugal) Grant No. SFRH/IF/01635/2015.

Appendix A Self-Focusing of Hermite-Gauss Beams

For completeness, we now consider the case where a system exhibits rectangular instead of cylindrical symmetry. The beam expansion is then most conveniently done in terms of Hermite-Gauss (HG) modes. The trial function is the $(m,n)$ HG mode,

[TABLE]

where $A$, $W$, $R$, and $\psi$ are as above, the transverse coordinates are Cartesian $(x,y)$, $m$ and $n$ are the $x$ and $y$ mode numbers, respectively, and give the number of nodes in each direction, and $H_{n}$ is a Hermite polynomial. The normalization constant $C_{mn}=1/\sqrt{2^{m+n}m!n!}$ is again chosen such that $P=A^{2}W^{2}$, irrespective of the mode.

Following the same procedure as in Section III, the integrated Lagrangian is found to be

[TABLE]

where the relevant integrals for the HG beams are defined by

[TABLE]

The $\alpha=1$ integrals can be readily evaluated using the orthogonality and recursion relations for Hermite polynomials of Appendix C:

[TABLE]

Variation of the Lagrangian with respect to $\psi$ and $W$ gives the same intermediate results as in Section III and the variation with respect to the spot size yields the HG critical power

[TABLE]

As in the cylindrical case, the lowest mode ($m=n=0$) recovers the Gaussian result as expected and for higher modes self-focusing is more difficult since the power threshold rises.

The interaction Lagrangian density is again the one of Eq. (11), which after integration in the transverse plane and variation with respect to $W$ leads to the interacting equation of motion for the spot sizes of each beam,

[TABLE]

where the integrals $K_{\beta\gamma}^{qr}$ are defined as

[TABLE]

Equating the term in square brackets to zero gives the linear system

[TABLE]

where $P_{i}^{0}$ is given by Eq. (21), whose solution is the set of interacting critical powers for interacting HG modes.

Appendix B Algebraic computations leading to the reduced LG Lagrangian

In this appendix we present some auxiliary computations leading to the Lagrangian of Eq. (4). We start with the first term in the Lagrangian density of Eq. (2), noting that $a\,\partial a^{*}/\partial z-a^{*}\partial a/\partial z=2i\,\mathrm{Im}(a\,\partial a^{*}/\partial z)$. Using the trial function of Eq. (3), we have

[TABLE]

For the gradient term, we must use the derivative property of associated Laguerre polynomials $(L_{p}^{\absolutevalue{\ell}})^{\prime}(x)=-L_{p-1}^{\absolutevalue{\ell}+1}(x)$, where the prime denotes differentiation with respect to the whole argument. Using the gradient in polar coordinates, $\bm{\nabla}_{\perp}=\bm{\hat{r}}\partial/\partial r+\bm{\hat{\varphi}}\,1/r\,\partial/\partial\varphi$, we obtain

[TABLE]

The last term is trivial, $a^{2}a^{*2}=\absolutevalue{a}^{4}$. Concerning the integration of $\mathcal{L}$ in the transverse plane, all integrals in $\varphi$ give $2\pi$, and the radial integrals are proportional to either $\int_{0}^{\infty}\absolutevalue{a}^{2}rdr$, $\int_{0}^{\infty}\absolutevalue{a}^{4}rdr$, $\int_{0}^{\infty}\absolutevalue{a}^{2}r^{3}dr$, $\int_{0}^{\infty}\absolutevalue{a}^{2}/rdr$, or $\int_{0}^{\infty}\absolutevalue{a}^{2}\,L_{p-1}^{\absolutevalue{\ell}+1}(2r^{2}/W^{2})/L_{p}^{\absolutevalue{\ell}}(2r^{2}/W^{2})rdr$, which can be brought to the form $I_{mns}$ of Eq. (5) using the change of variable $x=2r^{2}/W^{2}$. Adding all the above terms yields the Lagrangian of Eq. (4).

The derivation of the HG Lagrangian is identical, but uses the derivative of the Hermite polynomial, $(H_{n})^{\prime}(x)=2nH_{n-1}(x)$. In this case, the integrals can be separated in $x$ and $y$, and hence every term in Lagrangian of Eq. (18) has two $J$-integrals multiplying it.

Appendix C Useful properties of Special Functions

The following orthogonality and recursion relations of the associated Laguerre polynomials Abramowitz and Stegun (1972) were used in evaluating the $I$-integrals:

[TABLE]

Also used were the orthogonality and recursion relations for Hermite polynomials, given by

[TABLE]