Multidimensional localized states in externally driven Kerr cavities with a parabolic spatiotemporal potential: a dimensional connection

In this work, we study the bifurcation structures and the stability of multidimensional localized states within coherently driven Kerr optical cavities with parabolic potentials in 1D, 2D, and 3D systems. Based on symmetric considerations, we transform higher-dimensional models into a single 1D model with a dimension parameter. This transformation not only yields a substantial reduction in computational complexity, but also enables an efficient examination of how dimensionality impacts the system dynamics. In the absence of nonlinearity, we analyze the eigenstates of the linear systems. This allows us to uncover a heightened concentration of the eigenmodes at the center of the potential well, while witnessing a consistent equal spacing among their eigenvalues, as the dimension parameter increases. In the presence of nonlinearity, our findings distinctly reveal that the stability of the localized states diminishes with increasing dimensionality. This study offers an approach to tackling high-dimensional problems, shedding light on the fundamental dimensional connections among radially symmetric states across different dimensions, and providing valuable tools for analysis.


I. INTRODUCTION
In conservative physical systems, solitons or solitary waves are self-sustained localized wave packets, that maintain their shape unchanged while propagating through nonlinear dispersive or diffractive media.These nonlinear waves have been explored in a plethora of different physical contexts, such as Bose-Einstein condensates, plasmas, hadron matter, gravitation, and optics.[1][2][3][4].A completely different class of localized waves appears in open or dissipative systems, which experience a continuous energy exchange with the surrounding media.These localized dissipative states, hereafter LS, are also known as dissipative solitons [5].In dissipative optical systems, LS are under the balance between gain and loss, and manifest in systems of different dimensionality: e.g., 1D temporal LS involve a balance between nonlinearity and dispersion, while 2D LS form in the presence of nonlinearity and diffraction.The theoretical prediction of the emergence of 2D spatial LS [6] was experimentally confirmed using a large-area vertical cavity surface emitting laser with a coherent holding field [7].In parallel, the emergence of temporal cavity solitons, the counterparts of spatial cavity solitons, was predicted by substituting diffraction with dispersion [8,9].Subsequent experimental confirmations of temporal dissipative Kerr solitons in fiber optics [10] and microresonators [11] highlighted their connection with the generation of coherent optical frequency combs.These combs unlock a new era of technological possibilities, including optical buffering, optical clocks, dualcomb spectroscopy, and frequency synthesizers.[10,12].
Nevertheless, theoretical studies of 3D dissipative solitons [21][22][23][24][25][26][27][28][29] still face many fundamental and practical challenges.From a fundamental perspective, an important problem is the fragility of multidimensional states against any perturbation [3].A typical example of such fragility is the spatial or spatiotemporal wave collapse, that is commonly experienced by multidimensional solitons in optical materials with Kerr nonlinearity [3,21,30,31].To arrest, at least partially, these catastrophic instabilities, various strategies have been proposed, including dynamical regularization of collapse [32], the use of saturable, nonlocal, and competing nonlinearities [33][34][35][36][37][38], rapid longitudinal variations of the material parameters [39], the use of both static and twisted optical lattices [40][41][42][43][44] and other wave confinement methods [45][46][47].From the modeling perspective, a key challenge arises when considering the high-demanding computational load for the study of multidimensional LS.When the dimensionality increases, the computational complexity grows significantly larger, posing a bottleneck for numerical investigations of high-dimensional solitons.Typically, the workload grows quadratically for 2D problems and cubically for 3D problems, yielding an insufficient number of points in the discrete grids.For instance, when employing a thousand discretization points in 1D corresponds to a million points in 2D, and a billion samples in 3D [see Fig. 1(a)].Thereby, in high-dimensional problems, the demand for an excessively large number of grid points poses a key obstacle, since for practical reasons one is forced to use very low-resolution meshes, which may result in significant numerical errors.Moreover, the essential question that we try to address here is: how and why dimensionality affects the stability of LS?
In this work, we answer this question by performing an extensive bifurcation analysis of radially symmetric LS, in the context of externally driven Kerr cavities with a spatiotemporal parabolic confining potential [29].When considering this symmetry, one may reduce the description of multidimensional LS [see Fig. an arbitrary number of dimensions.By doing so, we solve the issue of the heavy numerical computation load for simulating the dynamics of highly dimensional states.Moreover, we are able to unveil why the latter are much more fragile than their 1D counterparts.Clearly, the dimensional reduction based on the radial symmetry inherently removes the possibility of studying other types of more complex states, which may emerge from the radial ones [48][49][50].On the one hand, this has the benefit of considerably facilitating the analysis; however, it leaves out the possible emergence of radially asymmetric states in the original systems.In any case, our procedure has the merit that it enables us to unveil the fundamental dimensional connection between radially symmetric states with different dimensions.This article is organized as follows.In Section II we introduce the mathematical nonlinear model describing the evolution of the optical field in the cavity.Section III explores the properties of the eigenmodes, which are associated with the linear problem.Later, in Section IV we study how, in the nonlinear regime, such modes lead to nonlinear LS.The bifurcation and stability analysis of these states is presented in Section V, where the dimensionality implications are carefully addressed.Finally, in Section VI we present a short discussion and draw our main conclusions.

II. THE RADIAL MODEL
In the mean-field approximation, the spatiotemporal field dynamics in externally driven Kerr cavities, in the presence of diffraction and/or anomalous dispersion, can be described in terms of the driven dissipative Gross-Pitaevskii equation where A(r,t) is the electric field envelope, r ∈ R D with D being the dimension, t is time, δ is the cavity phase detuning, P is the intensity of the external driving field or pump, and the losses are normalized to −1. ∇ 2 is the D-dimensional Laplacian and the term r 2 represents the parabolic potential Note that special cases of Eq. ( 1) describe the dynamics of temporal states A(r,t) = A(τ,t) in dispersive cavities when D = 1, of spatial structures in diffractive cavities when D = 2 [i.e., A(r,t) = A(x, y,t)], and for D = 3 of spatiotemporal states A(r,t) = A(τ, x, y,t), with diffraction and dispersion acting simultaneously.By polar or spherical coordinate transformations, radially symmetric solutions A(r,t) = A(r,t) of Eq. ( 1) in the region (0, R) obey (2) together with the boundary condition ∂ r A(r)| r=0, or, r=R = 0. Here, the term r −1 (D − 1)∂ r A, arising from the Ddimensional Laplacian ∇ 2 A, is the pivotal element for distinguishing the dimensionality of the problem.For D > 1, this term introduces an effective position-dependent field flow: in its presence, the field experiences an elevated flow rate when r → 0, which leads to an effective centripetal force, weighted by D, compelling the field to concentrate towards the center r = 0.
Equation ( 1) has been previously used to study 1D temporal LS in dispersive Kerr cavities with external phase modulation [51][52][53], and for analyzing the stabilization of highorder spherical (3D) light bullets, when diffraction and dispersion act simultaneously in the presence of a spatiotemporal potential [29].Similarly, a diffractive cavity with a spatial parabolic inhomogeneity can be described by Eq. ( 2) when D = 2. Building on these results, now we consider D as a continuous parameter.This permits us to explore the homotopic connection between radially symmetric LS across different dimensions, as well as the impact of the dimensionality on the stability of LS.
Physically, the radial localization around r = 0 is prompted by two main reasons.First, the presence of a parabolic potential breaks the translational invariance of LS, effectively forcing the field to be distributed around the center of the potential well.Second, the presence of a coherent and spatially uniform pump P effectively suppresses asymmetric field modes.This is because the energy pumping rate for each mode is intimately linked with the transverse integral of the modal profiles [52], thereby favoring the dynamics dominated by symmetric solutions.
All of the results presented here have been mutually validated through path-continuation (homotopic) techniques by using pde2path [54], as well as by means of full dimensional direct numerical simulations (DNS).Similarly, the LS stability was evaluated through both spectral analysis and DNS.

III. LINEAR EIGENMODES
Building upon the foundation laid by Eq. ( 2), let us start by analyzing the linear eigenvalue problem where ĤD is the linear operator and ψ D,n and δ D,n are the linear eigenmodes and eigenvalues, parametrized by the dimension parameter D and the mode number n.
In order to visually represent these radial modes, we mirror their profiles to negative r-values, so as to construct full 1D profiles.These eigenmode profiles in Fig. 2(a-c) are vertically shifted according to their respective eigenvalues for the 1D, 2D, and 3D scenarios, respectively.Even though the dimension of a physical system is discrete, we analyze the modification of the eigenvalues δ D,n by continuously varying the dimension parameter D [see Fig. 2(g)].
As it can be seen, the eigenvalues remain equidistant for different D. In other words, increasing the dimension does not affect the spacing between eigenvalues; instead, it uniformly shifts upwards all the eigenvalues in the potential as D increases.The mode separation g depends on the coefficients of diffraction (C 1 ) and dispersion (C 2 ) through the expression g = 2 4C 1 /C 2 .Note that, in our normalized Eq.( 2), C 1 = C 2 = 1.Thus, the relation these eigenvalues reads The mode profiles in Fig. 2(a-c) exhibit distinct characteristics: they manifest as symmetric Hermite-Gaussian (HG) modes 1D, and as symmetric Laguerre-Gaussian (LG) modes in 2D, and 3D spherical modes, respectively.
As the mode order n increases, the spatial profile of the eigenmodes grows wider within the potential.However, for higher dimensions, the modes become more tightly localized around the center of the potential well.For high-order modes, the highest intensity peaks are located at the potential boundaries in 1D, but they shift toward the center, and the central peak becomes dominant for 2D and 3D.This observation can be quantified through the calculation of the amplitude-weighted modal width which is shown in Fig. 2(h).The concentration of the modal profiles towards the center of the potential is a pivotal dimension-related signature.
The resonant modes for different eigenvalues are verified by numerically solving their associated equations within original number of dimensions (Eq.( 1)), in the absence of Kerr nonlinearity i|A| 2 A [see Figs.2(d, e, f (a, b, c), which confirms their appropriateness in describing the resonant modes of the linear cavity.For other pump values, the fields exhibit analogous distributions as they correspond to solutions in linear scenarios.

IV. BIFURCATION STRUCTURE OF NONLINEAR LOCALIZED STATES AND DIMENSIONAL CONNECTION
In order to elucidate the influence of Kerr nonlinearity on the system dynamics, let us compare the changes of linear and nonlinear LS when scanning the cavity detuning δ for a fixed pump amplitude P = 1.3.To do that, we may compare the variation of the norm N = |A(r)| 2 dr/L associated with the 1D radial profiles with detuning δ , where L = 30 is twice the radial domain size.These results are illustrated in the bifurcation diagrams of Figs.3(a,b,c) for D = 1, 2, 3.
In the absence of Kerr nonlinearity [see gray curves in Figs.3(a,b,c)], all solutions are stable, aligning linear resonance peaks precisely with their eigenvalues [Figs.2(a,b,c,g) for D = 1, 2, 3].The introduction of Kerr nonlinearity shifts the resonance peaks in Figs.3(a,b,c) towards positive δvalues, leading to both stable (depicted by black solid curves) and unstable solutions (depicted by red dashed curves).This shift is due to the fact that modes resonate with an effective nonlinear cavity detuning i(|A| 2 − δ D,n ) [see Eq. ( 2)], owing to the intensity-dependent phase shift i|A| 2 .In order to maintain the effective detuning nearly constant at each resonance, the linear detuning δ undergoes a substantial adjustment with increasing values of the intensity |A| 2 .In order words, higher intensities leads to larger values of δ , leading to the bending of resonance peaks and to the formation of fold bifurcations (FB).
Despite sharing an identical pump amplitude P = 1.3, there are noteworthy differences among the scenarios ob- served for different dimensions.In 1D [see Fig. 3(a)], a single pair of fold bifurcation points emerges at the first nance peak, and a slight tilting is observed.This leads to the coexistence of two different LS for the same δ -value (i.e., bistability).Two examples of such states are depicted in panels Figs.3(i,ii) for δ = 1.2 [see labels in Fig. 3(a)].
In 2D [see Fig. Notably, in 3D all solutions become unstable, except those on the lower branch of the first resonance peak; moreover, much stronger resonance tilts are observed in this case.Solution (v) [see Fig. 3(v)] becomes unstable, while solution (vi) [see Fig. 3(vi)] remains stable for δ = 1.2, as shown in Figs.3(c,v,vi).In this case, we also display the 3D reconstruction of the field A(x, y, τ) (see insets in Figs.3(v,vi)).To validate the results shown in Figs.3(i-vi), including their linear stability, we have performed DNS in the original 1D, 2D and 3D models.
Building upon this insight, we explored the modifications of the bifurcation structures across a continuous variation of the dimensionality of the system, by harnessing the dimension parameter D as a continuous variable, even though the dimension parameter is inherently a discrete quantity.To examine the validity of this approach, we chose a pair of LS solutions within the bistable region, for fixed values of P = 1.3 and δ = 1.2 [see the gray vertical lines in Figs.3(ac)], and performed their path-continuation in D. The outcome of these computations is shown in Fig. 3(d), where we plot the modification of the maximal norm of the radial profile |A(r)| as a function of D. In this diagram, the solutions (i-vi) at D = 1, 2, 3 correspond exactly to the solutions (i-vi) indicated with vertical gray lines in Figs.3(a-c) and Figs.3(ivi), respectively.In other words, one may effectively investigate the evolution of multidimensional solitons across their dimensionality, by simply exploiting a single dimension parameter D.
The spectral (i.e., linear) stability analysis carried out for this diagram, limited to the presence of radial perturbations, indicates that the LS states (ii,iv,vi) [see bottom branch in Fig. 3(d)] are stable.In contrast, in the upper branch, only the 1D and 2D LS are stable.After crossing D = 2, a Hopf bifurcation (HB) emerges, which destabilizes this branch.It is worth noting that despite obtaining solutions (i, iii, v) under identical parameters δ and P, their field amplitude peaks exhibit a substantial increase with dimensionality.In Fig. 3(d), we depict the maximum amplitude as max(|A|), rather than intensity max(|A| 2 ) for clarity: the resulting intensity difference is notably pronounced.This suggests the presence of substantial intensity-dependent phase modification (i.e., owing to nonlinearity) in the higher-dimensional systems.Consequently, one can anticipate that in these systems reducing the pump amplitude in these systems may enhance the stability of the LS.

V. BIFURCATION ANALYSIS THROUGH DIMENSIONS
Building on the previous findings, we perform a much more detailed bifurcation analysis of the LS for a continuous variation of the dimension parameter.For enhancing the comparison between different dimensions, we represent the diagram of the nonlinear LS shown in Figs.3(a, b, c) using the maximum amplitude max(|A|) in the bifurcation diagram presented in Fig. 4(a).
The LS solutions (i-vi) in Fig. 4(a) are illustrated in Figs.4(i-vi), through the extended r-profiles at δ = 1.2, −4, −9 for 1D and 2D, respectively.Subpanels in Figs.4(iv-vi) depict the 2D reconstructions of such profiles.The (i-vi) 1D and 2D states are stable, while 3D LS are mainly unstable [see Fig. 4(a)].Indeed, higher dimensions exhibit broader unstable and coexistence regions.These regions can be numerically tracked by path-continuing the fold and Hopf bifurcations in D and δ .As a result, we obtained the (D, δ )-phase diagram of Fig. 4(b), which offers important insights into the dynamical behavior of the system across the different dimensions.
In Fig. 4(b), within each FB curve F m , where m = 1, 2, 3, ..., we observe regions of coexistence of different LS, which are generated due to the tilt of the resonance peaks.Apparently, as the dimensionality increases, the system exhibits different numbers of pairs of FB points: one pair for D = 1, two pairs for D = 2, and many more pairs when D = 3.This phenomenon can be understood from the fact that higher dimensions concentrate much stronger field intensities around r = 0, resulting in larger intensity-dependent phase shifts.This, in turn, leads to much more pronounced tilts of the resonance peaks.More importantly, it is possible to map the system stability boundaries by path-continuing the HB points in Fig. 4(a), involving variations in both δ and D. The resulting HB curve, depicted in purple in Fig. 4(b), intersects the horizontal gray line at two points corresponding exactly to the HB points in Fig. 4(a).This curve organizes the LS according to their stability: above this curve, all states in the shallow purple region are unstable.
We perform a comparable bifurcation analysis for a larger pump.Specifically, for P = 3.2, the bifurcation diagrams corresponding to three dimensions are depicted in Fig. 4(c The unstable 2D LS (vii, viii, ix) in Fig. 4(c) correspond to the profiles depicted in Fig. 4(vii-ix).These states exhibit stronger localization around the origin r = 0 compared to those with the same detunings in Fig. 4(iv-vi).The 3D diagram in Fig. 4(c) mirrors the 2D case, but the LS exhibits a significantly higher peak intensity.
Importantly, the phase diagram in Figs.4(d) provides a clearer view, revealing that an increase in pump amplitude results in a more complex dynamic map.The fold bifurcation not only expands its region in detuning δ but also extends further into a larger region in lower dimensions.The unstable region, by the purple HB curves, enlarges its domain towards lower dimensions, reaching even D = 1.
For lower P values, the coexistence and breathing regions shrink, as shown in the bifurcation and phase diagrams of Fig. 4(e, f) for P = 0.75.Here, the stable 3D LS corresponding to positions (x,xi,xii) in Fig. 4(e) are shown in Fig. 4(x,xi,xii), along with their 2D and 3D representations.At such a low pump value, all nonlinear states in 1D and 2D are stable.In Fig. 4(f), only the initial four FB curves intersect at D = 3 whereas the Hopf bifurcation (HB) curves exhibit periodic intersections, creating a periodic alternation between stable and unstable intervals in Fig. 4(e).These results are in-line with our previous results [29].

VI. DISCUSSIONS AND CONCLUSIONS
We would like to emphasize the key advantages, as well as the limitations of the proposed approach.Firstly, considering radially symmetric solutions for reducing high-dimensional systems to lower ones permits to reduce complexity, which permits to gain physical insight, and to improve numerical accuracy thanks to the massively decreases computational time.Secondly, treating the system dimension D as a continuous and controllable parameter allows for a seamless exploration of solutions across their dimension, greatly facilitating the understanding of stability changes when transitioning from 1D to 3D structures.In the context of pattern forming systems, this has been proven to be extremely useful for analyzing the modification LS and their associated homoclinic snaking toward dimensions in the Swift-Hohenberg equation [55,56].
Nevertheless, reducing the problem complexity also leads to some trade-offs.In cases with a relatively weak pump, the radial symmetric simplification yields results which are virtually identical to those obtained in their full dimensional problems.However, as the pump intensity increases, a greater variety of instabilities may emerge.The radial symmetric simplification might remove any possibility of studying less constrained scenarios, involving more morphologically complex states.Consequently, in situations with a relatively high pump, the stability of the LS, derived via pathcontinuation of the radial profiles, may necessitate additional confirmations through DNS or linear stability analysis performed for the original problem across its full coordinates.
As a matter of fact, in systems with translation invariance (i.e., in the absence of the localizing potential), states with different morphology (e.g., localized stripe, hexagons and rhomboid patches) may arise from radial LS via symmetrybreaking bifurcations [48,49], which are absent in our case.This is also intrinsically related to stability properties of LS.
Here we only study how the system is stable against radial perturbations, thus neglecting any other degree of freedom, such as those associated with transverse or azimuthal perturbations, which are closely related to the geometric curvature of the nonlinear states [57][58][59].
To summarize, in this article, we have analyzed the formation of radial symmetric LS of different dimensions in externally driven optical cavities subject to the effect of a parabolic potential.We have transformed higherdimensional models into a single 1D model with a dimension parameter.This transformation provides an additional perspective: compared to the 1D system, 2D and 3D introduce an effective position-dependent field flow.This raised flow rate generates an effective centripetal "force" weighted by dimension parameter D, compelling the field to concentrate at the LS center.Linear eigensolutions reveal heightened concentrations of eigenmodes with consistent equal spacing among eigenvalues as the dimension parameter rises.Furthermore, we have performed a detailed bifurcation analysis of nonlinear solutions by using the dimension parameter D as a continuous variable.This has permitted us to gain essential information about the connection between LS of different dimensions: For the same pumping and detuning system parameters, localized states in higher-dimensional systems exhibit greater field concentration due to the introduced effective centripetal "force".This increased intensity results in broader bistable regions but diminished stable dynamical regions.The results are illustrated in (D, δ ) phase diagrams for different values of the pump strength.As the dimension parameter D or pump P increases, the system displays an expansion of the LS coexistence regions, albeit concurrently encountering a diminished stability range.This study introduces an approach to address high-dimensional problems, illuminating the fundamental dimensional connections among radially symmetric states across varying dimensions, and supplying valuable analytical tools.
In future research, we plan to extend the present work by including a spectral stability analysis of non-radial perturbations.This will provide the detection of symmetry breaking bifurcations, and point to new ways for studying the emergence of LS with non-radially symmetric morphology.

FIG. 1 .
FIG. 1.(a) Total number of sampling points required for physical models in 1D, 2D, and 3D models depend on the sampling points in a single direction.The shallow red region represents computational loads that are extremely heavy for a typical personal computer.Panel (b) shows radial symmetric multidimensional LS in 1D, 2D, and 3D models, featuring highlighted radial profiles represented by red curves.In (c), multidimensional LS in (b) can be analyzed in terms of their radial profile, by leveraging on their symmetry.

FBHBFIG. 3 .
FIG. 3. (a-c) Comparison of bifurcation diagrams of the field norm vs. δ in the absence (gray curve) and in the presence (black solid and red dashed curve) of Kerr nonlinearity for different dimensions: 1D in (a), 2D in (b), and 3D in (c), when P = 1.3.Vertical lines in (a-c) correspond to δ = 1.2.(i-vi) Individual solutions for the field amplitude distribution |A|, corresponding to the solution marked (ivi) in the bifurcation diagrams (a-d).(d) Bifurcation diagram of the maximum amplitude vs. the dimension parameter D, for P = 1.3 and δ = 1.2.
3(b)], two pairs of FB points appear around the first two resonance peaks.Between these peaks, 2D breathers [see the dashed curve] appear.Coexisting (i.e., bistable) LS are shown in Figs.3(iii,iv) for δ = 1.2.A reconstruction of the 2D field A(x, y) is obtained from the radial profiles and shown in Figs.3(iii,iv) inset.

1 FIG. 4 .
FIG. 4. Bifurcation diagrams for max(A) vs. δ and phase diagrams in the (D, δ ) plane for the three parameter values: P = 1.3 in (a,b), P = 3.2 in (c,d), and P = 0.75 in (e,f).In (a,c,e), stable solutions are shown as solid black curves, while unstable solutions are represented by red dashed curves.In (b,d,f), purple curves correspond to Hopf bifurcations (HB), with the remaining curves indicating fold bifurcations (FB) F n .(i-xii) Individual solutions for the field amplitude distribution |A|, corresponding to the solution marked solutions (i-iii) for 1D and (iv-vi) for 2D in diagram (a), (vii-ix) for 2D in diagram (c), and (x-xii) for 3D in diagram (e).
), and its associated (D, δ )-phase diagram is presented in Fig.4(d).In contrast to the previous case illustrated in Figs.4(a), the bifurcation diagram in Figs.4(c) displays a more pronounced tilt, featuring four pairs of FB in 1D and two pairs of HB.In the 2D scenario, a considerable portion of the bifurcation diagram is unstable, except for a confined region spanning approximately the interval −15 < δ < −13.