Dimensional reduction for optical beams with thermal nonlocal nonlinearity

Abstract

Nonlocal optical nonlinearities arising from the thermorefractive effect provide a long-range material response determined by heat diffusion and absorption. In graded-index media, this nonlocality fundamentally alters modal interactions, yet its accurate modeling remains computationally demanding when starting from the full spatial nonlinear Schrödinger equation. In this work, inspired by the nonpolynomial Schrödinger equation (NPSE) framework, we extend the dimensional reduction techniques to incorporate thermally mediated nonlocal nonlinearities. By coupling the optical field to an equation for the temperature-induced refractive index change, and employing a variational ansatz based on Laguerre–Gauss modes of the annular kind, of arbitrary azimuthal order, we derive explicit analytic expressions for the variational equations. The resulting effective model captures the dependence of the nonlinear interaction on mode order and degree of nonlocality, providing a tractable reduced description of the dynamics in thermal nonlocal media.

Graphical abstract

Keywords

nonlinear optics / beam propagation / spatial solitons

1 Introduction

Nonlocal optical nonlinearities represent an interesting self-interaction modality for waves propagating in media. They are based on the fact that the material optical response at a specific point is determined by the combined effect of the wave over an extended region of space [1]. This nonlocality drastically alters the dynamics of the optical field compared to local nonlinearity like the Kerr effect, often acting as a stabilizing mechanism for nonlinear waves. It has indeed been extensively predicted and experimentally demonstrated that nonlocal nonlinearities can arrest catastrophic wave collapse [2–4], suppress modulational instability in defocusing regimes [5], and support the formation of stable spatial solitary waves [6, 7]. Beyond stabilization, the nonlocal response profoundly affects soliton mobility, particularly in the presence of optical lattices [8]. Nonlocal interactions also arise in other systems, notably dipolar Bose–Einstein condensates [9, 10], where the long-range dipolar forces enable the emergence of supersolid states of matter [11].

Among the mechanisms generating optical nonlocality, thermal nonlinearities are especially important due to their ubiquity in absorptive media and their long-range character set by heat diffusion. Originating from light absorption and subsequent temperature redistribution [12–14], the thermal response is governed by a Poisson-like equation whose source is the absorbed optical power. Unlike rapidly decaying nonlocal responses, thermal nonlinearity is strongly influenced by geometry: boundary conditions shape the refractive-index profile and allow controlling soliton trajectories [15, 16]. This leads to a variety of nonlinear phenomena parametrized by the geometry of the medium [17–20]. Thermal diffusion can also sustain dispersive spatial shock waves [21, 22], offering a platform for studying extreme wave dynamics. A crucial feature of thermal nonlocality is its ability to stabilize nonlinear beams and spatial solitons. In local focusing cubic media, beams carrying orbital angular momentum undergo azimuthal instabilities and break into filaments [23]. Thermal diffusion suppresses this breakup, enabling stable vortex beams and ring solitons [24, 25]. This stabilization also influences vortex nucleation and the dynamics of vortex lines [26]. Parallel to these developments, interest in multimode optical fiber dynamics, and in particular in graded-index fibers (GRIN), has grown, motivated by the demand for higher data capacity in optical communications and the exploration of new nonlinear regimes [27, 28]. GRIN fibers with parabolic index profile have the additional advantage of being analytically treatable [29–32]. Recent work has shown that interacting Laguerre–Gauss beams in GRIN fibers exhibit modified self-focusing and distinct critical powers compared with single-mode excitation [33].

However, studying nonlocal thermal nonlinearities in a generic GRIN medium poses significant computational challenges, since fully spatial numerical simulations are costly, therefore, reduced dimensional models are desirable. Recently, we introduced a nonpolynomial Schrödinger equation (NPSE) [34, 35] for multimode fibers via a variational reduction of the full spatial nonlinear Schrödinger equation [36]. While effective for local Kerr nonlinearities, that model did not include nonlocal interactions. In this article, we discuss the extension of the NPSE framework to nonlocal thermal nonlinearities, deriving a coupled model consisting of a full spatial NLSE and a Poisson-type equation for the thermally induced index change. Using a variational ansatz based on Laguerre–Gauss modes with arbitrary azimuthal index and zero radial index [37], we obtain explicit expressions for the nonlocal contribution. We show that the standard NPSE ansatz is however insufficient to describe diffraction properly, and a radial phase term becomes necessary. The resulting effective one-dimensional model captures the dependence of the nonlinear interaction on mode order and transverse width, providing a generalized framework for thermally mediated interactions.

2 Fundamental equations

2.1 Spatial nonlinear Schrödinger equation in a graded-index medium

We consider the monochromatic electric field envelope \mathrm{\Phi }(x,y,z) propagating in an isotropic medium with refractive index n_0(x,y). Under the paraxial approximation, the propagation equation reads [30, 36]

(1)\begin{array}{rl}\mathrm{i}{\mathrm{\partial }}_z\mathrm{\Phi }=& -\mathrm{i}\frac{{\alpha }_0}{2}\mathrm{\Phi }-\frac{1}{2}{\mathrm{\nabla }}_{\perp }^2\mathrm{\Phi }\\ & +W(x,y)\mathrm{\Phi }+g{n}_{\mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{l}\mathrm{o}\mathrm{c}}(x,y,z)\mathrm{\Phi },\end{array}

where linear losses are modeled through the absorption coefficient {\alpha }_0; W(x,y) represents the GRIN potential,

(2)W(x,y)={\displaystyle \frac{1}{2}}{\mathrm{\Omega }}^2{r}^2,

with r^2=x^2+y^2, and the nonlinear index shift is written as

(3)n_{\mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{l}\mathrm{o}\mathrm{c}}(x,y,z)={\displaystyle \frac{\mathrm{\Delta }n(x,y,z)}{n_0}},

with n_0 a background refractive index. Unlike in Ref. [36], we consider a stationary field, neglecting the time dependence of the field \mathrm{\Phi }. The normalization units in Eq. (1) can be conveniently chosen as the beam waist w_0 on the transverse dimensions (x,y) and the corresponding Rayleigh range z_R=\pi w_0^2/{\lambda }_0 [38] on the axial dimension z, {\lambda }_0 being the wavelength associated to a bulk material with the refractive index present at the propagation axis.

2.2 Steady-state temperature field

The nonlinear shift in the refractive index (3) is given in our model by the thermo-optical effect, namely the refractive index change induced by a change of temperature of the optical medium. The change of temperature is itself induced by the local absorption, but due to thermal diffusion it will be nonlocal in space at the steady state. The field \mathcal{T}(x,y,z) describing temperature variation with respect to an ambient temperature, obeys the following steady-state heat equation [21, 22]

(4)-{\mathrm{\nabla }}^2\mathcal{T}=\eta |\mathrm{\Phi }{|}^2,

where \eta depends on the absorption {\alpha }_0 and the thermal diffusivity \kappa through

(5)\eta ={\displaystyle \frac{{\alpha }_0}{\rho c_p\kappa }},

with \rho the material density and c_p the specific heat at constant pressure. We will assume that the thermodynamic parameters \rho, c_p and \kappa are temperature independent, since we assume to have overall weak temperature variations, so we can use the parameters evaluated at a fixed temperature and pressure. The stationary thermal field is solved in the presence of the following boundary conditions: over the transverse directions x and y we have null boundary conditions at infinity, and over the longitudinal direction z we consider null boundary conditions at a flat input facet at z=0, and at a flat output facet at z=L. The refractive index change associated with the temperature variation is given by the local expression

(6)\mathrm{\Delta }n(x,y,z)=\beta \mathcal{T}(x,y,z),

where \beta =\mathrm{d}n/\mathrm{d}\mathcal{T} is the linear thermo-optic coefficient, which can be positive or negative, depending on the material.

In the presence of an optical beam propagating in the medium, \mathcal{T} varies slowly along z and the axial component of the Laplacian in Eq. (4) may be replaced by an effective confinement parameter 1/L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2 [1, 7, 19, 22, 39]. A convenient way to understand this approximation is to recall that the stationary temperature profile along z is expressed by the eigenmodes of the axial Laplacian, with the boundary conditions at z=0 and z=L (see also [39]). When analyzing the temperature profile distant from the edges, in the weak absorption regime, the full three-dimensional Green function expressed in Refs. [39, Eq. (13)], [[19], Eq. (5)] can be approximated so as to contain only the contribution from the fundamental mode along the z direction. The result of this approximation is that L_{\mathrm{e}\mathrm{f}\mathrm{f}}=L/\pi. The case of strong absorption can be handled using a similar approach but with a modified value of L_{\mathrm{e}\mathrm{f}\mathrm{f}} [21]. The resulting equation, valid sufficiently far from the sample boundaries, is therefore a 2D screened Poisson-type equation

(7)\left(-{\mathrm{\nabla }}_{\perp }^2+{\displaystyle \frac{1}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}}\right)\mathcal{T}(x,y,z)=\eta |\mathrm{\Phi }(x,y,z){|}^2.

The formal solution of Eq. (7) at fixed z is the transverse convolution

(8)\mathcal{T}(x,y,z)={\displaystyle \int {\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }\eta G({\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime })|\mathrm{\Phi }({\mathit{r}}_{\perp }^{\prime },z){|}^2,}

where G is the Green’s function, and {\mathit{r}}_{\perp }, {\mathit{r}}_{\perp }^{\prime } indicate transverse coordinates. In two transverse dimensions, this Green function is known in closed form and involves the modified Bessel function K_0 [24, 40]:

(9)G(\mathit{r})={\displaystyle \frac{1}{2\pi }}{K}_0\left({\displaystyle \frac{r}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\right).

At short distances K_0 reproduces the logarithmic singularity of the 2D Poisson kernel, while for r\gg L_{\mathrm{e}\mathrm{f}\mathrm{f}} it decays exponentially, reflecting the suppression of long-range diffusion by the heat transfer on the facets.

Within this model, the variation of the nonlinear index will satisfy

(10)g{n}_{\text{non-loc}}(\mathit{r},z)={\displaystyle \frac{\gamma }{2\pi }}{\displaystyle \int {\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }{K}_0\left({\displaystyle \frac{|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\right)|\mathrm{\Phi }({\mathit{r}}_{\perp }^{\prime },z){|}^2,}

with \gamma =g\beta {\alpha }_0/(\rho c_p\kappa ).

Combining the above elements, the optical propagation obeys the full spatial nonlocal NLSE

(11)\mathrm{i}{\mathrm{\partial }}_z\mathrm{\Phi }=-\mathrm{i}{\displaystyle \frac{{\alpha }_0}{2}}\mathrm{\Phi }-{\displaystyle \frac{1}{2}}{\mathrm{\nabla }}_{\perp }^2\mathrm{\Phi }+W(r)\mathrm{\Phi }+\gamma \mathrm{\Theta }[|\mathrm{\Phi }{|}^2]\mathrm{\Phi },

where the nonlocal operator is defined through the convolution

(12)\mathrm{\Theta }[|\mathrm{\Phi }{|}^2]({\mathit{r}}_{\perp },z)={\displaystyle \frac{1}{2\pi }}{\displaystyle \int {\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }{K}_0\left({\displaystyle \frac{|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\right)|\mathrm{\Phi }({\mathit{r}}_{\perp }^{\prime },z){|}^2.}

We remark that the 2D screened-Poisson model offers a substantial computational advantage over the full diffusion equation, as it integrates naturally into beam-propagation simulations where the nonlocal response is updated independently at each axial step. It also provides a convenient starting point for the variational dimensional reduction introduced below. Such a model is not unique to thermal media: screened-Poisson responses also describe reorientation nonlinearity in liquid crystals [19] and electrostrictive effects [7] and electron-ion interactions in plasmas [24]. In all these cases, the stationary transverse index profile is well approximated by a screened Poisson equation [21].

3 Effective action with thermal nonlocality

In order to perform a dimensional reduction, inspired by the NPSE formalism, in the presence of thermal nonlocality, it is convenient to rewrite Eq. (11) in variational form. The lossless version of Eq. (11) (i.e., setting {\alpha }_0=0), is the equation of motion associated with the action

(13)S={\displaystyle \int \mathrm{d}z{\displaystyle \int {\mathrm{d}}^2{\mathit{r}}_{\perp }\mathcal{L},}}

with Lagrangian density

(14)\begin{array}{rl}\mathcal{L}=& \frac{\mathrm{i}}{2}({\mathrm{\Phi }}^{\ast }{\mathrm{\partial }}_z\mathrm{\Phi }-\mathrm{\Phi }{\mathrm{\partial }}_z{\mathrm{\Phi }}^{\ast })-\frac{1}{2}|{\mathrm{\nabla }}_{\perp }\mathrm{\Phi }{|}^2-W(r)|\mathrm{\Phi }{|}^2\\ & -\frac{\gamma }{2}|\mathrm{\Phi }{|}^2\mathrm{\Theta }[|\mathrm{\Phi }{|}^2].\end{array}

The first line of Eq. (14) contains the local contributions due to diffraction and the GRIN potential, while the second line encodes the thermorefractive nonlocality through a nonlocal quartic term. With this choice, the Euler–Lagrange equation is Eq. (11), with the exception of losses, which we include phenomenologically a posteriori for consistency with the source of nonlocality. We assume the transverse field is a single Laguerre–Gauss (LG) mode with radial index p and azimuthal index m, and S=|m|,

(15)\mathrm{\Phi }(r,\theta ,z)=A(z)T_{pS}(r;\sigma (z))\exp \left[\mathrm{i}{\displaystyle \frac{b(z)}{2}}{r}^2\right]{\mathrm{e}}^{\mathrm{i}m\theta },

where \sigma (z) is a variational transverse width and T_{pS} is normalized as

(16){\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}rr|T_{pS}(r;\sigma ){|}^2=1.}

This implies that the power of the beam is simply expressed as P=2\pi |A{|}^2. The parameter b(z) regulates the phase profile over the radial direction. It is motivated by the known phase profile of paraxial beams, and was used in past literature for variational study of beam propagation [29]. The LG transverse profile is written as

(17)T_{pS}(r,\sigma )={\displaystyle \frac{1}{\sigma }}{C}_{pS}{\left({\displaystyle \frac{r}{\sigma }}\right)}^S{L}_p^S\left({\displaystyle \frac{r^2}{{\sigma }^2}}\right)\exp \left[-{\displaystyle \frac{r^2}{2{\sigma }^2}}\right],

with

(18)C_{pS}=\sqrt{{\displaystyle \frac{2p!}{(p+S)!}}},

where p is the radial number.

3.1 Effective Lagrangian: local terms

Inserting the ansatz into the full Lagrangian density and integrating over the transverse plane yields

(19)\begin{array}{rl}{L}_{\mathrm{l}\mathrm{o}\mathrm{c}}=& \mathrm{i}{A}^{\ast }\dot{A}\\ & -\frac{{\xi }_{pS}}{2}[{\sigma }^{-2}+({\mathrm{\Omega }}^2+b^2+\dot{b}){\sigma }^2]|A{|}^2,\end{array}

where the coefficients are

(20){\xi }_{pS}=S+2p+1,

as in Ref. [36]. The dot notation indicates differentiation with respect to the propagation variable z. These terms are identical to the local NPSE case [36, 41], except for the presence of the chirp term.

3.2 Effective action: nonlocal term

The thermal nonlocality enters the propagation Lagrangian in 1D through the last term of the Lagrangian density (14),

(21){\mathcal{L}}_{\text{non-loc}}=-{\displaystyle \frac{\gamma }{2}}|\mathrm{\Phi }({\mathit{r}}_{\perp },z){|}^2\mathrm{\Theta }[|\mathrm{\Phi }{|}^2]({\mathit{r}}_{\perp },z),

where the nonlocal thermal potential is defined as the convolution (12). The corresponding contribution to the effective one-dimensional Lagrangian L(z) is obtained by integrating over the transverse plane:

(22)\begin{array}{rl}{L}_{\text{non-loc}}=& \int {\mathrm{d}}^2{\mathit{r}}_{\perp }{\mathcal{L}}_{\text{non-loc}}\\ =& -\frac{\gamma }{4\pi }\int {\mathrm{d}}^2{\mathit{r}}_{\perp }|\mathrm{\Phi }({\mathit{r}}_{\perp },z){|}^2\\ & \times \int {\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }{K}_0\left(\frac{|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}\right)|\mathrm{\Phi }({\mathit{r}}_{\perp }^{\prime },z){|}^2,\end{array}

we remark that this expression is symmetric under the exchange of {\mathit{r}}_{\perp } and {\mathit{r}}_{\perp }^{\prime }. This is the general structure of the thermal nonlocal interaction at the Lagrangian level; using the LG ansatz (15) the nonlocal contribution becomes purely quartic in the axial amplitude A:

(23)L_{\text{non-loc}}=-{\displaystyle \frac{\gamma }{2}}|A{|}^4{\mathcal{I}}_{pS}(\sigma ),

where the transverse nonlocal integral is

(24)\begin{array}{rl}{\mathcal{I}}_{pS}(\sigma )=& \frac{1}{2\pi }\int {\mathrm{d}}^2{\mathit{r}}_{\perp }{\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }{K}_0\left(\frac{|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}\right)\\ & \times |T_{pS}(r;\sigma ){|}^2|T_{pS}(r^{\prime };\sigma ){|}^2.\end{array}

The dependence on the mode indices (p,S) and on the width \sigma enters only through {\mathcal{I}}_{pS}(\sigma ), which is the quantity to be evaluated explicitly.

We consider the nonlocal integral

(25)\begin{array}{rl}{\mathcal{I}}_{pS}(\sigma )=& \frac{1}{2\pi }\int {\mathrm{d}}^2{\mathit{r}}_{\perp }{\mathrm{d}}^2{\mathit{r}}_{\perp }^{\prime }{K}_0\left(\frac{|{\mathit{r}}_{\perp }-{\mathit{r}}_{\perp }^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}\right)\\ & \times {\rho }_{pS}(r;\sigma ){\rho }_{pS}(r^{\prime };\sigma ),\end{array}

where {\rho }_{pS}(r;\sigma )=|T_{pS}(r;\sigma ){|}^2 is radially symmetric. The integral can be simplified by using the Fourier convolution theorem in two dimensions (see Appendix A)

(26){\mathcal{I}}_{pS}(\sigma )=2\pi {\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}qq{\displaystyle \frac{L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}{1+L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2{q}^2}}{\tilde{\rho }}_{pS}(q;\sigma {)}^2,}

where

(27){\tilde{\rho }}_{pS}(q;\sigma )={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}rr{\rho }_{pS}(r;\sigma )J_0(qr),}

is the Hankel transform of \rho.

The effective Lagrangian in this case is obtained by combining the local part Eq. (19) with the nonlocal part Eq. (23), with the simplified value of the nonlocal contribution expressed implicitly by Eq. (26). The Euler–Lagrange equation for the axial amplitude follows from \mathrm{\partial }L/\mathrm{\partial }{A}^{\ast }=0, since L does not depend on {\dot{A}}^{\ast }. To have a consistent account of the thermal effects, we add the absorption in the equation of motion, taking into account the transverse shape of the mode

(28)\begin{array}{rl}\mathrm{i}\dot{A}=& -\mathrm{i}\frac{{\alpha }_0}{2}A\\ & +\frac{{\xi }_{pS}}{2}[{\sigma }^{-2}+({\mathrm{\Omega }}^2+b^2+\dot{b}){\sigma }^2]A\\ & +\gamma {\mathcal{I}}_{pS}(\sigma )|A{|}^{2}A.\end{array}

The first two terms correspond to the local contributions, while the third term is the nonlocal thermorefractive correction. The transverse width \sigma (z) enters the Lagrangian algebraically, so its equation of motion reduces to the condition \mathrm{\partial }L/\mathrm{\partial }\sigma =0. The local part yields

(29){\displaystyle \frac{\mathrm{\partial }{L}_{\mathrm{l}\mathrm{o}\mathrm{c}}}{\mathrm{\partial }\sigma }}={\xi }_{pS}|A{|}^2\left({\sigma }^{-3}-({\mathrm{\Omega }}^2+b^2+\dot{b})\sigma \right),

while the nonlocal part gives

(30){\displaystyle \frac{\mathrm{\partial }{L}_{\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{o}\mathrm{c}}}{\mathrm{\partial }\sigma }}=-{\displaystyle \frac{\gamma }{2}}|A{|}^4{\displaystyle \frac{\mathrm{d}{\mathcal{I}}_{pS}}{\mathrm{d}\sigma }}.

Collecting all terms, the condition \mathrm{\partial }L/\mathrm{\partial }\sigma =0 yields the algebraic equation

(31)\begin{array}{rl}& {\xi }_{pS}|A{|}^2\left[{\sigma }^{-3}-\left({\mathrm{\Omega }}^2+b^2+\dot{b}\right)\sigma \right]\\ & -\frac{\gamma }{2}|A{|}^4\frac{\mathrm{d}{\mathcal{I}}_{pS}}{\mathrm{d}\sigma }=0.\end{array}

Finally the Euler−Lagrange equation for the radial phase chirp parameter is

(32){\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}({\sigma }^2|A{|}^2)=2b{\sigma }^2|A{|}^2.

Equations (28), (31) and (32) define a dimensional reduction: the axial field evolves under a \sigma-dependent nonlinearity, while the transverse width is determined by the balance between the local and nonlocal thermal contributions. These equations constitute the main result of the work. While they are not immediately written with an algebraic nonlinearity, similarly to the NPSE approach [34], closed forms are available for specific choices of the transverse distributions, like the Gaussian one. We remark that, unlike the previously studied systems, where it was possible to express explicitly the value of \sigma and substitute it back into the equation for the axial amplitude, in this case, since {\mathcal{I}}_{pS} is a complicated function of \sigma, this is no longer possible. However, it is remarkable that closed-form expressions can be found for {\mathcal{I}}_{0S}, i.e., the case with p=0, as shown in Appendix B.

3.3 Special case: Gaussian beam

For the fundamental Gaussian mode (p=0, m=0, hence S=0), the transverse profile reduces to

(33)T_{00}(r;\sigma )={\displaystyle \frac{\sqrt{2}}{\sigma }}\exp \left(-{\displaystyle \frac{r^2}{2{\sigma }^2}}\right),

so that the radial intensity is

(34){\rho }_{00}(r;\sigma )=|T_{00}(r;\sigma ){|}^2={\displaystyle \frac{2}{{\sigma }^2}}\exp \left(-{\displaystyle \frac{r^2}{{\sigma }^2}}\right).

Its Hankel transform is again a Gaussian,

(35){\tilde{\rho }}_{00}(q)={\mathrm{e}}^{-{\sigma }^2{q}^2/4}.

Inserting this into Eq. (26) yields

(36){\mathcal{I}}_{00}(\sigma )=2\pi L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2{\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}q{\displaystyle \frac{q{\mathrm{e}}^{-{\sigma }^2{q}^2/2}}{1+L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2{q}^2}}.}

This integral can be performed in terms of the exponential integral \mathrm{Ei} [40]

(37){\mathcal{I}}_{00}(\sigma )=-\pi \exp \left({\displaystyle \frac{{\sigma }^2}{2L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}}\right)\mathrm{Ei}\left(-{\displaystyle \frac{{\sigma }^2}{2L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}}\right),

its derivative reads

(38){\displaystyle \frac{\mathrm{d}{\mathcal{I}}_{00}}{\mathrm{d}\sigma }}={\displaystyle \frac{\sigma }{L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}}{\mathcal{I}}_{00}(\sigma )-{\displaystyle \frac{2\pi }{\sigma }}.

The corresponding integrals for S\ge 1 are less straightforward and can be obtained as discussed in Appendix B.

We remark that, if the thermal nonlocal term (37) is considered in the limit of L_{\mathrm{e}\mathrm{f}\mathrm{f}}\to 0, the equations for the beam evolution reduce to the ones derived in the variational model of Ref. [29] (the limit is also discussed in the derivation of [21, Eq. (7)]). Unlike the treatment in Ref. [29], in the present model it is not possible to simply generalize the results to pulses, as the response time of the thermal properties would make the analysis much more involved, even neglecting chromatic dispersion. Finally, in the special case of a linear evolution, Eq. (31) becomes the well-known equation for beam width evolution (see also Ref. [42]).

4 Numerical simulations

To validate the reduced model, we performed direct simulations of the full two-dimensional paraxial diffraction equation with thermal nonlocality using a standard split-step Fourier method. The linear propagation was computed in the transverse momentum domain, while the nonlinear thermorefractive potential was updated at each axial step through convolution with the kernel discussed in Eq. (10). This approach follows established numerical strategies for nonlocal nonlinear Schrödinger equations [5, 24], and allows for accurate benchmarking of beam evolution, including dynamical regimes where the dimensionally reduced model is expected to perform poorly, like the ones of dispersive shock waves in the self-defocusing regime [21], and symmetry-breaking azimuthal instability [24]. The reduced variational model represented by Eqs. (28, 31, 32) is solved using a standard fourth-order Runge−Kutta method.

4.1 Stationary beam profiles

Starting from the equations of motion, we consider the conservative case {\alpha }_0=0. In this regime Eq. (28) imposes that the beam power P=2\pi |A{|}^2 is conserved. Using this information and substituting in Eq. (32), we obtain that the derivative \dot{\sigma } is proportional to b: to have a stationary value, we require b=0. A particularly interesting analysis is then to solve Eq. (31) for the stationary value of \sigma. In the following, we set |A{|}^2=1. The aforementioned Eq. (31) in the static case can be interpreted as the condition of stationarity of the following effective potential

(39)V_{\mathrm{e}\mathrm{f}\mathrm{f}}(\sigma )={\displaystyle \frac{{\xi }_{pS}}{2}}\left({\displaystyle \frac{1}{{\sigma }^2}}+{\mathrm{\Omega }}^2{\sigma }^2\right)+{\displaystyle \frac{\gamma |A{|}^2}{2}}{\mathcal{I}}_{0S}(\sigma ).

The structure of the effective potential V_{\mathrm{e}\mathrm{f}\mathrm{f}}(\sigma ) changes significantly with the degree of thermal nonlocality, as shown in Fig. 1. In the defocusing regime [Fig. 1(a)], increasing L_{\mathrm{e}\mathrm{f}\mathrm{f}} shifts the minimum of the potential toward larger widths, indicating weaker transverse confinement and a broader stationary beam. In the focusing regime [Fig. 1(b)], the nonlocal contribution changes the potential more dramatically: the attractive term induced by diffusion becomes important at small \sigma, and the depth of the potential well increases with L_{\mathrm{e}\mathrm{f}\mathrm{f}}.

To assess the validity of the dimensional reduction, we first compare the stationary states obtained via imaginary time propagation [43–45] with the predictions of the variational method. Figure 2 compares the radial intensity profiles |\mathrm{\Phi }(\mathit{r}){|}^2 in a defocusing regime, \gamma =0.5, with a graded-index potential, \mathrm{\Omega }=1. In the case of short-range nonlocality of panels (a) and (b), having L_{\mathrm{e}\mathrm{f}\mathrm{f}}=0.5, the variational prediction shows excellent agreement with the exact numerical solutions for both vortex orders S=1 and S=2. However, as the thermal response becomes more nonlocal [panels (c) and (d), with L_{\mathrm{e}\mathrm{f}\mathrm{f}}=5.0], the physical beam profile deforms and departs from the functional form assumed in the variational ansatz. A similar analysis is done in the focusing case and without confining potential in Fig. 3, with a parameter \gamma =-0.22. In this analysis the nonlinear parameter and the effective length L_{\mathrm{e}\mathrm{f}\mathrm{f}} were chosen to prevent the solver from collapsing to a filamented state. Among the simulations with L_{\mathrm{e}\mathrm{f}\mathrm{f}}=1.2 [panels (a) and (b)], the one with high vorticity in panel (b) displays a poor matching of the variational results. This is interpreted by considering the fact that this choice of parameters leads to a situation similar to the one represented in panel (b) of Fig. 1: when diminishing the effective length, the curvature of the effective potential diminishes more and more. This makes the system much more sensitive to small variations in the radial shape and deviations from the Laguerre–Gauss shape. The variational model still reproduces the qualitative broadening of the beam in exact solutions in the strongly nonlocal regime L_{\mathrm{e}\mathrm{f}\mathrm{f}}=5.0 [panels (c) and (d)].

4.2 Beam propagation

Figure 4 shows the dynamical evolution of a vortex beam S=1 propagating with losses set to \alpha =0.2 with strong thermal nonlocality L_{\mathrm{e}\mathrm{f}\mathrm{f}}=40. The axial amplitude |A(z){|}^2 decays monotonically due to linear absorption [panel (a)], while the transverse width \sigma (z) undergoes damped breathing oscillations constrained by the graded-index potential [panel (b)]. The correspondence between the reduced one-dimensional variational model (solid red lines) and the full two-dimensional split-step Fourier simulations (blue markers) is substantial. This demonstrates that, despite the moderate discrepancies observed in the static mode profiles, the variational method accurately captures the overall beam evolution dictated by the interplay of diffraction, thermal diffusion, and the parabolic trapping potential, especially for small propagation distances. We remark that the inclusion in the ansatz of the radial chirp term as described in Eq. (15) is important. This can be deduced by analyzing Eq. (31), where, in the absence of the variable b, any non-dissipating beam would have a fixed width for the entire propagation distance.

5 Conclusions

In this work we have examined a reduced one-dimensional description of beam propagation in graded-index media with thermal nonlocal nonlinearities, derived from the full spatial nonlinear Schrödinger equation coupled to a screened-Poisson model. By extending the variational methodology to include a nonlocal convolution term, we obtained an effective action in which the thermal interaction appears through a mode- and width-dependent nonlocal energy functional. The use of an annular Laguerre–Gauss trial function enabled closed form expressions for the Hankel transform of the intensity and yielded a compact formulation of the reduced dynamics. The variational reduction is performed for the conservative system, and linear absorption is then reintroduced phenomenologically. The resulting evolution equations show that thermal nonlocality alters the dependence of the nonlinear interaction on the transverse width and mode order. The analysis also indicates that a single-width ansatz is insufficient to capture the dynamics in nonlocal media, and that a radial phase-chirp parameter is needed.

The approach may be extended to multimode or partially coherent fields, where thermal diffusion mediates effective intermodal coupling. Further generalizations may include saturable responses or time-dependent thermal relaxation, to broaden the applicability of the method to pulsed systems as well, where the interplay with local electronic Kerr effect and chromatic dispersion would be interesting.

6 Appendix A: Fourier−Bessel convolution

To reduce this expression we use the Hankel transform of the kernel,

(40)\tilde{K}(q)={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}rr{K}_0\left({\displaystyle \frac{r}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\right)J_0(qr),}

which provides the inversion formula

(41)K_0\left({\displaystyle \frac{|\mathit{r}-{\mathit{r}}^{\prime }|}{L_{\mathrm{e}\mathrm{f}\mathrm{f}}}}\right)={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}qq\tilde{K}(q)J_0(q|\mathit{r}-{\mathit{r}}^{\prime }|).}

Substituting this into (25) and exchanging the order of integration yields (we drop the subscripts of \rho for brevity)

(42){\mathcal{I}}_{pS}(\sigma )={\displaystyle \frac{1}{2\pi }}{\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}qq\tilde{K}(q)J(q),}

(43)J(q)={\displaystyle \int {\mathrm{d}}^2\mathit{r}{\mathrm{d}}^2{\mathit{r}}^{\prime }\rho (r)\rho (r^{\prime })J_0(q|\mathit{r}-{\mathit{r}}^{\prime }|).}

Using the addition theorem for Bessel functions together with azimuthal integration [9, §8.530],

(44){\displaystyle {\int }_0^{2\pi }\mathrm{d}\theta {\displaystyle {\int }_0^{2\pi }\mathrm{d}{\theta }^{\prime }{J}_0(q|\mathit{r}-{\mathit{r}}^{\prime }|)=(2\pi {)}^2{J}_0(qr)J_0(q{r}^{\prime }),}}

we obtain

(45)\begin{array}{rl}J(q)=& (2\pi {)}^2\left[{\int }_0^{\mathrm{\infty }}\mathrm{d}rr\rho (r)J_0(qr)\right]\\ & \times \left[{\int }_0^{\mathrm{\infty }}\mathrm{d}{r}^{\prime }{r}^{\prime }\rho (r^{\prime })J_0(q{r}^{\prime })\right].\end{array}

Using the definition of the Hankel transform of the radial intensity, in Eq. (27), the nonlocal integral becomes

(46){\mathcal{I}}_{pS}(\sigma )=2\pi {\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}qq\tilde{K}(q)|{\tilde{\rho }}_{pS}(q){|}^2.}

Finally, inserting the explicit kernel transform

(47)\tilde{K}(q)={\displaystyle \frac{L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}{1+L_{\mathrm{e}\mathrm{f}\mathrm{f}}^2{q}^2}},

we retrieve the compact expression (26).

7 Appendix B: Closed-form Hankel transform for Laguerre−Gauss intensities

We compute the Hankel transform {\tilde{\rho }}_{pS}(q) for the LG intensity. Using the explicit expressions defining (15), one has

(48){\rho }_{pS}(r)={\displaystyle \frac{C_{pS}^2}{{\sigma }^{2S+2}}}{r}^{2S}{\left[L_p^S\left({\displaystyle \frac{r^2}{{\sigma }^2}}\right)\right]}^2\exp \left(-{\displaystyle \frac{r^2}{{\sigma }^2}}\right).

Accordingly,

(49){\tilde{\rho }}_{pS}(q)={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}rr{\rho }_{pS}(r;\sigma )J_0(qr).}

A particularly compact expression is obtained for the vortical modes with p=0 and S\ne 0, for which

(50){\rho }_{0S}(r)=B_{0S}{r}^{2S}\exp \left(-{\displaystyle \frac{r^2}{{\sigma }^2}}\right),

with

(51)B_{0S}={\displaystyle \frac{2}{{\sigma }^{2(S+1)}S!}}.

In this case the Hankel transform can be expressed using confluent hypergeometric functions

(52){\tilde{\rho }}_{0S}(q)={}_1{F}_1\left(1+S;1;-{\displaystyle \frac{q^2{\sigma }^2}{4}}\right),

where {}_1{F}_1 denotes the confluent hypergeometric function of the first kind. Equation (52) provides a compact analytic expression for the Hankel transform of the annular (p=0) LG vortices, and is particularly useful when evaluating the nonlocal energy {\mathcal{I}}_{0S} through Eq. (26). Indeed, performing the integration for a fixed value of S one finds that the quantity {\mathcal{I}}_{0S} admits an explicit analytical expression in terms of polynomials in the variable \sigma /L_{\mathrm{e}\mathrm{f}\mathrm{f}} and expressions containing the exponential integral function \mathrm{E}\mathrm{i}. In the case S=0 this procedure reproduces the closed form reported in Eq. (37). For S>0, the derivative d{\mathcal{I}}_{0S}/d\sigma required in Eq. (31) is likewise available analytically once the value of {\mathcal{I}}_{0S} is determined.

The scaling of the nonlocal coefficient is analyzed in Fig. A1. The S=0 case (solid blue line) defines an upper bound, as increasing the azimuthal index S systematically decreases {\mathcal{I}}_{0S}. The dependence on \sigma contrasts with the variational analysis obtained with local nonlinearity, where interaction energy scales as 1/{\sigma }^2 independently of S. In the nonlocal case, there is no simple algebraic decay valid for all values of \sigma.

References

[1] C. Conti , M. Peccianti , and G. Assanto , Observation of optical spatial solitons in a highly nonlocal medium, Phys. Rev. Lett. 92(11), 113902 (2004)

[2] O. Bang , W. Krolikowski , J. Wyller , and J. J. Rasmussen , Collapse arrest and soliton stabilization in nonlocal nonlinear media, Phys. Rev. E 66(4), 046619 (2002)

[3] W. Królikowski , O. Bang , N. I. Nikolov , D. Neshev , J. Wyller , J. J. Rasmussen , and D. Edmundson , Modulational instability, solitons and beam propagation in spatially nonlocal nonlinear media, J. Opt. B 6(5), S288 (2004)

[4] S. K. Turitsyn , Spatial dispersion of nonlinearity and stability of multidimensional solitons, Theor. Math. Phys. 64(2), 797 (1985)

[5] W. Krolikowski , O. Bang , J. J. Rasmussen , and J. Wyller , Modulational instability in nonlocal nonlinear Kerr media, Phys. Rev. E 64(1), 016612 (2001)

[6] D. Suter and T. Blasberg , Stabilization of transverse solitary waves by a nonlocal response of the nonlinear medium, Phys. Rev. A 48(6), 4583 (1993)

[7] C. Conti , G. Ruocco , and S. Trillo , Optical spatial solitons in soft matter, Phys. Rev. Lett. 95(18), 183902 (2005)

[8] Z. Xu , Y. V. Kartashov , and L. Torner , Soliton mobility in nonlocal optical lattices, Phys. Rev. Lett. 95(11), 113901 (2005)

[9] G. Gligorić , A. Maluckov , L. Hadžievski , and B. A. Malomed , Collapse instability of solitons in the nonpolynomial Schrödinger equation with dipole–dipole interactions, J. Phys. At. Mol. Opt. Phys. 42(14), 145302 (2009)

[10] A. Gammal , L. Tomio , and T. Frederico , Critical numbers of attractive Bose-Einstein condensed atoms in asymmetric traps, Phys. Rev. A 66(4), 043619 (2002)

[11] L. Tanzi , E. Lucioni , F. Fama’ , J. Catani , A. Fioretti , C. Gabbanini , R. N. Bisset , L. Santos , and G. Modugno , Observation of a dipolar quantum gas with metastable supersolid properties, Phys. Rev. Lett. 122(13), 130405 (2019)

[12] A. G. Litvak , Self-focusing of powerful light beams by thermal effects, ZhETF Pisma Redaktsiiu 4, 341 (1966)

[13] P. Brochard , V. Grolier-Mazza , and R. Cabanel , Thermal nonlinear refraction in dye solutions: A study of the transient regime, J. Opt. Soc. Am. B 14(2), 405 (1997)

[14] S. J. Sheldon , L. V. Knight , and J. M. Thorne , Laserinduced thermal lens effect: A new theoretical model, Appl. Opt. 21(9), 1663 (1982)

[15] C. Rotschild , O. Cohen , O. Manela , M. Segev , and T. Carmon , Solitons in nonlinear media with an infinite range of nonlocality: First observation of coherent elliptic solitons and of vortex-ring solitons, Phys. Rev. Lett. 95(21), 213904 (2005)

[16] A. Minovich , D. N. Neshev , A. Dreischuh , W. Krolikowski , and Y. S. Kivshar , Experimental reconstruction of nonlocal response of thermal nonlinear optical media, Opt. Lett. 32(12), 1599 (2007)

[17] B. Alfassi , C. Rotschild , O. Manela , M. Segev , and D. N. Christodoulides , Boundary force effects exerted on solitons in highly nonlocal nonlinear media, Opt. Lett. 32(2), 154 (2007)

[18] A. Alberucci and G. Assanto , Propagation of optical spatial solitons in finite-size media: Interplay between nonlocality and boundary conditions, J. Opt. Soc. Am. B 24(9), 2314 (2007)

[19] A. Alberucci , C. P. Jisha , N. F. Smyth , and G. Assanto , Spatial optical solitons in highly nonlocal media, Phys. Rev. A 91(1), 013841 (2015)

[20] A. Alberucci , M. Peccianti , and G. Assanto , Nonlinear bouncing of nonlocal spatial solitons at the boundaries, Opt. Lett. 32(19), 2795 (2007)

[21] N. Ghofraniha , C. Conti , G. Ruocco , and S. Trillo , Shocks in nonlocal media, Phys. Rev. Lett. 99(4), 043903 (2007)

[22] G. Marcucci , D. Pierangeli , S. Gentilini , N. Ghofraniha , Z. Chen , and C. Conti , Optical spatial shock waves in nonlocal nonlinear media, Adv. Phys. X 4(1), 1662733 (2019)

[23] M. S. Bigelow , P. Zerom , and R. W. Boyd , Breakup of ring beams carrying orbital angular momentum in sodium vapor, Phys. Rev. Lett. 92(8), 083902 (2004)

[24] A. I. Yakimenko , Y. A. Zaliznyak , and Y. Kivshar , Stable vortex solitons in nonlocal self-focusing nonlinear media, Phys. Rev. E 71(6), 065603 (2005)

[25] D. Briedis , D. E. Petersen , D. Edmundson , W. Krolikowski , and O. Bang , Ring vortex solitons in nonlocal nonlinear media, Opt. Express 13(2), 435 (2005)

[26] V. Biloshytskyi , A. Oliinyk , P. Kruglenko , A. Desyatnikov , and A. Yakimenko , Vortex nucleation in nonlocal nonlinear media, Phys. Rev. A 99(4), 043835 (2019)

[27] F. Poletti and P. Horak , Description of ultrashort pulse propagation in multimode optical fibers, J. Opt. Soc. Am. B 25(10), 1645 (2008)

[28] A. Antikainen , L. Rishøj , B. Tai , S. Ramachandran , and G. P. Agrawal , Fate of a soliton in a high order spatial mode of a multimode fiber, Phys. Rev. Lett. 122(2), 023901 (2019)

[29] M. Karlsson , D. Anderson , and M. Desaix , Dynamics of self-focusing and self-phase modulation in a parabolic index optical fiber, Opt. Lett. 17(1), 22 (1992)

[30] P. Parra-Rivas , Y. Sun , and S. Wabnitz , Dynamics of three-dimensional spatiotemporal solitons in multimode waveguides, Opt. Commun. 546, 129749 (2023)

[31] O. V. Shtyrina , M. P. Fedoruk , Y. S. Kivshar , and S. K. Turitsyn , Coexistence of collapse and stable spatiotemporal solitons in multimode fibers, Phys. Rev. A 97(1), 013841 (2018)

[32] A. A. Aguilar-Cardoso , C. Li , T. J. B. Luck , M. F. Ferrer-Garcia , J. Upham , J. S. Lundeen , R. W. Boyd , M. F. FerrerGarcia , J. Upham , J. S. Lundeen , and R. W. Boyd , Tailoring spatial modes produced by stimulated parametric down-conversion, Phys. Rev. A 112(4), 043541 (2025)

[33] L. Sá and J. Vieira , Self-focusing of multiple interacting Laguerre−Gauss beams in Kerr media, Phys. Rev. A 100(1), 013836 (2019)

[34] L. Salasnich , A. Parola , and L. Reatto , Condensate bright solitons under transverse confinement, Phys. Rev. A 66(4), 043603 (2002)

[35] L. Salasnich , A. Parola , and L. Reatto , Effective wave equations for the dynamics of cigar-shaped and diskshaped Bose condensates, Phys. Rev. A 65(4), 043614 (2002)

[36] F. Lorenzi and L. Salasnich , Variational approach to multimode nonlinear optical fibers, Nanophotonics 14(6), 805 (2025)

[37] E. Karimi , R. W. Boyd , P. de la Hoz , H. de Guise , J. Rehacek , Z. Hradil , A. Aiello , G. Leuchs , and L. L. Sánchez-Soto , Radial quantum number of Laguerre−Gauss modes, Phys. Rev. A 89(6), 063813 (2014)

[38] B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, 2 Volume Set, John Wiley & Sons, 2019

[39] A. Alberucci , A. Piccardi , M. Peccianti , M. Kaczmarek , and G. Assanto , Propagation of spatial optical solitons in a dielectric with adjustable nonlinearity, Phys. Rev. A 82(2), 023806 (2010)

[40] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Vol. 55, Courier Corporation, 1965

[41] F. Lorenzi and L. Salasnich , Atomic soliton transmission and induced collapse in scattering from a narrow barrier, Sci. Rep. 14(1), 4665 (2024)

[42] H. R. Lewis , Classical and quantum systems with time-dependent harmonic-oscillator-type Hamiltonians, Phys. Rev. Lett. 18(13), 510 (1967)

[43] B. Jackson, Vortices in Trapped Bose Einstein Condensates, Doctoral Thesis, Durham University, 2000

[44] I. Tikhonenkov , B. A. Malomed , and A. Vardi , Vortex solitons in dipolar Bose−Einstein condensates, Phys. Rev. A 78(4), 043614 (2008)

[45] W. Bao and Q. Du , Computing the ground state solution of Bose–Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput. 25(5), 1674 (2004)

[46] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, 2014

RIGHTS & PERMISSIONS

Higher Education Press