J. Sánchez-Curto 2 and P. Chamorro-Posada 2

We consider arbitrary angle interactions between spatial solitons and the planar boundary between two optical materials with a single power-law nonlinear refractive index. Extensive analysis has uncovered a wide range of new qualitative phenomena in non-Kerr regimes. A universal Helmholtz-Snell law describing soliton refraction is derived using exact solutions to the governing equation as a nonlinear basis. New predictions are tested through exhaustive computations, which have uncovered substantially enhanced Goos-Hanchen shifts at some non-Kerr interfaces. Helmholtz nonlinear surface waves are analyzed theoretically, and their stability properties are investigated numerically for the first time. Interactions between surface waves and obliquely incident solitons are also considered. Novel solution behaviours have been uncovered, which depend upon a complex interplay between incidence angle, medium mismatch parameters, and the power-law nonlinearity exponent.


Introduction
A light beam impinging on the interface between two dissimilar dielectric materials is a fundamental optical geometry [1][2][3][4][5][6][7][8][9][10][11][12].After all, the single-interface configuration is an elemental structure that facilitates more sophisticated device designs and architectures for a diverse range of photonic applications.The seminal work of Aceves et al. [6,7] some two decades ago considered perhaps the simplest scenario, where a spatial soliton (i.e., a self-trapped and self-stabilizing optical beam) is incident on the boundary between two different Kerr-type materials.Their intuitive approach reduced the full complexity of the electromagnetic interface problem to something far more tractable, namely, the solution a scalar equation of the inhomogeneous nonlinear Schrödinger (NLS) type.The development of an equivalent-particle theory [3][4][5][6] provided an enormous level of insight into the behaviour of scalar solitons at material boundaries.The adiabatic perturbation technique developed by Aliev et al. [13,14] provides another toolbox for analyzing interface phenomena (e.g., light incident on the boundary between a linear and a nonlinear medium).Photorefractive [15] and quadratic [16] materials have also been considered.
A recurrent feature of the waves at interfaces literature is the appearance of the paraxial approximation, which combines the assumptions of broad (predominantly transversepolarized) beams and slowly varying envelopes [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].The adoption of this ubiquitous mathematical device can impose some strong physical constraints that should be borne in mind when modelling precisely these types of angular geometries.Indeed, the class of problem at hand is inherently nonparaxial since impinging beams may be arbitrarily oblique with respect to the interface.External refraction (where the refracted beam deviates away from the interface) provides a specific context where beam refraction cannot be described using conventional approaches.Paraxial wave optics must be applied with care since, in potentially off-axis regimes, it holds true only where angles (in the laboratory frame) of incidence, reflection, and refraction with respect to the reference direction are negligibly (or nearnegligibly) small.
Recently, we proposed the first scalar model of spatial solitons at interfaces that is valid across the entire angular range [17,18].By respecting the essential role played by Helmholtz diffraction [19][20][21][22][23][24], the angular restriction was lifted while retaining an intuitive and manageable envelope equation.Preliminary analyses considered bright [17,18] and dark [25,26] spatial solitons incident on the boundary between dissimilar Kerr-type materials.They focused on establishing and developing the propagation aspects of our Helmholtz interfaces approach.By enforcing appropriate continuity conditions at the interface, a Snell's law for Kerr spatial solitons was derived whose validity was tested and confirmed by extensive numerical computations.Here, we take the first steps in a systematic study of the materials aspects of nonlinear beam-interface interactions.The simplest non-Kerr system one might consider is a class of host media whose refractive index n NL (E) has a generic powerlaw dependence on the (complex) electric field amplitude E [27][28][29]: where α is a positive coefficient, n 0 is the linear index (at the optical frequency), and the exponent q lies within the range 0 < q < 4. Typically, the nonlinear response of the medium is assumed to be weak so that αE q 0 /n 0 O (1), where E 0 is the peak field amplitude.
Power-law models have played a key role in the theory of nonlinear waves for the past three decades [30,31].Indeed, [32] provides an excellent review of the fundamental importance of model (1) in photonics contexts.Materials that fall into this broad category include some semiconductors (e.g., InSb [33] and GaAs/GaAlAs [34]), doped filter glasses (e.g., CsS x Se x−2 [35,36]), and liquid crystals [32].One expects non-Kerr regimes (where q deviates from the value of 2) to give rise to a diverse range of new quantitative and qualitative phenomena.The physical basis for this expectation lies in the idealized nature of the Kerr response.In a range of materials, one can often find higher-order nonlinear effects coming into play.Perhaps the most obvious example of the breakdown of Kerr-type behaviour is optical saturation, where the refractive index change becomes bleached in the presence of sufficiently high-intensity illumination.In such cases, model (1) with q < 1 can be used to describe generic leadingorder corrections from a saturable (dispersive) nonlinearity [35,36].
In this paper, a detailed account is presented of arbitraryangle refraction of spatial solitons at the interface between dissimilar power-law materials.Also of intrinsic interest are nonlinear surface waves (i.e., localized modes that travel along the boundary).This fundamental class of excitation has been the subject of previous power-law studies involving a single interface [35][36][37][38][39] and guided waves in multilayer structures (e.g., slab waveguides) [40][41][42][43].Stability characteristics have been inferred from inspection of powerpropagation constant solution branches.However, to the best of our knowledge, direct verification of such predictions [37][38][39][40][41][42][43] (e.g., through numerical solution of the underlying nonlinear Helmholtz equation) has been absent from the literature to date.Rather, computational studies of surface waves tend to have been in the limit of slowly varying envelopes and nonlinear Schrödinger-type models, typically of the diffusive-Kerr [44,45], thermal [46], or saturable [47] type.Here, we investigate the stability of exact analytical Helmholtz surface waves through direct numerical calculations.As a fairly stringent test of solution robustness, we also report on some key findings concerning arbitraryangle interactions between surface waves and solitons.In beam-refraction and surface-wave contexts, simulations have uncovered strikingly distinct behaviours as the exponent q is varied and across a range of quasi-paraxial and fully nonparaxial angular regimes.
The layout of this paper is as follows.In Section 2, we propose a governing equation for scalar optical fields in two adjoining power-law materials with dissimilar medium coefficients.Exact analytical bright solitons are presented for both media, and these solutions are used as a nonlinear basis to derive a generalized Helmholtz-Snell law.In Section 3, extensive computations test predictions of the new refraction law over a range of system parameters.We also extend our first calculations of the Goos-Hänchen (GH) shifts [48] in the Helmholtz angular regime [49] with powerlaw nonlinearities.Nonlinear surface waves are derived in Section 4, and simulations provide what appears to be the first full investigation of the stability properties of this new class of Helmholtz solution.We conclude, in Section 5, with some comments about the impact of our results.

Helmholtz Model for Scalar Soliton Refraction
The formalism of Helmholtz soliton theory [23,24] is now deployed to develop a framework for describing refraction phenomena in wider classes of nonlinear optical materials.This type of modelling approach is valid when the beam waist w 0 is much broader than its free-space carrier wavelength λ, such that ε ≡ λ/w 0 O (1).Ultranarrow beam corrections to the governing equation, typically obtained from single-parameter (i.e., ε-based) order-of-magnitude analyses of fully-nonlinear Maxwell equations [50][51][52][53][54][55], are unnecessary in such regimes.

Governing Equation.
Within the scalar approximation [19][20][21][22][23][24], we consider an electric field of the form (2) which is time harmonic with angular frequency ω.The laboratory space and time coordinates are (x, z) and t, respectively.In medium j (where j = 1 and 2), it is well known that the complex spatial field E(x, z) satisfies the Helmholtz equation where c is the vacuum speed of light.The refractive index distribution n j (E) on either side of the boundary is introduced through n 2 j (E) ≡ n 2 0 j + α j |E| q , where n 0 j is the linear index at frequency ω and α j is a nonlinearity coefficient.To facilitate comparison with our earlier work [17,18,25,26], we look for travelling-wave solutions to (3) of the form E(x, z) = E 0 u(x, z) exp(ik 1 z).Here, E 0 is a (real) scale factor determining electric-field units, u(x, z) is the dimensionless envelope, and exp(ik 1 z) biases the solution in the forward longitudinal direction (taken to be z), where k 1 ≡ n 01 ω/c is the (linear) propagation constant of the carrier wave in medium 1.It then follows that u satisfies the inhomogeneous equation where h(x, z) is a Heaviside function that is equal to zero (unity) in the domain of medium 1 (medium 2).Equation (4) may be normalized with respect to the parameters in medium 1, in which case the following governing equation may be derived without further approximation [17,18,56,57]: ( The dimensionless longitudinal and transverse coordinates are ζ = z/L D and ξ = 2 1/2 x/w 0 , respectively, where L D = k 1 w 2 0 /2 is the diffraction length of a reference (paraxial) Gaussian beam.The inverse beam width is quantified by O (1), where ε ≡ λ/w 0 , and the field amplitude scales with is supplemented by the mismatch parameters [17,18,25,26] which determine how the linear and nonlinear refractive properties of the system change as one traverses the boundary.Equation ( 5) allows one access to material scenarios where Δ < 0 (i.e., configurations with n 02 > n 01 ) [17,18,25,26].By contrast, the scalings deployed in classic paraxial theory [8,9] restrict those models to consideration of regimes with Δ > 0. It is also apparent that setting κ ≈ 0 in an attempt to recover the paraxial model is going to lead to complications when handling the linear mismatch term Δ/4κ.The physical and mathematical difficulties of interpreting the paraxial approximation as the single-parameter limit κ ≈ 0 have been discussed at length elsewhere [23,24]; it is particularly well illustrated by interface geometries.

Solitons as a Nonlinear
Basis.When investigating the refraction of nonlinear light beams at material boundaries, it is essential to have an appropriate set of basis functions with which to formulate the problem.Such a basis is provided by the underlying exact analytical Helmholtz solitons [56].In the following analysis, we align the interface along the z axis so that it is located at transverse position x = 0. Medium 1 (the domain of the incident beam, where h = 0) is taken to occupy the half-plane −∞ ≤ x < 0, while medium 2 (the domain of the refracted beam, where h = +1), occupies 0 ≤ x ≤ +∞.
In medium 1, the governing equation ( 5) becomes Sufficiently far from the interface, (7) admits exact analytical solitons of the form [56] where η 0 is the peak amplitude of the beam, a = q[η q 0 /(2 + q)] 1/2 determines the (inverse) solution width, and quantifies nonlinear phase shift through the (typically small) quantity 4κβ 0 .The ± sign flags evolution in the forward/backward longitudinal direction.The propagation angle of the beam in the laboratory (i.e., the (x, z)) frame, denoted by θ inc and measured with respect to the z axis, is related to the transverse velocity parameter V inc through tanθ inc = (2κ) 1/2 V inc [23,24].In medium 2, u satisfies and one may derive similar families of solitons, Note that the connection between transverse velocity V ref and propagation angle θ ref , that is, unaffected by the (additional, linear) term Δ/4κ in (9) or by the nonlinear coefficient α.The geometry of these solitons, and their inherent stability against perturbations to the local beam shape, was explored in detail in [56].

Phase Continuity and Refraction.
In recent analyses, we have shown that arbitrary-angle refraction is well described by anticipating that the phase distribution of the light be continuous across the interface [17,18,25,26].Matching the phases of solutions (8a) and ( 10) at x = 0 leads to the requirement that Hence, continuity is possible if and only if the incident and refracted solitons share a common longitudinal sense (i.e., both must be in either the forward or backward directions).By rearranging (11), one can show that V ref is related to V inc through Expressed in this way, (12) provides a helpful form "V 2 ref = V 2 inc + deviation," where the sign of the deviation can be analysed separately.It is then instructive to define a net mismatch parameter δ as [17,18] This parameter can be interpreted as the sum of linear and nonlinear mismatches in material parameters.Its sign fully characterizes beam refraction.When δ > 0, one has that This regime is referred to as internal refraction, and it corresponds to the situation where the beam in medium 2 is deviated toward the interface (see Figure 1(a)).Conversely, δ < 0 implies that This external refraction regime corresponds to the beam in medium 2 being bent away from the interface (see Figure 1(b)).The special case of δ = 0 is the transparency condition, where linear and nonlinear index mismatches oppose each other exactly so that The interface is thus essentially transparent to the incident beam (see Figure 1(c)), which experiences no net change in dielectric properties as it crosses the boundary.It is interesting to note that the absence of an interface provides a parameter subset (i.e., Δ = 0 and α = 1) that satisfies the transparency condition identically.

The Helmholtz-Snell
Law for Spatial Solitons.By recognizing the rotational symmetry inherent to Helmholtz spatial solitons [23,24,56], it becomes clear that "forward" and "backward" designations are arbitrary.The only physical distinction between the two families is the propagation direction relative to the observer.By deploying trigonometric identities to eliminate velocities V inc and V ref , the forward and backward solutions in each medium may be written as and In this representation, the angles are bounded by −180 • < θ inc, ref ≤ +180 • with respect to the z-axis.By matching the solution phase at ξ = 0, one can obtain a compact Helmholtz-Snell refraction law involving laboratory-frame angles: where It is worthwhile noting that (15a) has a form which is almost exactly identical to that encountered when studying the classic electromagnetic problem of plane wave refraction at the boundary between different linear dielectrics.Thus, the single correction factor γ captures the interplay between finite-waist beams (through the appearance of κ) and discontinuities in both the linear and nonlinear properties of the adjoining media.The exponent q appears implicitly through β 0 .
When a beam encounters the boundary with a rarer medium, there is little penetration of light across that boundary until the incidence angle exceeds a critical value, denoted by θ crit .At criticality, where θ inc = θ crit , the trajectory of the incident beam is deviated so that, in principle, the outgoing beam travels along the interface (i.e., θ ref = 0).Applying this condition to law (15a) and (15b) leads to an analytical prediction for θ crit in terms of the mismatch parameters Δ and α and also the solution parameter 4κβ 0 : In practice, one rarely finds the refracted soliton travelling along the interface boundary since other effects tend to appear (we will return to this point later).

Universal versus Specific Representations.
There is clearly a universal flavour about ( 12), ( 13), (15a), (15b) and ( 16).For instance, there is no explicit mention of the system nonlinearity so that refraction is fully described by the mismatch parameters Δ and α and the beam parameter 4κβ 0 .These equations are, in fact, more general than they first appear; for instance, laws of exactly the same structure govern the refraction of plane waves in power-law materials: a wave with real amplitude u 0 has β 0 ≡ u q 0 (it is noteworthy that the refraction analysis for plane waves does not capture the modulational instability of such solutions in the single power-law context [58]).
The power-law nature of the problem becomes apparent after one substitutes for β 0 from (8b).The γ factor (c.f.(15b)) then becomes while the relation for the critical angle (c.f. ( 16)) is given by and the net mismatch parameter (c.f. ( 13)

Simulations of Solitons at Power-Law Interfaces
The Helmholtz type of off-axis nonparaxiality demands that the inequalities κ O( 1) and 4κβ 0 O(1) are always met, which is equivalent to the simultaneous requirements of broad beams with moderate intensities, respectively [23,24,56].Here, attention is restricted to configurations where the mismatch parameters are relatively small, typically α = O(1) and |Δ| O (1).We now proceed with a threestage analysis.The simplest case to consider first is that of linear interfaces.We then move on to investigate nonlinear interfaces and conclude by noting the dependence of GH shifts [48,49] on the nonlinearity exponent q.Stable solitons of the homogeneous power-law Helmholtz model tend to exist in the continuous interval 0 < q < 4 [27,56].For definiteness, we consider here only three discrete values: q = 1 (sub-Kerr), 2 (Kerr), and 3 (super-Kerr).
The fact that smaller values of κ give rise to better theorynumerics agreement, despite the increased magnitude of the linear-interface perturbation term at Δ/4κ, invites comment.We suspect that one possible explanation may lie in the origin of the Helmholtz-Snell law, whereby one matches solution phase (but not amplitude) at the boundary: the matching condition thus takes no account of amplitude curvature.In the laboratory frame, broader beams (i.e., those characterized by smaller κ values) tend to have lower amplitude curvature, and the corresponding spatial solitons (which play the role of nonlinear basis functions) thus map much more consistently onto the inherent assumptions of the analytical approach.
Upon crossing the interface, the refracted soliton may suffer small oscillations (in its amplitude, width, and area) reminiscent of those reported in previous studies [56], and be accompanied by a radiation pattern.Computations [59] have verified the effective independence of the refraction angle θ ref with respect to the incident amplitude η 0 .Accordingly, the curves in Figure 2 are essentially insensitive to q; they are quantitatively very similar to those obtained for q = 2 [10] and (when θ inc is sufficiently above θ crit in internal-refraction regimes) for q = 3.Any interaction between a spatial soliton and an interface generally involves three distinct components: a reflected beam, a refracted beam (sometimes more than one), and radiation.The way in which the incident energy is distributed amongst these components depends on a complicated interplay between the interface and beam parameters, and also the incidence angle.At very small angles (e.g., θ inc < 1 • ), the interaction can be highly inelastic and nonadiabatic (especially in external refraction regimes).Crucially, the single refracted soliton (as predicted in Section 2) dominates as θ inc approaches even modest nonparaxial angles, with reflected and radiation components hardly excited at all.The Helmholtz-Snell law embodied by (15a) and ( 15b) is, of course, most valid in such regimes.

Solitons at Nonlinear Interfaces.
Nonlinear interface effects dominate beam refraction when 4κβ 0 |1 − α| |Δ|.Without loss of generality, we isolate such effects by setting Δ = 0 so that the net mismatch parameter is given by δ = 4κβ 0 (1 − α) = 8κη q 0 (1 − α)/(2 + q).Refraction thus becomes far more sensitive to κ in nonlinear regimes (compare this to linear regimes, where δ = Δ is independent of κ).Theoretical predictions are shown in Figure 3.While there is generally good agreement with numerics for both κ = 2.5 × 10 −3 and κ = 1.0 × 10 −4 when α ≈ 1.0, the fit becomes less reliable for α = 2.0 and α = 0.3.For such parameters, the nonlinear refractive index change across the boundary is no longer small: one cannot expect to find such a close match because of strong nonlinear effects (e.g., beam splitting and radiation phenomena).While the fit is clearly better for broader beams (κ = 1.0 × 10 −4 ), the Helmholtz-Snell predictions for narrower beams (κ = 2.5 × 10 −3 ) are still in good qualitative agreement.
Detailed attention is first paid to regimes with α > 1 (external refraction, since δ < 0), where the nonlinearity is stronger in the second medium.Since the width of the refracted soliton is proportional to α −1/2 , it follows that the beam must become narrower as it crosses the interface.In this type of material regime, the incident soliton always has sufficient energy flow to excite a self-trapped soliton-like state in medium 2.
Simulations have shown that nonlinear external refraction tends to induce stronger oscillations in the parameters (amplitude, width, and area = amplitude × width) of the outgoing beam than in the linear case.Such oscillations are not captured by the adiabatic analysis in Section 2 (which anticipates a stationary state), but one expects their appearance intuitively.Qualitatively different effects can appear at quasi-paraxial incidence angles as the exponent q is varied; an illustrative example is shown in Figure 4 for θ inc = 3 • when α = 2.0.A unit-amplitude soliton exhibits a pronounced splitting phenomenon when q = 1 (see Figure 4(b)), whereby the field distribution in the second medium is shared between a dominant externally refracted beam (as predicted by analysis) and a weaker internally refracted component (there is also a low-amplitude reflected  component in the form of radiation modes).Since the internally refracted beam carries away some of the momentum of the incident beam, it follows that the dominant refracted beam travels at a smaller angle than that predicted by (15a) and (15b).This type of splitting is not present for unit-amplitude solitons with q = 2 (see Figure 4(b)), though it may appear for incident solitons with higher peak intensities [60].In such cases, the properties of the daughter solitons may be quantified with recourse to inverse scattering techniques.Splitting is also absent at q = 3 (see Figure 4(c)), though one finds quite a complicated radiation ripple pattern in the second medium.
Refraction in nonparaxial regimes tends to be a much cleaner process, with little radiation generated by the beaminterface interaction in comparison with quasi-paraxial regimes.Even at modest angles (e.g., θ inc = 30 • ), where the interface perturbation is distributed over a relatively short interaction length, the quantitative characteristics of the outgoing beam depend crucially on the power-law exponent.Both the depth of modulation and (longitudinal spatial) frequency of the oscillations tend to increase with q, as shown in Figure 5(a).When q = 2, the oscillations tend to vanish in ζ; for q = 1 and 3, they survive in the long-term evolution (this is also true for the oscillations shown in Figure 4(a)).A more detailed comparison of how the q affects beam refraction is presented in Figures 5(b)-5(d).
For material combinations with α < 1 (internal refraction, since δ > 0), the nonlinearity is weaker in the second medium.In that case, one should expect a critical angle to exist (in accordance with (17b)).If the incident soliton survives the interaction with the interface, then the refracted beam may be expected to undergo self-reshaping oscillations in its parameters, with the overall trend being toward an increase in solution width.Simulations have confirmed this to be the case, with diffractive broadening generally accompanied by a reduction in peak amplitude (see Figure 6(a))-these oscillations are reminiscent of those uncovered previously for perturbed initial-value problems [56].
Computations have uncovered a range of q-dependent effects, an illustrative sample of which is shown in Figure 6 for beams with κ = 2.5 × 10 −3 , a nonparaxial incidence angle θ inc = 30 • , and a nonlinear mismatch of α = 0.5.The (longitudinal spatial) frequency of the reshaping oscillations tends to decrease with increasing q (c.f., the increase with q when α > 1).Also at higher q values (e.g., for q = 3), a threshold phenomenon can appear whereby the energyflow [56] of the incident soliton may not be great enough to excite a refracted beam (if the energy flows associated with solutions (8a) and (10) are denoted by W inc and W ref , respectively, then it can be shown that W ref ≈ W inc /α 1/2 ).This instability is shown in Figure 6(d): upon colliding with the interface, the beam breaks up into radiation (this scenario is also present at quasi-paraxial incidence angles above the critical angle θ crit ).angle must depend on q (a prediction supported by simple inspection of Figures 4,5,and 6).At this point, it also becomes instructive to consider the trajectory of refracted beams more carefully.Detailed numerical calculations reveal that at quasi-paraxial incidence angles, the beam in the second medium tends to follow a straightline path.Such a simple notion of refraction, founded upon intuition from plane wave theory, is illustrated in Figure 7(a) for a nonlinear interface with α = 2.0 and a beam with θ inc = 3 • and κ = 2.5 × 10 −3 .However, if the incidence angle is increased into the nonparaxial domain (e.g., θ inc = 30 • ), a qualitatively different picture emerges.Now, the straightline path ξ − V ref ζ = 0 predicted by solution (10) defines an average trajectory about which the refracted beam "snakes."Snaking is more apparent with sub-Kerr nonlinearities (i.e., where q < 2), and it increases for narrower beams (i.e., larger values of κ) at a fixed amplitude (see Figure 8(a), where η 0 = 1.0).Beams with larger amplitudes also exhibit snaking, but oscillations tend to be more rapid in the longitudinal direction (see Figure 8(b)).The snaking phenomenon is most pronounced in regimes with α > 1, where the nonlinearity is stronger in the second medium.There is also an intrinsic dependence on θ inc that can be seen in Figure 7.For small angles of incidence, the incoming soliton experiences an interface perturbation that is distributed over a relatively long distance.The refracting beam is able to accommodate the inhomogeneity in the system since changes in focusing properties of the host medium occur gradually in the longitudinal direction.For larger-incidence angles, the effective beam-interface interaction length may be much shorter.Solitons impinging on the boundary then exhibit a sharp (rather than a gradual) perturbation whose action is to induce sustained oscillations.

Goos-Hänchen Shifts at Power-Law Interfaces.
Recently, GH shifts [48] have been investigated within the context of Helmholtz spatial solitons at Kerr-type material interfaces [49].These shifts describe the translation in the trajectory of a reflected beam relative to its position as predicted by geometrical optics.Extensive numerical investigations considered the interplay between incidence angle θ inc , material mismatches (Δ, α), and the nonparaxial parameter κ.Radiation-induced trapping was found to play a key role in determining the magnitude of the shift.Also uncovered were giant external GH shifts (in regimes with δ > 0 but where the second medium has a weaker nonlinearity (i.e., α < 1)).While a similar investigation of GH shifts in the power-law context is certainly outside our current scope, a small selection of results will now be presented to illustrate how they depend upon the nonlinearity exponent q.We begin by considering linear interfaces and unitamplitude incident solitons with κ = 2.5 × 10 −3 .According to (16), interfaces with Δ = 0.0025 have a theoretical critical angle of θ crit ≈ 2.86 • (this value depends only very weakly on q). Figure 9(a) gives a representative set of results.Inspection shows that, for any θ inc , the magnitude of the shift is generally greater for systems with q = 1 than for q = 2 or q = 3.The true critical angle (which can only be found through full simulations) is also slightly greater than that predicted by theory (for q = 1 and q = 2, θ crit ≈ 3.016 • and θ crit ≈ 3.030 • ; both angles exceed their theoretical values of θ crit ≈ 2.857 • and θ crit ≈ 2.859 • , respectively).While the qualitative behaviour of systems with q = 1 and q = 2 is largely very similar, strong qualitative differences have been uncovered in the case of q = 3.As θ inc approaches the theoretical critical angle, the incident soliton often becomes unstable against the interface perturbation.Large amounts of radiation tend to be generated by the interaction (c.f. Figure 9(d)), so that there is essentially no reflected or refracted beam and a GH shift is thus not easily quantifiable (or even meaningful).However, when θ inc is sufficiently above θ crit , the refraction angle is still well described by theory.
GH shifts at nonlinear interfaces have also been analyzed; results are presented in Figure 10 for α = 0.7 and where system nonlinearity has been augmented by considering incident solitons with η 0 = 2.0.Regimes with Δ = −0.001and Δ = −0.0025are associated with linear external refraction, while (13) shows that δ > 0 (i.e., for these parameter sets, net refraction is internal so that a critical angle should still exist).One general trend to emerge is that the true critical angle is slightly less than the theoretical value (c.f.linear interfaces of Figure 9, where the true critical angle slightly exceeds theory).However, it is worth noting that the qualitative behaviour predicted by ( 16), namely that θ crit increases with q, is supported by numerics.Close to the (true) critical angle, simulations show that there is a strong divergence in the GH shift (which becomes highly sensitive to θ inc ).Two other general trends are that (i) GH shifts are larger (sometimes notably) for q = 1 than for q = 3; (ii) in nonlinear regimes, the GH shifts depend more strongly on q than for the case of linear interfaces (compare Figure 10 to Figure 9(a)).Figure 10(b) reveals new types of behaviour at powerlaw interfaces when q / = 2.In particular, for q = 3 one enters a regime wherein the GH shift no longer increases monotonically with θ inc ; instead, there is a marked decrease in the shift before the divergence at θ inc ≈ θ crit sets in.These results provide clear evidence that one can, quite reasonably, expect to find new qualitative phenomena in material regimes that deviate from the idealized Kerr-type response.

Helmholtz Nonlinear Surface Waves
Surface waves are well known in nonlinear photonics, being stationary localized light states that travel along the interface between different media.The transverse mode profiles are typically asymmetric due to the differences in dielectric properties defining the interface.We now derive the surface modes of model ( 5) using solitons (8a) and ( 10) as a nonlinear basis.These new solutions are most conveniently parameterized by β, which is related to the propagation constant in paraxial theory [27,56].

Exact Analytical Solutions.
To proceed, one seeks solutions to (5) where k ζ is the propagation constant and F(ξ − ξ j ) is the (real) envelope profile that is centred on ξ j .After substituting for u and defining κk 2  ζ − 1/4κ ≡ β, it can be shown that in medium 1

u(ξ, ζ)
while in medium 2 For a nonlinearity exponent q, the surface waves associated with any given interface are parameterized solely by β.
The displacements ξ 1 and ξ 2 , as yet undetermined, can be found by considering the auxiliary equations that arise from respecting continuity of u and its normal derivative (here Goos-Hänchen shift Figure 9: Demonstration of the GH shift for a unit-amplitude (η 0 = 1.0) spatial soliton at a linear interface with Δ = 0.0025 and when κ = 2.5 × 10 −3 .(a) Variation of the GH shift with changing nonlinearity exponent q (the q = 3 results (inset) closely follow those for q = 2 until radiation effects come into play more strongly).(b), (c), and (d) show the full numerical solution |u(ξ, ζ)| of (5) when q = 1, 2 and 3, respectively (note that, over longer propagation lengths, the solution in (d) breaks up into radiation).The incidence angle in (b), (c), and (d) is θ inc = 3.016 • , which exceeds the (almost q-independent) critical angle θ crit ≈ 2.86 • .∂u/∂ξ or, equivalently, dF/dξ) across the interface.These conditions lead to respectively.After some algebraic manipulation of (19a) and (19b), one finds that where the parameters δ and μ are given by δ

Solution Families and Wave Power.
For both forwardand backward-propagating surface waves, there exist two solution families.The origin of this duality lies in solving simultaneous equations (19a) and (19b), where one is eventually obliged to find the roots of quadratic equations.Figure 11 reveals that, for fixed (Δ, α, β), the profile depends strongly on the nonlinearity exponent q.That is, the peak amplitude, width, and area all decrease with increasing q.The difference between the two peak amplitudes and the distance of each solution peak from the interface also decrease with increasing nonlinearity exponent.
Since the surface wave profiles differ, it is plausible that the two families will not share the same stability properties.We begin an analysis of Helmholtz solutions (18a) and (18b) by considering the power P, where as a function of the free parameter β for different values of the nonlinearity exponent q.The energy-flow invariant W [56] is related to P through W(β) = ±(1 + 4κβ) 1/2 P(β), where the ± sign here corresponds to forward-or backwardpropagating envelopes (being distinct from the sign choice in (20a) and (20b)).A representative set of curves is shown in Figure 12, where it can be seen that P(β) comprises two branches.In regime 1 (where Δ > 0 and α > 1), the lower (upper) branch corresponds to the −(+) sign in (20a) and (20b).This situation is reversed for regime 2 (where Δ < 0, 0 < α < 1), in which the lower (upper) branch corresponds to the +(−) sign (see Figure 11).We note that for lowerbranch solutions, the peak of the surface wave always resides in whichever medium has the lower linear refractive index.
Global trends in the parameter dependence of the modes profiles can be readily identified and discussed in the context of the two solution branches.For instance, one might fix Δ, β, and κ and consider the effect of varying α.In regime 1, one finds that upon increasing α, the upperbranch solutions tend to retain their shape while the lowerbranch solutions experience a decrease in amplitude, width, and area.The separation between the pair of solutions also becomes greater, with each localized wave moving away from the interface.As α is increased in regime 2, the lower-branch solutions tend to retain their shape while the upper-branch solution exhibits decreases in amplitude, width, and area.Also, the separation between the solutions tends to decrease with increasing α (so that the solutions move toward the boundary).

Surface Wave Stability.
Except near the intersection point (where β ≈ β min ), both P(β) branches satisfy the classic Vakhitov-Kolokolov (VK) criterion for stability; namely, dP/dβ > 0 [61].Extensive simulations have revealed that lower-branch solutions always tend to remain self-trapped within the vicinity of the interface (so long as dP/dβ > 0) evolving with a stationary profile over arbitrarily long distances.
Upper-branch solutions tend to display a spontaneous instability in finite ζ.A set of typical results is shown in Figure 13 for regime 1 with Δ = 0.005 and α = 2.0, where the input wave is localized predominantly in medium 1 (compare with Figure 11(a)).The initial stages of evolution appear to be stationary, but instability sets in after a finite propagation length.The unstable solution deviates spontaneously into medium 2, crossing the boundary and shedding radiation in the process.The beam in medium 2 undergoes narrowing since α > 1.For fixed interface and solution parameters, the instability growth rate clearly increases with q.However, the angular deviation of the (reshaping) daughter beam relative to the interface is largely insensitive to q.
Qualitatively different effects appear in regime 2 with Δ = −0.005and α = 0.5; this time, the input wave is localized predominantly in medium 2 (compare with Figure 11(b)).After a finite propagation length, the surface wave bends smoothly away from the interface and is deflected deeper into medium 2. There is relatively little radiation shed in this process, and the localized wave suffers only a very small change to its shape (largely because the beam remains always on the same side of the interface, so does not encounter changes in refractive index).In common with regime 1, the instability growth rate increases with q. 4.5.Interactions between Solitons and Surface Waves.The stability of lower-branch surface waves is now investigated by considering their resilience against interactions with spatial solitons.Only a brief summary is presented here since the primary motivation is to uncover qualitatively new effects that depend upon the exponent q (detailed quantitative analyses are reserved for future works).For definiteness, we present simulation results for collisions between a unitamplitude (η 0 = 1.0) soliton and surface waves in regimes 1 (Δ = 0.005, α = 2.0) and 2 (Δ = −0.005,α = 0.5) with β = 2.0 and κ = 2.5 × 10 −3 .
Regime 1 is considered first for a quasi-paraxial incidence angle of θ inc = 3 • (see Figure 14).When q = 1, the two distinct beams persist after the interaction.The path of the outgoing soliton has been deflected relative to its ingoing trajectory.The surface wave, on the other hand, survives as a localized spatial structure but can no longer be interpreted as a "surface wave" per se since it travels obliquely to (not along) the interface.This picture is qualitatively different for q = 2 and 3; there, the interaction results in the coalescence of the soliton and surface wave, producing a single higherintensity narrow filament travelling obliquely to the interface (narrowing is to be expected for medium combinations with α > 1).It is noteworthy that the propagation angle of the filament, relative to the interface, increases with q.Also, as one might expect, nonlinear beams interacting at quasiparaxial angles tend to shed a large amount of radiation.The qualitative behaviour can change dramatically at nonparaxial angles; a representative set of simulations for θ inc = 30 • is shown Figure 15.We have not observed coalescence phenomena; instead of this, individual beams retain their separate identities and can be clearly resolved.While the soliton often survives intact (and experiences a narrowing effect due to α > 1), the evolution of the surface wave depends strongly on the nonlinearity exponent: (i) for q = 1, it acquires slow modulations in its shape but remains localized within the vicinity of the interface (i.e., it remains essentially a surface wave after the interaction); (ii) for q = 2, its path is deviated by the interaction so that it no longer travels along the interface (this obliquelyevolving self-trapped structure is, by definition, not a surface wave); (iii) for q = 3, the collision destroys it completely.It is interesting to note the general trend that largerinteraction angles generate far less radiation than their paraxial counterparts [62].
We now turn our attention to similar interaction scenarios in regime 2. For a quasi-paraxial incidence angle of 3 • , the behaviour is strikingly different from that uncovered for the same angle in regime 1 (compare Figures 16 and  14, respectively).When q = 1, the soliton survives the interaction and the surface wave remains quasi-bound to the interface (but exhibiting a longitudinal "skimming" effect).5) when the nonlinearity exponent is q = 1 (surface wave follows interface), 2 (surface wave deflected), and 3 (surface wave destroyed), respectively.
For q = 2 and 3, the interaction deflects the surface wave away from the boundary (i.e., the surface wave becomes an obliquely-evolving beam).However, the behaviour of the soliton is different for q = 2 and 3: it survives intact in the former case and breaks up into radiation in the latter (this effect is related to the threshold phenomenon discussed in Section 3.2 and is not a consequence of the interaction with the surface wave).

Conclusion
We have presented, to the best of our knowledge, the first investigation of the way spatial solitons behave at the planar interface between dissimilar materials whose refractive index has a power-law dependence on the electric field amplitude.This analysis has thus extended arbitrary angle refraction considerations beyond the ubiquitous Kerrtype case [17,18,25,26].Exact analytical solitons have been deployed as a nonlinear basis [56], permitting the derivation of a generalized Helmholtz-Snell law.Extensive numerical computations have tested its predictions, which are most accurate in regimes where only the linear refractive index changes across the boundary.
A range of new quantitative and qualitative effects that depend strongly upon the exponent q has been identified.For example, simulations have found that, at linear interfaces with Δ > 0 and where q = 1 or 2, there is generally a well-defined transition (as θ inc increases) from soliton reflection, through GH shifting, to soliton refraction.In contrast, systems with q = 3 are often far more complex: the reflection-to-refraction transition is generally obscured by radiation effects over a finite band of incidence angles around the (theoretical) critical angle: solitons interacting with the interface may collapse into low-amplitude diffracting waves, with GH shifts becoming difficult to interpret or quantify in the absence of a well-defined reflected beam.However, strong supporting evidence has been obtained to confirm the validity of our Helmholtz-Snell modelling in arbitrary-angle non-Kerr regimes.In this way, the first steps have been taken towards understanding how (fully 2D) diffraction/nonlinearity interplays govern spatial soliton refraction in a much wider class of systems.
Nonlinear surface waves of model ( 5) have been derived, and we have performed the first numerical analysis of these types of solutions.Simulations have addressed the stability properties of the new surface waves, which tend to lie on one of two possible branches of the classic (β, P) curves.Solutions lying on the lower branch are predicted to behave as stable robust entities, while solutions on the upper branch are inherently unstable.Extensive computations have lent direct numerical support for this stability prediction in the more general Helmholtz context, and the growth rate of the upper-branch instability has been found to increase with q.
The stability properties of lower-branch Helmholtz surface waves have been further investigated by considering collisions with obliquely incident spatial solitons.A rich variety of behaviours, which depend crucially on both the nonlinearity exponent and the interaction angle, has been discovered.Finding analytical descriptions (e.g., through a perturbation theory [62]) of these phenomena seems a remote possibility since much of the behaviour is clearly nonadiabatic.Hence, computer simulations play a fundamental role in investigating solitons, surface waves, and their interactions in non-Kerr regimes.
The research presented in this paper provides a clear indication that deviating from the ideal Kerr-type nonlinearity (q = 2) can give rise to novel, interesting, and potentially exploitable phenomena.Each component of this paper (testing the Helmholtz-Snell law, calculating GH shifts, analyzing surface wave stability, and studying soliton-surface wave interactions) is a problem for detailed investigation in its own right.Our findings unpin analyses of other types of optical (and nonoptical) contexts involving solitons and surface waves where the power-law type of nonlinearity takes centre stage.One can expect other distinct classes of surface wave to exist when the interface comprises combinations of focusing/defocusing power-law nonlinearities [42,43,63]; the stability properties of these waves can, quite reasonably, be expected to differ from those reported here.Furthermore, the validity of our Helmholtz-Snell modelling in powerlaw regimes suggests that it may also be applicable to other material configurations, for example, to single-and multiinterface problems with cubic-quintic [64][65][66][67] and saturable [68][69][70] nonlinearities.Research is currently underway that investigates the generality of our findings in these other contexts, and preliminary results do suggest wider applicability.

Figure 1 :
Figure 1: Schematic diagram illustrating (a) internal (θ ref < θ inc ) and (b) external (θ ref > θ inc ) refraction in the laboratory frame.The transparency condition (θ ref = θ inc ) is shown in part (c).External refraction regimes tend to be highly angular and cannot be adequately described by the paraxial approximation.

Figure 7 (
b) quantifies this snaking effect for the external refraction simulations shown in Figures 5(b)-5(d).

Figure 13 :
Figure 13: Spontaneous instability of nonlinear surface waves lying on the upper solution branch of Figure 12(a), where κ = 2.5 × 10 −3 and β = 2.0 (interface mismatch parameters are Δ = 0.005 and α = 2.0).(a) Evolution in ζ of the peak amplitude |u| m of the beam.(b), (c), and (d) show the full numerical solution |u(ξ, ζ)| of (5) when the nonlinearity exponent is q = 1, 2, and 3, respectively.Note that the profiles of the input waves in (b) and (d) correspond to the upper-branch solutions shown in Figure 11(a).