Solitary wave billiards

PHYSICAL REVIEW E 107, 034217 (2023) Solitary wave billiards Jesús Cuevas-Maraver Grupo de Física No Lineal, Departamento de Física Aplicada I, Universidad de Sevilla, Escuela Politécnica Superior, Calle Virgen de Africa 7, 41011 Sevilla, Spain and Instituto de Matemáticas de la Universidad de Sevilla (IMUS), Edificio Celestino Mutis, Avenida Reina Mercedes s/n, 41012 Sevilla, Spain Panayotis G. Kevrekidis and Hong-Kun Zhang Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-4515, USA (Received 31 March 2022; accepted 8 March 2023; published 31 March 2023) In the present work we explore the concept of solitary wave billiards. That is, instead of a point particle, we examine a solitary wave in an enclosed region and examine its collision with the boundaries and the resulting trajectories in cases which for particle billiards are known to be integrable and for cases that are known to be chaotic. A principal conclusion is that solitary wave billiards are generically found to be chaotic even in cases where the classical particle billiards are integrable. However, the degree of resulting chaoticity depends on the particle speed and on the properties of the potential. Furthermore, the nature of the scattering of the deformable solitary wave particle is elucidated on the basis of a negative Goos-Hänchen effect which, in addition to a trajectory shift, also results in an effective shrinkage of the billiard domain. DOI: 10.1103/PhysRevE.107.034217 I. INTRODUCTION In recent years, the study of solitary waves has been a topic of widespread appeal in a broad range of fields. For instance, relevant coherent structures arise in the exploration of electrical field dynamics within optical fibers [1,2], as well as in the study of nonlinear effects in plasmas [3]. They emerge in atomic clouds of Bose-Einstein condensates (BECs) [4,5], as well as in so-called rogue waves in the ocean [6]. Of central role within most of the above studies and fields has been the envelope wave model of the nonlinear Schrödinger equation [7–10], which can be used to describe the electric field of light, the wave function of BEC atoms, or the water wave elevation. Indeed, the resulting nonlinearity may stem from different sources, such as the so-called Kerr effect in optics [1,2] or an effective mean-field nonlinearity due to contact interactions between bosons [4,5], yet its impact is similar across fields in producing robust solitary waves of either bright [2,10], dark [11], or rogue [6]form. On the other hand, the notion of billiards is one that has attracted significant attention both at the classical, as well as at the quantum side. Here, a point particle is reflected within a (typically) enclosed two-dimensional (2D) domain Q(of different possible shapes) [12,13](seealso[14] for an extensive set of references). The particle moves freely aside from its interaction with the boundary ∂Q, by elastic reflections without loss of speed. Depending on the nature of the billiard (e.g., square or circular or instead an elliptical stadium [15,16] or a square with an enclosed circle [12,17]), the outcome of the reflections could be a closed trajectory associated with integrable motion, or an ergodic trajectory associated with a chaotic billiard. While the above aspects involve classical billiards, we do note in passing the consideration of quantum billiards. The quantum mechanical billiard is given by the linear Schrödinger equation typically with Dirichlet boundary condition. It describes a wave function associated with a probability density |ψ|2within the domain while the wall features an infinite potential wall. A relevant summary of some of the activity in the context of experiments (and theory) in the theme of microwave billiards and the connections to quantum chaos can be found in [18], where a number of references to this theme are summarized. The dynamics of a quantum billiard are determined by the Hamiltonian equations of motion on the domain Q. The wave function ψof a single quantum particle with unit mass obeys the Schrödinger equation i∂ψ ∂t=−∇ 2ψ+V(r)ψ, r=(x,y)∈int(Q)(1) with, e.g., Dirichlet boundary conditions ψ(x,y)=0, for (x,y)∈∂Q. Usually one can assume in such a setting that V(x,y)=0for(x,y)∈int(Q), and V(x,y)=∞,for(x,y)∈ ∂Q. Most research has been concentrated, for such quantum billiards, on the search for quantum chaos. Even though it was proved that if a classical billiard is ergodic, then the corresponding quantum billiard is quantum ergodic, socalled quantum “scars” are observed in most cases, which demonstrate that quantum billiards are more complicated. The Bunimovich stadium [15,16] has recently become a very popular quantum billiard model, while the study of quantum chaos in different billiards is generally a topic that is gaining considerable traction [19–21]. A less studied billiard is the “soft” version, where the boundary of the billiard table Qis defined by a finite potential 2470-0045/2023/107(3)/034217(11) 034217-1 ©2023 American Physical Society CUEVAS-MARAVER, KEVREKIDIS, AND ZHANG PHYSICAL REVIEW E 107, 034217 (2023) wall V, with V(x,y)=0for(x,y)∈int(Q) and V(x,y)∈ (0,∞)for(x,y)∈∂Q. In such a setting, for a Sinai (dispersing) billiard [12,17], the existence of a stable island around a periodic orbit which is tangent to the billiard’s wall (or near the corner) was found numerically in [22]; see also [23,24]for other examples of such soft billiards. In the present work our aim is to explore a notion that, in a sense, interweaves the above concepts while paving an avenue involving the interplay of nonlinearity, billiards and possibly quantum mechanics (but also perhaps more concretely in our examples below, nonlinear optics) in the form of solitary wave billiards. In particular, the connections are quite natural: Solitary waves are typically thought of as effective particles (indeed, hence the name) [2,7–11]. Therefore, it is natural to examine these effective particles in a billiard setting. At the same time, nonlinear Schrödinger equation models appear in a natural sense to generalize the concept of linear quantum billiards in the presence of different types of nonlinearity, stemming from nonlinear optics [2], plasmas [3], or atomic BECs [11]. In that light too, it is useful to adapt some of the ideas of linear quantum billiards to such more complex cases where the localized wave function is held together because of the interplay of nonlinearity and dispersion while slowly losing some of its initial kinetic energy through the (potential) interactions with the wall. In that vein, the notion of a soft billiard is a natural one as the solitary wave interacts with a finite potential wall. However, the interaction may be inelastic even in the context of an infinite potential due to the nonintegrability of the partial differential equation (PDE) model. The latter is the typical scenario in the associated 2 +1-dimensional spatiotemporal models. It is worthwhile to mention here that the idea of billiards in settings such as BECs has not only been theoretically explored [25,26], but also experimentally realized in the form of an optically induced confinement potential (i.e., billiard) for the ultracold atoms [27]. However, exploring such solitary wave billiard ideas requires some extra care. One needs to avoid the hurdle of collapse (arising, e.g., in 2 +1-dimensional cubic nonlinearity models [10]). Instead one needs to consider a nonlinearity that allows for a robust 2 +1-dimensional standing wave; here the saturable nonlinearity of relevance to photorefractive crystals is thus selected [28] as a physically relevant example. While in such systems pattern formation has been much studied earlier [29], here we have in mind a different setting well summarized in the authoritative review of [30]. In particular, as is explained therein (see, e.g., Fig. 3.4 and associated discussion), optical beams of interfering plane waves in the ordinary polarization cause an effective crystal lattice potential for the extraordinary polarization. Rather, our proposal (which, while not realized to our knowledge, seems well within the reach of such experiments) is to utilize suitable beams in the context of the ordinary polarization in order to form, e.g., the billiard configurations proposed below. Then, the solitary wave can be boosted using a Galilean transformation and the billiard dynamics will be accordingly initiated. In what follows, we discuss our numerical observations and an emerging qualitative understanding of the relevant phenomenology. We also note in passing that, to our knowledge, no Hamiltonian variant of this type has been considered in the literature. A dissipative (but occurring over a narrow region of the associated parameter space) analog of such concepts has been presented in [31]; see also [32]. A pioneering effort in a related direction involved the experimental demonstration of the interaction of the so-called walkers (a droplet bouncing on vibrating bath) and their interactions with barriers, as shown, e.g., in the work of [33]. Indeed, the associated field of hydrodynamic quantum analogs has seen a wide range of developments that have now been summarized in reviews such as [34]. As we will see below, the present setting will offer a family of (stationary) solitary waves in terms of their amplitude and/or width as we will showcase below, but also the freedom to boost these solitary waves at essentially arbitrary speeds based on the above mentioned Galilean transformation. We consider these features to be key to the potential interest in the setting that we propose herein. Our presentation of the results is structured as follows. In Sec. II, we present the mathematical setup of the problem. In Sec. III, we discuss the numerical observations and their comparison to the point billiard model. Finally, in Sec. IV, we summarize our findings and present some conclusions and directions for future study. II. MATHEMATICAL SETUP A. Soliton wave billiard We start from a two-dimensional nonlinear Schrödinger (NLS) equation with (saturable) photorefractive nonlinearity [28], iψt=L(ψ),L(ψ)=−∇ 2ψ+V(r)ψ−2|ψ|2 1+|ψ|2ψ, r=(x,y)∈Q(2) with V(r) being the potential where the relevant geometry of the billiard will be introduced on Q, and ∇2ψ=ψxx +ψyy. Here, ψ(r,t) represents the envelope of the electromagnetic field for the optical problem of interest, with the density ρ= |ψ|2corresponding to the light intensity. While the particular nonlinearity is less directly relevant to atomic BECs, one can envision other nonlinearities there (e.g., competing ones as in the recently budding field of quantum droplets [35]) that can lead to a similar phenomenology. In the latter case ψ(r,t) plays the role of the condensate wave function, while |ψ|2 represents the atomic density. In the optical setting the potential V(r) stems from the refractive index profile [2], while in BECs the ability to “paint” arbitrary (and, hence, also billiard) potentials has been demonstrated [36]. The above dynamical equation can be derived from a Hamiltonian formulation as iψt=δH[ψ] δψ∗,(3) where H[ψ] is the (conserved during the dynamics) Hamiltonian functional: H[ψ]=r∈Q [|∇ψ|2+(V(r)−2)|ψ|2 +2log(1+|ψ|2)]dr.(4) 034217-2 SOLITARY WAVE BILLIARDS PHYSICAL REVIEW E 107, 034217 (2023) The initial condition for our dynamical simulations will be a stationary solitary wave set into motion by means of a Galilean boost. This is a desirable feature of the present model as the traveling waves herein can be generically achieved and with arbitrarily selected speeds due to the Galilean invariance. The stationary state is found by introducing the ansatz ψ(x,y;t)=exp(iωt)φω(x,y)(5) into Eq. (2), and a fixed point algorithm along with a finitedifference discretization has been used to this aim. It follows that φω, if it exists, is a nonlinear eigenfunction of F: −ωφω=F(φω)=−∇ 2φω+V(x,y)φω−2|φω|2 1+|φω|2φω. (6) Once this stationary wave is found, we fix ω, and the initial condition of the simulation is attained from the Galilean boost of the stationary structure centered at (x,y)=(0,0), ψ0(x,y)≡ψ(x,y;0)=φω(x,y)exp[i(kxx+kyy)] =φω(x,y)exp[ik(xcos θ+ysin θ)],(7) leading to a wave moving with an initial velocity v0=2|k| that forms an angle θwith the xaxis, and k=(kx,ky)=|k|(cos θ,sin θ)=k(cos θ,sin θ).(8) Consequently, the energy imparted to the moving coherent structure (in comparison to its stationary state) is given by E≡E[ψ0(x,y)] −E[φω(x,y)] =d2r[|∇ψ0|2−|∇φω|2]=N|k|2(9) and is naturally associated with the excess kinetic energy, with Nbeing the squared L2norm or the mass, as it is often referred to, of the stationary state N=Q |φω(r)|2dr.(10) Assigning N≡2mto an effective mass, we can express this excess energy as E=2m|k|2=1 2mv2 0.(11) If we suppose that the solitary wave is moving in a vanishing potential landscape, it can be considered as a quasiparticle that freely moves with an initial kinetic energy equal to E, and this quasiparticle can be interpreted as possessing an effective mass m. Note that at each time t⩾0, the solitonic wave function ψ(x,y,t), which starts from ψ0(x,y), defines a reference probability measure ν(t)ontheσalgebra of the billiard table Q, with density function ρ(r,t)=|ψ(r,t)|2 |ψ(r,t)|2dr. Using this reference measure, we can investigate some physical quantities (given by the expected position, variance, and more general moments) of the solitary wave billiard. To study this billiard ψ(x,y,t), we first consider its center of mass trajectory γ(t), which describes the expected position of the coherent structure. More precisely, the expected position of the localized pattern with respect to the reference measure ν is defined as γ(t):=Eνt(r)=Q (x,y)ρ(r,t)dr. Since r=(x,y), we can also write γ(t)=(X(t),Y(t)) with X(t)=Eνt(x)=xdνt,Y(t)=Eνt(y)=ydνt, as the expected xand yposition of the solitary wave, respectively. Note that (X(t),Y(t))can also be interpreted as the center position of the relevant pattern. In the following sections, we will mainly focus on the study of r(t)=(X(t),Y(t)) numerically. To assess how well the expected position of the solitary wave within the billiard describes its motion, we also calculate the variance of rwith respect to the reference measure ν(t): Varνt(r):=Eνt(r−γ(t)2) =r∈Q ([x−X(t)]2+[y−Y(t)]2)ρ(r,t)dr This quantity characterizes the width of our solitary wave. For narrower solitary waves, it is enough to study the center of mass orbit γ(t), corresponding to the solitonic wave function ψ(x,y,t). III. NUMERICAL RESULTS A. Numerical simulation setup of the dynamical orbit within the solitary wave billiard In our simulations, to simplify the calculations, we first locate the solitary wave at the point of maximal probability density ρ(r,t), which is denoted as ˜γ(t)=(˜x(t),˜y(t))such that, for any t⩾0, ρ(˜γ(t))=max r∈Qρ(r,t). Our computations suggest that the curve ρ(˜γ(t))is well defined, i.e., for each t, there is only one (principal, i.e., above a selected cutoff) local maximum for the density function. We subsequently consider a window of size 2δ×2δcentered on ˜γ(t), and denote it as Uδ(t)⊂Q; this is suitably adjusted when the wave is in the vicinity of the boundary. In our simulation, we select δ=4. The conditional probability measure ˆνt:=ν|Uδ(t)with density function ˆρcan be calculated as ˆρ(r,t)=|ψ(r,t)|2IUδ(t) Uδ(t)|ψ(r,t)|2dr. Moreover, the simulated expected position of the solitary wave within the billiard with respect to the conditional measure ˆνcan be calculated as (ˆ X(t),ˆ Y(t)):=Eˆνt(r)=Q (x,y)ˆρ(r,t)dr. This also provides a smoother and more accurate evaluation of the center of mass, in comparison to the argmax over the numerical grid. Furthermore, the variance associated with the 034217-3 CUEVAS-MARAVER, KEVREKIDIS, AND ZHANG PHYSICAL REVIEW E 107, 034217 (2023) FIG. 1. Dependence of the norm and the width of the solitary wave with respect to its frequency for the relevant family of solutions. solitary wave position is given by W(t):=r∈Q[x−ˆ X(t)]2+[y−ˆ Y(t)]2ˆρ(r,t)dr. In Fig. 1, one can see that the narrowest solitary wave is found for ω=0.95, and that the associated frequency lies in ω∈(0,2) interval, with the lower limit corresponding to the bifurcation of the solitary wave from the linear modes (continuous spectrum) band, while the upper limit reflects a divergence of its mass. Indeed, the continuous spectrum of the problem lies in the interval (−∞,0), hence the avoidance of resonances therewith leads to solitary waves with ω>0. On the other hand, the large intensity limit of the nonlinearity leads the nonlinear term to the asymptotic value of −2, hence Eq. (6) accordingly becomes (in the large intensity limit) (ω−2)φω=∇ 2φω, which requires ω<2. As a result, solutions are expected to be found for 0 <ω<2 in line with our numerical results. In our study, we mainly consider waves with frequency ω=0.5. In the case examples that follow, we will either consider infinite potentials (implemented via Dirichlet boundary conditions) or finite potentials that are strictly higher than the above-mentioned kinetic energy of the solitary wave, so as to ensure a genuine billiard situation in which the coherent structure will be reflected from the potential walls. B. Goos-Hänchen shifts of reflection upon potential walls in billiards of constant potential Before delving into the systematics of the solitary wave trajectories within the entire billiard, we briefly highlight a key feature of the billiard reflection of the wave upon collision with the potential walls, for the case when the potential Vis finite at the boundary of Q. When describing total reflection of a solitary wave in such a billiard potential setting, we need to identify expressions for the phase shift that occurs for each such event between the incident and reflected waves as a function of angle of incidence. A way to describe the phase shift is as being due to the orbit actually traveling a small distance into the potential wall before being reflected. An incident orbit behaves as if it were laterally displaced upon reflection. This feature has been previously recognized as being analogous to the so-called Goos-Hänchen shift in optics, occurring also to solitary waves during their interaction with an interface; see, e.g., Ref. [37] and the discussion therein. We will discuss this in further quantitative detail in what follows below. It is interesting to note here that the Goos-Hänchen effect (GHE) has been recently encountered not only in the context of solitary waves interacting with external potentials but also of linear waves interacting with solitary ones [38]. We first consider a solitary wave placed in a billiard involving a square potential. It is known that classical square billiards are completely integrable [39]. Moreover, the orbits are periodic if tan θis rational, for the collision angle θ. When studying the evolution of solitary waves, we first considered a square barrier potential of amplitude α, i.e., V(x,y)=αif |x|>3L/4or|y|>3L/4, 0 otherwise. (12) The domain is taken with L=40. This potential is shown on the left panel of Fig. 2. The evolution of dynamical equations has been found by means of numerical integration using afixedstepδt=0.0625 fourth-order Runge–Kutta algorithm with a spatial discretization δx=0.5. The latter value is reasonable enough as the soliton width is several times that value. In fact, we have used smaller values such as 0.2 and 0.1 and have corroborated that the results do not change in any observable manner (although the latter cases utilize far larger computational resources and hence are considerably slower, hence the selection of δxused for our presentation herein). We have also checked that, for this choice of δt and δx, the relative error for the norm, i.e., |N(t)/N(0) −1|, extending the definition of Eq. (10) to time dependent states as N(t)=|(r,t)|2dr, does not grow above 10−6for the time horizons considered. In what follows, we have chosen cos θ=3/5 for the initial collision angle. The right panel of Fig. 2depicts the path of a classical point particle in this square potential. One can indeed see, in line with the theoretical expectation, that the dynamics is periodic for this choice of θ. On the other hand, in the left panel of Fig. 3, the path of the solitary wave with v0=0.02 and α=4 is considered. One can see that, in contrast to the classical billiard case, the dynamics is not periodic. In addition, the finite size of the solitary wave particle does not allow its center to reach the edge of the potential barrier 034217-4 SOLITARY WAVE BILLIARDS PHYSICAL REVIEW E 107, 034217 (2023) FIG. 2. Left: Form of the square potential barrier given by (12). Right: The trajectory of a classical point particle (right), in this potential. (at ±3L/4 in each direction), but rather it effectively gets reflected “sooner.” This is caused by the reflection rule of the trajectory of the solitary wave upon a constant potential (similar to what is observed in [31]). Nevertheless, it is important to recognize that the solitary wave is not an infinitesimal particle mass. Instead, it features a finite width, and for that reason this effective particle is reflected while its center of mass is still at a considerable distance from the domain boundary. There is some tunability of the width of the solitary wave (controlled by ω), however, the constraint of avoiding the cubic nonlinearity and associated collapse phenomenology does not allow for the freedom of (effectively) arbitrarily selecting the solitonic width, as, e.g., can be done in one-dimensional such systems. Rather, here we are practically considering not only a soft potential (unless we examine the Dirichlet boundary condition case), but also a nonvanishing size particle that features a finite width and can, in principle, manifest internal modes capable of storing energy. This “deformability” of the solitary wave can also be detected in Fig. 3at the positions of wall interactions where the trajectory of the solitary wave curves around in the reflection process, a manifestation of the GHE discussed above (see also details below). As an interesting additional case example along this vein, in the right panel of Fig. 3, we have considered the example of the structure with the smallest width, namely ω=0.95. Interestingly, as we can see in that case the trajectory does not close, indicating the subtle nature of the corresponding phenomenology. It can be observed that the dynamical evolution does not purely depend on the angle of incidence, but also on other features of the soft particle, such as its width, its potential excitation of internal modes (or dispersive radiation modes, as will be seen below), etc. This can also be confirmed in the consideration of settings with different α(such as α= 10) or different boundary conditions (such as Dirichlet), for FIG. 3. Evolution of the center-of-mass of the finite-width solitary wave with α=4andv0=0.02, for ω=0.5(left)andω=0.95 (right). The final time of both simulations is t=2×105. Notice that, in this figure and the other ones showing square barriers, the axis limits correspond to the barrier boundaries. See also companion movies movie1.gif (corresponding to ω=0.5) and movie2.gif (corresponding to ω=0.95) in the Supplemental Material [40] depicting the evolution of the solitary wave dynamics, in which the whole domain is shown with the barrier limits being depicted as white lines. 034217-5 CUEVAS-MARAVER, KEVREKIDIS, AND ZHANG PHYSICAL REVIEW E 107, 034217 (2023) FIG. 4. Evolution of the center of mass of the solitary waves with v0=0.1 (left panel) and v0=0.8 (right panel), in the square barrier potential (12). In both cases, α=4andω=0.5. The final times of the simulations are t=4×104and t=5×103, respectively. See also companion movies movie3.gif and movie4.gif in the Supplemental Material [40] depicting the evolution of the solitary wave dynamics. which, as we have checked, the soliton follows a very similar trajectory to that of the α=4 case. We have also explored how the above dynamics is altered when the velocity is increased and the barrier height is kept fixed; see, e.g., left panel of Fig. 4for the solitary wave with v0=0.1 and α=4. The resulting motion in this case as well is not periodic. If the wave moves more rapidly (v0=0.8, as shown in the right panel of Fig. 4), it deforms nontrivially in the collisions with the barrier, and some radiation is emitted. The evolution of the coherent structure motion can be observed in more detail in the companion movies (Supplemental Material [40]), where the density (|ψ(x,y;t)|2) and its logarithm is tracked. One can see how the solitary wave approaches the barrier and is reflected after each collision. This process is central to our observation of the trajectories of the solitary waves (and their nonclosed form), hence we focus on this further in what follows. However, before doing so, we point out the substantial difference between this high speed case example and the ones shown earlier, as concerns the substantial emission of dispersive wave radiation. This results in the solitary wave moving within a “turbulent” small amplitude background, which, in turn, weakly affects (via scattering processes) the solitary wave motion. While this aspect is not pursued further herein, we find it an interesting direction for further exploration in the future. Let us remark that in this particular case, given the faster growth of the error in the fourth-order Runge-Kutta integrator, we have numerically integrated the evolution by means of an adaptive step-size Dormand-Prince algorithm, keeping the relative error below 10−7during the whole integration time. An explanation for the fact that the wave’s center of mass does not reach the barrier is based on the intriguing observation of a negative Goos-Hänchen effect (see, e.g., [41]for such an effect in one-dimensional optical systems, as well as [37] for a solitonic example). Figure 5shows the evolution of the integrated norm with respect to the xcoordinate, i.e. |ψx(y,t)|2=dx|ψ(x,y,t)|2together with a scheme defining the Goos-Hänchen shift (GHs), , for the collision of the solitary wave with a wall barrier of height α. As is well known, the GHE occurs when a beam of finite extent is incident on a medium with smaller refraction index with an angle higher than the critical angle for total reflection [recall that −V(x,y) plays the role of refraction index detuning in optics]. The shift is caused by the field penetrating the medium for a small distance, forming a nonuniform plane wave that is evanescent in the direction normal to the interface and propagating also along it. The shift is positive in the usual case of plane waves flowing parallel to the interface in the forward direction. There are, however, cases when the wave reflects from a point in front of the contact point in the geometric optics picture; in these cases, observed for plane waves at interfaces with metals and predicted for metamaterials, the plane waves flow parallel to the interface in the backward direction and the GHE is negative; in addition, part of the beam can cross the interface despite the negative GHE, as is the case in metamaterials. All the aforementioned behavior holds for linear waves. In the case of nonlinear waves, in order to provide an analogy with geometric optics, one can trace straight lines tangent to the trajectory of the soliton’s center of mass at the initial point, in a fashion similar to what is done in Fig. 5. Notice that, because of the GHs, the motion of the wave’s center of mass when interacting with the wall can be approximated by a classical particle in a domain shrunk by a distance σ=(cot θ)/2. An important point is that in the square barrier, the impact angle of the solitary wave with the vertical sides of the square is complementary to the angle of impact with the horizontal sides. Consequently, there will be two different GH shifts to account for, namely the one corresponding to an initial angle θ(denoted as ) and the one for θc=π/2−θ[denoted as c, with the corresponding shrinkage σc=(ctan θ)/2] in the perpendicular direction. Then, the trajectory of the center of mass of the solitary wave can be approximated by a classical particle in a rectangle bounded by |x|=3L/4−σcand |y|=3L/4−σ. Figure 6 shows the numerically computed dependence of the domain shrinking caused by the GHE with respect to v0for fixed ω=0.5, and also the dependence with ωfor v0=0.1. One can see that, in the latter case, the minimum values of σand σc 034217-6 SOLITARY WAVE BILLIARDS PHYSICAL REVIEW E 107, 034217 (2023) FIG. 5. Left panel: Evolution of the projected solitary wave density when interacting with a wall barrier of height α=4 located at y=0; the initial velocity of the coherent structure is v0=0.1 and the initial angle θ=cos−13/5. Right panel: Path followed by the center of mass of the wave together with a scheme indicating some relevant parameters like the Goos-Hänchen shift or the domain shrinking σ. are not attained for the narrowest solitary wave, i.e., that with ω=0.95, but rather around ω=1.05 (the particular value depends on the barrier height). The above description based on the GHE enables a systematic understanding of the solitary wave interaction with a square billiard. Figure 7shows the evolution of a classical particle in a shrunk domain by the GHs corresponding to a solitary wave in the presence of a barrier with α=4. This evolution is quantitatively compared with that of the center of mass of the coherent structure. Notice that, until the second collision, the two trajectories overlap; the separation between them starts to subsequently appear because the wave interacts without following a perfectly straight line, but rather has a curved part in its trajectory, culminating in the GHE and resulting in an evolution of the particle in an effectively shrunk domain. This weak effect compounds itself over time, eventually leading to the nonclosure of the orbit. We have also checked that the relevant situation does not improve for larger α, or for Dirichlet boundary conditions. While this procedure allows for a highly accurate description of individual collisions, slight deviations from straight motion of our deformable effective particles ultimate lead to the nonclosed trajectory, generic phenomenology observed herein. FIG. 6. Dependence of the domain shrinking caused by the Goos-Hänchen effect for θ=cos−1(3/5) (left panels) and θ=cos−1(4/5) (right panels) for fixed frequency ω=0.5 (top panels) and fixed initial speed v0=0.1 (bottom panels). Dirichlet boundary conditions are denoted by α→∞. 034217-7 CUEVAS-MARAVER, KEVREKIDIS, AND ZHANG PHYSICAL REVIEW E 107, 034217 (2023) FIG. 7. Trajectory of a classical particle in a shrunk domain with a Goos-Hänchen shift corresponding to a solitary wave in a square barrier potential of height α=4 with the domain shrinking corresponding to v0=0.02 and ω=0.5. The left panel shows the evolution until t=2×105, whereas right panel captures the early stage up to t=34880 and compares it with the evolution of the the center of mass of the solitary wave (full blue line). One can observe that up to the second collision with the barrier, the trajectories overlap, while at later times, they slightly deviate from each other. The fact that the domain shrinking caused by the GHE is different for θand θcleads the solitary waves billiard to behave effectively as if it is in a rectangular domain. In that case, closed orbits in classical billiards can only occur when (a/b)tanθ∈Q, with aand bbeing the rectangle sides. In our case, if tan θ∈Qand supposing that the solitary waves billiard behaves as a classical particle in a shrunk domain of sides 3L/2−2σand 3L/2−2σc, closed orbits could only be found in some isolated points within the parameter space. However, there is a case for the classical billiard where closed orbits can occur in the shrunk domain, namely this in which θ=θc. This is possible when θ=45◦. Because of this, we have launched a soliton with this initial angle, but displaced from the domain center; in particular, we have initialized the soliton from the point (−3L/8,3L/8)whichwouldgive rise to a square trajectory in the case of a classical particle. Figure 8shows the trajectory for a soliton launched with an initial velocity v0=0.1. This trajectory is not a square because the initial point is not at the center of the second quadrant of the shrunk domain. Notice that the closed trajectory can be found for v00.25, as for larger values the radiation caused by the soliton impacting on the barrier prevents from the trajectory closure. C. Some other prototypical solitary wave billiards We end this section by showing the outcome of the solitary wave dynamics in a Bunimovich stadium and a Sinai billiard (see the left top and bottom panels, respectively, of Fig. 9). Here, we are motivated by the chaotic nature of these billiards even for classical point particles and we seek to observe the dynamics of our deformable effective particle therein. In the first case, the barrier (composed by a square whose sides have a length 48 attached to two semicircles of radius 28) has a potential height of α=4 and is embedded into a domain [−64.5,64.5] ×[−32.25,32.25], and the solitary wave is launched from (0,0). In the second case, the domain size is [−75,75] ×[−75,75], encompassing a square barrier of side length 112.5 and a circular barrier of radius 2, both centered at (0,0), that have a potential height of α=4. Here, the solitary wave is launched from the point (−20,−20). The path followed by the center of mass of the wave of initial velocity v0=0.02 is depicted in the right panels of Fig. 9.In both cases, we can observe the chaotic nature of the dynamics which is reminiscent of the ergodic nature of the point particle case, yet again with a significant difference borne out of the GHE present in the solitonic case. The latter effect naturally, in this case as well, precludes the accessibility of a substantial region of the configuration space, similarly to the above presented square billiard case. IV. CONCLUSIONS AND FUTURE CHALLENGES In the present work, we have proposed a new paradigm for the dynamics of solitary waves, namely the study of solitonic billiards. We have judiciously selected a model that, on the one hand, is physically relevant, yet on the other hand avoids well-known pathologies of higher-dimensional NLS models, such as the presence of self-focusing and wave collapse [10]. By choosing the saturable nonlinearity model, we present a setting of relevance, e.g., to photorefractive optical crystals, that possesses a Hamiltonian nature (in line with classical point-particle billiards), but which also has other key advantages, such as the ability to Galilean boost the coherent structures. For this model, we identify the full branch of fundamental solitary wave solutions as the frequency of the standing waves is varied. We identified some important differences that this solitonic billiard features in comparison to the point particle case. Indeed, the finite width of the solitary particle and the nonintegrable nature of the model present new possibilities including that of the collision storing some of the kinetic energy into internal mode oscillations, as well as that of the collision not being perfectly elastic. Both of these features have been observed in direct numerical simulations with 034217-8 SOLITARY WAVE BILLIARDS PHYSICAL REVIEW E 107, 034217 (2023) FIG. 8. Left: Evolution of the center-of-mass of the finite-width solitary wave with α=4, v0=0.1andω=0.5 launched from (−15,15) with an initial angle θ=45◦. See companion movie movie5.gif in the Supplemental Material [40], in which the launching point is marked by a white dot in order to observe that the solitary wave comes back to that point. Right: Trajectory of a classical point particle launched from the same site and initial angle. sufficiently high speeds and/or involving sufficiently wide solitary waves. At the same time, we have observed the presence of a solitonic analogue of the (negative) Goos-Hänchen effect which also leads to significant deviations from the point particle case, including the effective shrinkage of the billiard domain. These features combined render even integrable FIG. 9. Left panels: Form of the Bunimovich stadium and Sinai billiard potentials. Right panels: Evolution of the center of mass of the solitary wave with v0=0.02 in a Bunimovich stadium (left panel) and a Sinai billiard (right panel). See also companion movies movie6.gif and movie7.gif in [40] depicting the evolution of the billiard dynamics. The final time of both simulations is t=3.2×105. 034217-9