Interactions of accessible solitons with interfaces in anisotropic media: the case of uniaxial nematic liquid crystals

We investigate, both theoretically and experimentally, spatial soliton interaction with dielectric interfaces in a strongly anisotropic medium with non-locality, such as nematic liquid crystals. We throw light on the role of refractive index gradients as well as optic axis variations in both voltage- and self-driven angular steering of non-local solitons. We specifically address and then employ in experiments a suitably designed electrode geometry in a liquid crystalline cell in order to define and tune a graded dielectric interface by exploiting the electro-optic response of the material through the in-plane reorientation of the optic axis in two distinct regions. We study both refraction and total internal reflection as well as voltage controlled steering of spatial solitons.


Introduction
Self-localization of light has been a core topic in nonlinear optics since its early days [1,2]. The first theoretical studies focused on the simplest nonlinearity, i.e. the local optical Kerr effect where the nonlinear perturbation of the refractive index is proportional to the intensity of light [3]. It was soon recognized that Kerr spatial optical solitons, i.e. nonlinear wavepackets conserving their shape and size in propagation owing to the balance of linear diffraction and self-focusing, are unstable in (2+1)D geometries [4]. Experiments carried out with pulsed laser beams led to the observation of a rich phenomenology, depending on the pulse duration and the specific features of the nonlinearity [5]. Light self-trapping with continuous wave lasers was initially observed in lead glasses via thermo-optic effects [6] and in sodium vapors via atomic transitions [7]. The stabilization of solitons in these experiments was ensured by saturation and non-locality in the nonlinear response of the material [4,8]. Despite the initial attention mainly devoted to bell-shaped bright spatial solitons, in the last two decades several studies have dealt with other families of spatial solitons (for a recent overview, see [9]).
In the limit of a high non-locality, spatial solitons are also called accessible as their propagation is governed by the same equation of the quantum harmonic oscillator [10]. Highly non-local responses characterize nematic liquid crystals (NLC) [11] and lead-doped glasses [12], their nonlinearities based on molecular reorientation and thermal effects, respectively. Spatial optical solitons in NLC are often referred to as nematicons [13,14]. Noteworthy, in both thermo-optic and reorientational media the nonlinearity is ruled by Poissonlike equations [15], with non-local range determined by the sample boundaries; the latter strongly affects soliton-soliton dynamics and allows long-range interactions [16,17], with potential applications to all-optical signal processing [14,[18][19][20]. Another type of nonlinearity is associated with the photorefractive response, where non-locality relates to the diffusion of light-induced carriers. In photorefractive media the nonlinear index well is not single-humped for a bell-shaped input and the transition from self-focusing to self-defocusing can be achieved via crystal rotation [21,22].
Investigations of spatial solitons have not been limited to homogeneous media: soliton interactions with linear inhomogeneities were theoretically addressed in nearly integrable systems [23], including a particle-like model for interactions with nonlinear interfaces [24,25]. Soliton emissions from waveguides featuring a nonlinear cladding [26] or from an amplifying potential trap [27] were numerically analyzed; more recently, efforts were devoted to the impact of the paraxial approximation in soliton reflection and transmission [28], as well as the interplay between nonlinear and Anderson localization in disordered media [29][30][31]. Experimentally, soliton scattering by local photonic potentials [32], surface solitons [33,34], soliton reflection and refraction at dielectric interfaces [35][36][37][38], soliton tunneling through a linear barrier [39] and solitons escaping a trapping potential [40,41] were demonstrated. In this context, through their unique properties, NLC allowed the observation of soliton equivalents of well-known optical phenomena, such as non-specular reflection [42] and the Goos-Hanchen shift [43]. To describe nematicon dynamics in the presence of photonic potentials, a modulational theory was recently developed [44][45][46].
An interesting property of a sub-class of spatial solitons is their ability to self-bend when their power changes. While the first proposed mechanism for the self-steering of solitary waves was based on phase modulation of their profile [47], experimental demonstrations were carried out in parametric media through cascading [48], in thermo-optic glass [49], photorefractives [50] and in NLC [51][52][53].
In this paper, we investigate in detail the interaction of the lowest-order nematicons (i.e. single-hump bell-shaped solitons) with planar dielectric interfaces. We address theoretically and numerically the role of walk-off in the interaction of a light beam with the interface, yielding non-specular reflection and power-dependent changes of the trajectory (i.e. selfturning). Experimentally we show how nematicons survive, as single entities, the interaction with an electro-optically defined interface when either transmitted (refracted) or reflected, encompassing in the latter case a markedly non-specular reflection with respect to the wave vector. We emphasize the tunability of the beam trajectory with the bias applied on each of the regions across the interface, demonstrating a larger angular span than in homogeneous NLC samples (where deflection is solely due to walk-off changes [54]) and an optimized deflection range using in-plane steering configurations.

Basic equations
We consider an NLC planar cell of thickness L along x and infinitely extended in the plane yz. The optical properties of a uniaxial NLC are completely established by the two scalar parameters n and n ⊥ , the refractive indices for electric fields parallel and normal to the optic axis (or molecular directorn), respectively. It is customary to introduce the optical anisotropy a = n 2 − n 2 ⊥ . We assume that the director, forming an angle θ withẑ, lies in the plane yz everywhere in the sample. In the absence of intense light, the NLC molecules are oriented with a one-dimensional profile versus y in the mid-plane θ (x = L/2, y, z) = θ 0 (y) (see figure 1(a)); Figure 1. (a) Top: sketch of the distribution of θ 0 versus y according to θ 0 = θ 1 + tanh(y/d) (solid black line) and piecewise linear (dashed red line) with a non-zero slope between −d/2 and d/2, with d the extent of the interface along y. Bottom: corresponding director distribution in the plane yz. Inset: definition of the director orientation θ. (b) Extraordinary refractive index versus θ (solid black line) and its approximation with a straight line according to n app e = n e (θ = 45 • ) + (dn e /dθ)| θ =45 • (θ − π/4) (red dashed line) for n || = 1.7 and n ⊥ = 1.5. Inset: relative error (n e − n app e )/n e . (c) Comparison between the exact formula n 2 e /n 2 b = n 2 e /n 2 b − 1 (black solid line) and the approximation n 2 e = 2 (n e /n b − 1) (dashed red line). Inset: corresponding relative error. Numerically computed profile of (d) θ and (e) the optical perturbation ψ versus y for θ in = θ 1 = 45 • and θ 2 = 65 • ; the NLC is the E7 mixture excited by a Gaussian beam in y = 0 (the interface center) with a waist w = 5 µm. The lines correspond to P 2D = 0.1 mW (blue), 1 mW (green), 5 mW (red) and 10 mW (cyan) from bottom to top, respectively.
in actual samples such a configuration can be obtained and controlled with interdigitated electrodes at the boundaries in yz [38]. We define ψ as the nonlinear contribution to θ due to the reorientational torque induced by extraordinary-polarized light, i.e. θ = θ 0 + ψ. We also recall the expressions for the walk-off angle δ = arctan { a sin (2θ ) [ a +2 ⊥ + a cos (2θ )] } and the extraordinary refractive index n e = ( cos 2 θ n 2 ⊥ + sin 2 θ n 2 ) −1/2 ; the former expressions of n e and δ are valid for plane waves with wave vectors parallel toẑ. In the above conditions and in the highly non-local limit, a simplified two-dimensional (2D) model can describe soliton propagation in the cell mid-plane x = L/2 [55]. Defining by δ b (z) the beam walk-off at a fixed z, in the paraxial approximation the nonlinear propagation of an extraordinarily-polarized beam, launched with a wave vector alongẑ, is governed by [56][57][58] where A is the magnetic field linearly polarized along x, , D y is the diffraction coefficient along y, k 0 is the vacuum wavenumber, n b is the extraordinary index of the carrier and n 2 e = n 2 e (θ ) − n 2 b is the refractive index landscape, including nonlinear (ψ) and linear effects (θ 0 ). To get a reasonable approximation, n b can vary with z according to n b (z) = n e (θ 0 ), i.e. following the linear changes in director distribution. In equation (2), we neglected the second-order partial derivatives of θ 0 along y and z; equation (1) is derived for negligible variations in walk-off δ across the beam profile and a wave vector only deflected in the plane yz. In agreement with a wave vector compelled to lie on the plane yz, hereafter we take θ 0 = θ in + ρ (y), with ρ a linear perturbation of the director orientation. Finally, by invoking paraxiality we ignore the variations of δ and n e with wave vector tilt in equation (1); hence, for large deflections we can expect discrepancies between the solutions of (1) and the exact fields computed from the complete Maxwell equations. Nevertheless, equation (1) permits fast and efficient numerical simulations while yielding consistent results with respect to exact solutions in most cases of interest.

Nonlinearity in homogeneous samples
The system of equations (1) and (2) can be easily treated in the homogeneous case θ 0 = θ in and in the perturbative limit (ψ θ 0 ) to get an effective non-local Kerr coefficient n 2H describing the dependence of the nonlinear response on the initial director orientation and on the physical parameters of the NLC. Setting n 2 e = n 2H P 2D g(y, z) we obtain [56] with g(y, z) = G 2D (y − y , z − z )|A| 2 dy dz /P 2D (G 2D is the Green function of equation (2)), δ 0 = δ (θ 0 ) and P 2D the effective beam power in the simplified 2D problem. For small anisotropies equation (3) provides n 2H = γ a sin [2(θ 0 − δ 0 )] sin (2θ 0 ).
For increasing powers the nonlinearity does not only affect the beam waist via selffocusing, but even the beam trajectory through nonlinear changes in walk-off [52,59]. Another scalar parameter, c δ = dδ b /dP 2D , can be introduced to quantify how self-deflection depends on the director rest angle θ 0 , providing [60]

Soliton path
The application of Ehrenfest's theorem to equation (1) provides a single ordinary differential equation which predicts the soliton trajectory in the plane yz via the function y s (z) = yp(y) dy where we defined the normalized wavefunction p(y) = |A| 2 / |A| 2 dy. This yields having defined the centered momenta µ m = p(y) (y − y s ) m dy. The right-hand side (rhs) of equation (5) is an equivalent force F ruling the evolution of the beam path. In particular, the first and second terms on the rhs of equation (5) depend on the gradient of the refractive index and the longitudinal (i.e. along z) variations in walk-off, respectively. For highly non-local solitons the term for m = 0 in equation (5) is much larger than the others, thus providing where θ b is the angle corresponding to the peak intensity and is given by θ b = θ in + ρ 0 + ψ 0 , with ρ 0 the value of ρ on the beam axis and ψ 0 the maximum of the nonlinear perturbation ψ; we also set F δ = d tan δ b /dz, with F δ being the contribution to F stemming from longitudinal variations in walk-off, both of linear (θ 0 ) and nonlinear (ψ) origin. We are interested in the trajectory dependence on the beam power, in both the perturbative (ψ θ 0 ) and highly nonlinear (ψ comparable with θ 0 ) regimes [61]. Qualitatively, we expect that the soliton can change its own path by means of two distinct mechanisms: changes in index gradient and nonlinear changes in walk-off related with the optical reorientation of the director.

Force induced by an index gradient
Hereby we discuss the role of the first term on the rhs of equation (5), focusing on the interplay between the linear distribution θ 0 and the nonlinear perturbation ψ. In the general case, the mth derivative of n 2 e with respect to y can be computed knowing the shape of θ along y and applying the chain rule; accordingly, for m > 0 the terms containing θ 0 and ψ are interconnected due to nonlinear terms in θ . The two contributions can be separated in the limit of small nonlinear effects, allowing the linearization of n 2 e (θ), as we will address in the following.

Perturbative regime.
As anticipated, a drastic simplification is allowed if small variations in n e are assumed (thus setting n 2 e = 2n b n e ; see figure 1(b)) and for a linear relationship between n e and the director orientation θ (thus setting n e = dn e /dθ (ρ + ψ) with dn e /dθ independent of y; see figure 1(c)). Under these conditions the force stemming from an index gradient consists of two contributions F ρ and F ψ , the former stemming from the linear director distribution θ 0 and the latter from the nonlinear perturbation ψ. Neglecting the contribution from walk-off in the rhs of equation (5), the soliton trajectory is given by where we defined the two force components as In the highly non-local limit, i.e. for very narrow beams with respect to both ψ and ρ, only the term for m = 0 survives in the rhs of equations (8) and (9). The overall force acting on the soliton is given by the index gradient computed at the beam peak, i.e. the soliton behaves like a particle and obeys geometrical optics [32,62]. We wish to stress that the term F ψ is non-zero due to the nonlinear relationship between the intensity profile |A| 2 and the optical reorientation ψ, as stated by equation (2). In fact, F ψ is zero in all nonlinear media where light-induced index perturbations are ruled by a linear equation (e.g. Poisson's equation in thermal media) and left-right symmetry is satisfied (not necessarily true in the presence of asymmetric boundary conditions [49]).
In order to study the behavior of F ρ and F ψ with power, we take a Gaussian ansatz for the beam intensity profile, i.e. |A| 2 = κ 0 (P 2D /w 2 ) exp [−2 (y − y s ) 2 /w 2 ], where κ 0 = 4n e cosδ b Z 0 /π (see [58,63] for details on field normalization). At this point, we neglect the waist dependence on power through self-focusing and consider a z-invariant profile, i.e. ∂ψ/∂z = 0. The latter condition corresponds to an infinite cell along z. Figures 1(d) and (e) present the overall director angle θ and the perturbation ψ versus the transverse coordinate y, computed by numerically solving equation (2). From figure 1(e) the nonlinear perturbation ψ conserves the even parity and the shape at each power provided the beam waist is much smaller than the barrier width d. In addition, the magnitude of the nonlinear perturbation saturates at high powers, i.e. when the director and electric field of the beam become parallel. It is apparent that, in the perturbative regime, F ρ depends only on the waist w and is independent of power P 2D , whereas the term F ψ depends on P 2D via the reorientation ψ, i.e. equation (2).
In order to investigate the behavior of F ρ with waist w, we consider a fixed distribution θ 0 (shown in figure 1(a)) and allow w to vary. Figure 2 plots F ρ versus w: for w d the force tends to an asymptotic value given by the index gradient at the beam peak, in agreement with equation (8); for larger w, higher-order terms (starting from m = 1) in the series have to be accounted for, thus providing a power dependence of the soliton trajectory given that the beam waist depends on the input power via self-focusing. Nonetheless, such dependence is rather weak for w/d 1 and thus can usually be neglected, particularly when dealing with nematicons. Finally, figure 2 shows the dependence of F ρ on the shape of the interface computed via equation (8): the calculated variations are small and do not significantly affect the propagation of the beam.
with the Green function G 2D = L/(2π) exp (−π |y − y | /L). In writing equation (10), we ignored boundary effects by assuming a cell infinitely extended along y. In the case of a homogeneous sample (i.e. ρ = 0 everywhere), F ρ is clearly zero and the nonlinear perturbation is even with respect to the beam axis in y = y s . Therefore, all the odd derivatives of ψ are zero; moreover, the beam retains its parity in propagation due to the symmetry of the system, making all the µ m zero for odd m. Summarizing, in homogeneous samples F ψ is zero and, therefore, a nematicon propagates along a straight line due to the absence of a mean force acting on it.
Let us now turn our attention to inhomogeneous samples. Equation (10) suggests that the magnitude of F ψ strongly depends on the ratio between the widths of |A| 2 , G 2D and ρ, determined by w, L and d, respectively. In the following, we assume d larger than w and L, consistently with the conditions used in deriving equation (2).
To study the behavior of F ψ , let us write ρ = ρ 0 + ρ(y) and set ρ ≈ ρ 1 (y − y s ), with ρ 0 = ρ (y s ) and ρ θ in + ρ 0 ; such a condition provides sin Equation (11) is valid in the nonlinear perturbative regime, when the relationship between the intensity |A| 2 and the perturbation ψ is linear and lends itself to the use of the Green's function formalism. Figure 3(a) shows the numerically computed ψ even versus y in the general case: for w L the perturbation ψ even follows the Green function G 2D except in the proximity of y = 0, where the derivative of ψ exists even if G 2D is not differentiable. In fact, for light beams much narrower than the Green function (i.e. the highly non-local limit), ψ even can be approximated as where ψ 0 is the peak of the even portion of the perturbation ψ which depends on the boundary conditions (i.e. the thickness L of the sample), corresponding to ψ(y = y s ) for a bell-shaped beam (graphed versus w at a given power in figure 3(b)). The net effect of ρ on ψ even is that of modulating the amplitude of the nonlinear index well, but not its shape. Let us now focus on ψ odd and its action on the beam trajectory. We take y s = 0 without loss of generality; hence, ρ 0 = ρ(0). We can find a closed form expression for ψ odd : In the highly non-local limit w → 0, equation (14) provides with u 0 the Heaviside step-function. Figure 3(c) plots ψ odd versus y, whereas figure 3(d) graphs the peak of ψ odd versus the beam waist w: the peak of the nonlinear perturbation ψ odd and its width go to zero for w → 0, i.e. when |A| 2 tends to a Dirac's delta function. Figure 3(d) displays the local force F local (proportional to ∂ψ odd /∂ y) for various w: as w decreases, the force distribution gets narrower, but at the same time the peak of F local becomes higher: in the limit w → 0 the local force F local tends to a Dirac's delta function. By substituting equation (14) into equation (9) and truncating the series at the first term (corresponding to m = 0), we obtain As expected, the nonlinear force is proportional to the optical power P 2D and the derivative ρ 1 of ρ with respect to y. The sign of the force depends on cos the nonlinear index well pushes the beam toward positive (negative) y for ρ 1 > 0 (ρ 1 < 0), whereas for π/4 < θ in + ρ 0 − δ b < π/2 the trend is opposite: F ψ is negative and deflects the beam toward negative (positive) y for ρ 1 > 0 (ρ 1 < 0). Noteworthy, for Figure 3(f) shows the exact behavior (i.e. directly from the general expression of Ehrenfest's theorem) F ψ versus w. In the limit w → 0, a non-zero force acts on the beam even in the absence of the odd part of the perturbation ψ odd (see equation (15) and figure 3(d)). We can also compute F ψ from the approximate solution (16) When using this approach, F ψ is non-zero in the limit w → 0, but taking a different value from the exact solution (see figure 3(f)) although equation (9) appears to predict the survival of just the term for m = 0 when w → 0. This discrepancy can be easily solved by noting that, for w → 0, the derivatives of ψ odd diverge due to the discontinuity of the first derivative of G(y) in y = 0; this implies that, in systems featuring a differentiable Green function, F ψ goes to zero as w → 0.
To elucidate the weight of F ψ in comparison with F ρ , we introduce the ratio r = F ψ /F ρ ; retaining only the terms for m = 0 in equations (8) and (9), in the highly non-local limit w → 0 we find A closed-form solution of equation (10) can also be derived in the local case when the beam profile is much wider than the Green function G 2D , i.e. w L; in this regime the nonlinear perturbation can be written as ψ(y) = γ L 2 π 2 sin [2 (θ in + ρ (y) − δ b )] |A(y)| 2 if ρ varies slowly across the Green function G 2D , as well. The force F ψ can be easily obtained from the expression of ψ using equation (9) truncated at m = 0, i.e. assuming a linear perturbation ρ slowly varying with respect to the intensity profile |A| 2 . The ratio r thus reads Noticeably, equation (18) tells us that r goes to zero for large w/L, in agreement with figure 3(f): in such a condition the nonlinear perturbation is negligible as the light intensity is too low to appreciably contribute to the overall distribution of the refractive index.
Concluding this section, in both local and highly non-local regimes the force F ψ stems from the y-derivative of the nonlinear coefficient n 2H given by equation (3), as apparent from the coefficient cos (17) and (18) (the changes of n e and δ with θ in were neglected). The non-locality, fixed by the ratio w/L, affects the size of F ψ with two asymptotic expressions, a constant value and (w/L) −2 , in the limits w/L → 0 (equation (17), highly non-local case) and w/L → ∞ (equation (18), local case), respectively.

Power dependence of the index gradient in the highly nonlinear regime.
Up to now we have addressed the force acting on the soliton in the perturbative regime, i.e. when ψ θ 0 . However, new relevant phenomena take place when the power excitation is high enough to appreciably modify the orientation of the NLC molecular director. In fact, while the component F ψ remains negligible in typical experimental conditions, the component F ρ , proportional to the linear index gradient (i.e. ∂ρ/∂ y), is affected by the director distribution along the beam path θ b (including optical contributions given by ψ 0 ), as stated by equation (6). The behavior of θ b for several input parameters w and P 2D is shown in figure 4(a). From a physical point of view, even if the gradient of ρ is fixed, the relationship between the director orientation θ and the corresponding extraordinary index n e (plotted in figure 1(b)) changes due to the optical perturbation in θ , leading to a power dependence in the gradient of n e (linked to ∂ρ/∂ y).
Let us start with a homogeneous sample of NLC; in the non-perturbative and highly nonlocal regime the optical perturbation ψ at first order can be approximated by [60] with ψ 0 (P 2D ) the peak of ψ for a beam power P 2D . Equation (19) is not valid in the proximity of the beam peak, in complete analogy with equations (12) and (13). Figure 1(e) displays the numerically computed profiles of ψ and confirms the validity of solution (19). Thereby, its shape is as in the perturbative case, but with a peak ψ 0 (P 2D ) obeying the equation dψ 0 /dP 2D = const × sin [2 (θ in + ψ 0 − δ b )] and undergoing saturation in agreement with the reorientation equation (2) [60] (see figure 4(a)). If a linear profile ρ much wider than L is added, it is clear that the shape, and thus the symmetry of ψ, will not change substantially in the non-perturbative regime when compared with the low-power case. Thus, F ψ remains negligible with respect to F ρ , even at high powers. It is convenient to define the nonlinear contribution to the force F nlin , computed by averaging the nonlinear local force in its complete form, i.e. F local ∝ ∂(n 2 e (θ) − n 2 e (θ 0 ))/∂ y, across the beam profile. Physically F nlin accounts for the variations in the force as induced by the perturbation ψ. In the low-power limit F nlin becomes the quantity F ψ defined by equation (9). We also introduce the linear force F lin , computed from the linear local force F local ∝ ∂n 2 e (θ 0 )/∂ y; such a contribution clearly does not depend on power (see section 2.4.1) and, for small index increments, coincides with F ρ defined by equation (8) and plotted in figure 2. Figure 4(b) shows the ratio r = F nlin /F lin versus beam waist w for three different powers. For P 2D = 0.1 mW the system is in the perturbative regime, with |r | increasing as w tends to zero in agreement with  (6). The powers are as (a, b), the changes with w being proportional to the beam power: specifically, the optical perturbation moves θ b away from the initial value of 52.5 • , thus decreasing dn e /dθ| θ b ( figure 1(b)). In these simulations, we took d = 100 µm, θ 1 = 45 • and θ 2 = 60 • .
figure 3(f), reaching a peak around 0.03 in agreement with equation (17). Noteworthy, as shown in figure 4(b), r is negative since θ b ≈ 52.5 • , thus making the cosine term in equation (17) negative. As the power increases, the nonlinear reorientation ψ is no longer negligible and r changes its behavior versus w. For any fixed power the modulus of the nonlinear force is maximum for small waists w owing to the larger reorientation, see figure 4(a). Finally, figure 4(c) compares the exact value of the overall force F and its approximation (given by equation (6)) versus w. The two sets of curves are nearly overlapping, demonstrating the validity of our theoretical predictions.

Highly non-local limit: particle-like nematicons and the role of anisotropy
Let us focus on the highly non-local limit, i.e. when w L, the case corresponding to most reported experiments (see section 4). In this limit solitons have a Gaussian profile and can be described by their position y s (z) and waist w(z) (proportional to √ µ 2 ). As stated by equation (5), the beam path is curved by two different forces stemming from the index gradient and the longitudinal variations in walk-off. Figure 4 shows that the nonlinear component F ψ can where we neglected terms proportional to µ m with m > 1. In equation (20) the propagation coordinate z and the refractive index profile −κ α n 2 e play the roles of an equivalent time and a potential, respectively, with κ α = D y /(2n 2 b ). Conversely, the second term on the rhs of equation (20) is an equivalent electric field, given by the temporal derivative of an equivalent vector potential proportional to tan δ b . We define the effective velocity v s = dy s /dz; geometrically, if α is the angle between the Poynting vector and z, we obtain v s = tan α. Moreover, after calling β the angle formed by the local wave vector and axis z, we introduce also v k = tan β. In the paraxial limit v s ≈ v k + tan δ. Using the identity D y = n 4 e (n 2 ⊥ n 2 ), multiplying equation (20) for v s and integrating between the initial (z = 0) and final states (labeled by z), we obtain where U = n 4 e /(4n 2 ⊥ n 2 ), and v s (0) = v 0 and tan δ(0) = tan δ 0 are the initial velocity (i.e. the initial angle between the Poynting vector and the z-axis) and the initial walk-off, respectively. Hence, we can define equivalent kinetic and potential energies as T eq = (1/2)v 2 s and V eq = −U − 1 2 tan 2 δ, respectively; thus, the equivalent soliton energy is E eq = T eq + V eq = (1/2)(v 2 s − tan 2 δ) − U . Figure 5(a) plots U and V eq versus θ . The integral term I = δ b δ 0 v k d tan δ in equation (21) represents the change in total equivalent energy due to the z-dependent equivalent electric field stemming from longitudinal walk-off variations.

Incident nematicon with a wave vector parallel to the interface.
We study the outcome of equation (21) by considering a self-confined beam impinging on the interface with the wave vector parallel to it. It is straightforward to derive v 0 = tan δ 0 : hence, owing to walk-off and at variance with propagation in isotropic media, the soliton interacts with the interface even if its phasefronts are normal to it. Equation (21) transforms into v s = ± [n 4 e (z) − n 4 e (0)] (2n 2 ⊥ n 2 ) + tan 2 δ + I . Since the leading terms of both v k and δ are proportional to a , the leading term of I is proportional to 2 a : thus, when v s is evaluated, in first approximation I can be neglected with respect to the differences in the refractive index. At this point, we can address the beam path when θ 1 is fixed and θ 2 is free to vary (see figure 5(b)). We start from θ 2 > θ 1 : the refractive index is larger in the second medium (i.e. past the interface); thus v k increases by the action of the index gradient. As θ 2 decreases, the barrier height n 2 e and the overall deflection reduce, with the beam undergoing a net motion toward the interface. When θ 2 = θ 1 (intersection between solid and black dashed lines in figure 5(b)) is satisfied, the sample is homogeneous and v s = tan δ 0 everywhere. If θ 2 decreases further, v k at the output section z becomes negative: despite v k < 0, the beam can still cross the barrier owing to walk-off. To evaluate the error introduced by neglecting the integral I in the equivalent energy conservation law, figure 5(b) graphs the angle α computed for a plane wave: a very good agreement is apparent between the two approaches, with slight differences arising for large deflection angle and in the proximity of θ 1 = θ 2 . The transition between refraction and total internal reflection (TIR) takes place for v s (θ 2 ) = 0, corresponding to θ 2 ≈ θ 1 − δ(θ 1 )/2. For further reductions in θ 2 , the beam still bounces back in the incidence region 1, conserving the same v s (z) (see the flat trend of α in figure 5(b)). Figure 5(c) compares the two approaches to the computation of the output angle α versus θ 1 . The solid line corresponds to −tan δ(θ 1 ) as predicted by the conservation of E eq , whereas the dashed line stems from a plane wave model. Slight differences stem from the n e dependence on β, a subtlety not accounted for in our energy-based model. In fact, by applying equation (21) and since the two quantities n e and δ do not change between initial and final states, we get v s (z) = −v 0 when TIR occurs, i.e. the Poynting vector undergoes specular reflection in our simplified model. We stress, however, that specularity in the Poynting vector implies strong non-specular reflection in the wave vector, with a change in modulus equal to twice the walk-off δ 1 .
Finally, let us name z p as the inversion point for the beam trajectory, that is, the z-coordinate corresponding to v s = 0. We want to demonstrate that the integral I in equation (21) is null in this case, i.e. the integrations from z = 0 to z p and from z p to z provide the same result in modulus, but with the opposite sign. The variations of v k are proportional to the gradient of n e and thus take the same value in both integrals, whereas the longitudinal variations in walk-off δ are equal in modulus, but with opposite signs, proving that our assumption I = 0 for TIR was correct. (21) can also be used when the input wave vector is at a finite angle with z. Figure 6 illustrates the results obtained for the output velocity and the comparison with a plane wave model. The overall steering angle is larger than for the case k ↑↑ẑ, mainly because of larger TIR angles. In fact, from (21) we easily get v s ≈ −v 0 following the interaction with the interface, with v 0 accounting for both the input tilt of k and the initial walk-off δ(θ 1 ). The larger deflections for tilted input pinpoint some inaccuracy in the particlelike model owing to paraxiality invoked in deriving equation (21). Additionally, with respect to the case k ↑↑ẑ, the error in imposing the conservation of E eq grows due to a larger I .  Finally, we stress that the nematicon trajectory can be computed by integrating v s , with sign determined by continuity in the ray trajectory. At the same time, it is important to note that our formalism is valid in the adiabatic approximation, i.e. in the absence of abrupt changes in the dielectric properties of the medium.

Numerical simulations of nematicons at planar interfaces
In this section, we will show and discuss numerical simulations of the (1+1)D system made by equations (1) and (2). We employed a code based on the beam propagation method (BPM) for the optical propagation and on an over-relaxed Gauss-Seidel algorithm for the director distribution [56]. The numerical window was chosen to be large enough across y to avoid spurious effects due to the interaction with the boundaries [51]. Noteworthy, including the derivatives of θ along z enables us to correctly model two important aspects of nematicon propagation: nonlinear changes in walk-off and longitudinal non-local effects [56]. Figure 7 displays the calculated evolution for an input wave vector parallel toẑ. The beam penetrates the barrier despite its phasefronts being normal toẑ, the latter behavior occurring thanks to the walk-off. The interface width d, taken equal to 100 µm as in experiments [38], is large enough to ensure adiabaticity, i.e. that the beam survives as a single entity after interacting with the linear potential ρ, when it is either reflected (negative jump in θ) or transmitted (positive jump in θ ).

Beam propagation method simulations
Let us first address the perturbative nonlinear regime. When θ 1 > θ 2 (figures 7(a)-(c)), the gradient of the refractive index is negative and the beam is repelled by the barrier. Therefore, after the interaction with an interface with a negative jump in the director orientation θ, the wave vector is always bent toward negative y. Importantly, this condition is not sufficient to ensure total internal reflection. In fact, even if k is tilted toward negative y, the Poynting vector can be tilted positively, the latter corresponding to a beam entering region 2. Finally, the optical penetration depth inside the barrier (i.e. the maximum y s ) depends on the shape of ρ for fixed θ 1 and θ 2 : the output position of the beam varies with d, even if the slope does not change in agreement with geometric optics.
When θ 1 < θ 2 (figures 7(e)-(g)) the force is positive; hence, k is tilted toward positive y and the beam undergoes gradual refraction, its pointwise trajectory depending on the distribution of θ .
The above results apply to both linear (figures 7(a) and (e)) and nonlinear cases (figures 7(b), (c), (f) and (g)). The main difference between the two regimes stems from diffractive spreading at low powers, as the beam undergoes strong deformations (or even splitting) when its size becomes comparable with the interface width d [62]. Conversely, when the power is high enough to ensure self-trapping, the beam does not suffer significant losses when interacting with the linear barrier. Moreover, under self-trapping the simulations demonstrate that its transverse profile conserves parity and is very close to a Gaussian.  D y (θ b ). At higher powers the beam self-focuses and eventually self-traps; as predicted, the soliton waist periodically oscillates along z, with the period getting smaller and smaller with excitation. Noteworthy, in the transmitted/refracted case a change in breathing period is observed due to different nonlinearities across the barrier, see equation (3). Figure 8 shows the profile of θ corresponding to figure 7. In the linear regime, θ ≈ θ 0 (figures 8(a) and (e), P 2D = 0.1 mW). For P 2D = 1 mW the nonlinear perturbation becomes appreciable (figures 8(b) and (f)) and promotes self-trapping. For P 2D = 5 mW the optical perturbation ψ becomes comparable with θ 0 (figures 8(c) and (f)); thus light propagates in the highly nonlinear regime. Finally, from figure 8 it is clear that the nonlinear index well retains an even parity even at large powers, consistently with section 2.4.2. Moreover, the transverse size of the nonlinear perturbation ψ remains equal to the cell thickness L, independently of the power P 2D and the discontinuity jump in θ . Figures 8(d) and (h) graph the angle perceived by the beam θ b versus z. When the beam is reflected, the optical reorientation decreases inside the linear barrier due to the lower nonlinearity (see equation (3)), whereas when light is transmitted such a difference is overcome by the variations of ρ. We note how θ b is a smooth function along z because we accounted for the non-locality in the direction longitudinal to the soliton propagation, whereas fast variations would have been observed if ∂ 2 θ/∂z 2 was neglected in equation (2).

Nematicon self-steering through linear interfaces
To complete the description of soliton propagation through linear interfaces, we need to address the changes in the nematicon trajectory with input power. Figures 9(a)   (TIR) than in refraction. The dependence of the nematicon path is consistent with nonlinear walk-off changes: when θ 1 = 30 • ( figure 9(b)) the walk-off increases with power, moving the Poynting vector toward positive y; when θ 1 = 60 • (figure 9(c)) the walk-off diminishes with power and the beam is deflected toward negative y; finally, if θ 1 = 45 • (figure 9(d)) the walk-off is close to its maximum versus θ: the soliton path does not vary appreciably for powers P 2D as large as 5 mW.
To estimate the relative roles of nonlinear walk-off and of its interplay with the force stemming from the index gradient on self-deflection, we simulated the beam evolution neglecting walk-off variations with power, i.e. setting δ b = δ 0 everywhere. Figures 9(e)-(g) plot the calculated trajectories for the same set of excitation parameters used in panels (a)-(d). To understand the interplay among these two forces we can examine figure 10, where the angles θ and θ correspond to the relative maxima of dn e /dθ 0 and δ, respectively. In the refracted case (θ 1 < θ 2 ) the index-based force F nlin is always positive with magnitude depending on dn e /dθ 0 . As illustrated in figure 10(b), for θ 0 lower than θ the force F nlin increases for increasing powers (dF nlin /dP 2D > 0), whereas for θ 0 > θ F nlin decreases (dF nlin /dP 2D < 0). Moreover, dF δ /dP is positive when θ 0 < θ and negative when θ 0 > θ , see figure 10(b). Summarizing, F nlin and F δ have similar trends versus θ when the soliton crosses the barrier. Nematicon refraction is illustrated in figures 9(a) and (e) and confirms our findings: the beam self-deflects toward positive y in both cases, the dominant contribution being F δ (in (e) self-deflection is not appreciable with the employed scale). When the soliton undergoes total internal reflection (θ 1 > θ 2 ) the force F nlin is negative and depends on dn e /dθ 0 , in analogy to the refraction case; hence, dF nlin /dP 2D changes sign with respect to the attracting barrier, a trend verified in figures 9(f) and (g). Conversely, when compared with the transmissive case, the force F δ retains the same sign and dependence on θ 0 . Figures 9(b) and (c) demonstrate that dF δ /dP 2D is dominant in the nematicon path, as the nematicon bends positively when θ < π/4 and negatively in the opposite case.
Finally, we also simulated light propagation when the equivalent potential acting on light is equal to the nonlinear perturbation ψ, i.e. setting n 2 e = const × ψ in equation (1). The results (not shown) indicate that self-steering is negligible: the output position of the beam varies by a few microns for the same set of parameters used in figure 9, mainly because of the longitudinal non-locality. Such small power-dependent displacements are consistent with the predictions from section 2.4.2.

Cell geometry and tunability with external bias
For the experimental study, we employed a sample as sketched in figure 11: a layer of the NLC mixture E7 is sandwiched between two glass slides held parallel to one another by Mylar spacers and separated by L = 100 µm. The inner surfaces of the slides are mechanically treated (rubbed) to provide director alignment at an angle θ R = 10 • with respect to the z-axis. Two pairs of comb-patterned ITO (indium tin oxide) electrodes with fingers parallel toẑ are deposited on the same surfaces, each pair controlling independently the director orientation in the region (volume) beneath, with the overall structure working effectively as a graded dielectric interface with a tunable barrier height. To clarify the operating mechanism of such a tunable interface, we can first consider a region far away from its center, conveniently located at y = 0. To a first approximation, an external voltage applied between the two combs forces the director reorientation in the plane yz, effectively producing the orientation of the optic axis controlled by the applied bias [56,64]. If two pairs of interdigitated electrodes are realized on each slide and aligned to each other, two regions with director orientations θ 1 and θ 2 can be formed outside a monotonically graded interface (of extension d) by biasing each pair at distinct voltages V 1 Figure 11. Sketch of the double comb-like configuration of the electrodes. The two combs generate a dominant y-component of the low-frequency electric fields E 1 and E 2 (corresponding to voltages V 1 and V 2 ), able to independently reorient the NLC molecules in regions 1 and 2, respectively. and V 2 , respectively (see figure 11). If the two comb electrodes on each slide are separated by less than the non-locality range (that is, the NLC thickness L) along y, then the interface has a quasi-linear distribution of θ (from θ 1 to θ 2 ) along y [57]. In our sample, we employed two comb electrodes symmetrically placed with respect to y = 0, with a third comb connected to the common ground potential and its center finger at y = 0 (figure 11). In order to ensure an electric disturbance much narrower than the Green's function of the reorientational response [56], each finger has a width of /4 = 15 µm and the finger-to-finger separation is /4 = 15 µm (figure 11); the resulting electric field has a periodicity while reorientation is periodic with /2 because of its dependence on the field squared [65]. A graded interface centered in y = 0 is therefore created with width d given approximately by the thickness L.
Finally, a third glass slide (not shown in figure 11) was placed orthogonally to the other two to realize the input interface, with rubbing along y to favor the excitation of an extraordinarily polarized field in the NLC layer; the input interface also avoids the insurgence of a meniscus which could, in turn, introduce undesired depolarization effects.

Experimental results
Light coupling in the cell was performed with an infrared beam (1064 nm) of waist w 0 ≈ 3 µm launched in the mid-plane x = L/2, avoiding interactions with the upper and lower boundaries [15,49]; its wave vector was orthogonal to the input interface, i.e. k//ẑ, and its polarization parallel to y to excite the extraordinary wave in the NLC layer. The input power was 5 mW to guarantee self-confinement throughout the whole range of parameters (orientation, incidence, etc). With the beam launched at y ≈ −250 µm, using a microscope and a CCD camera we monitored (through out-of-plane light scattering) its evolution in the sample versus the difference V = V 1 − V 2 between applied voltages and thus the slope of the linear barrier. For V 1 = 0 V, the effect of V 2 is to generate a linear barrier 'attracting' light, as n e (θ 1 ) < n e (θ 2 ). Entering region 2 (figure 11), the beam undergoes refraction and, after the interaction, its Poynting vector emerges at an angle α as large as 30 • for 0 V < V 2 < 3 V (see figure 12), depending on both the refractive index and the walk-off mismatches in regions 1 and 2, respectively (see section 2.5).
Conversely, when V 2 = 0 V and V 1 varies, the refractive index in region 1 is higher than in region 2; thus nematicons are repelled by the interface. It is worth highlighting that changes of V 1 introduce different walk-off angles δ 1 (θ 1 ), thus altering the incidence angle of the Poynting vector impinging on the interface. Starting from V 2 ≈ 0.9 V (for lower V 2 , despite the negative angle withẑ, the nematicon can cross the barrier because of walk-off, see section 2.5), n e (θ 1 ) becomes large enough (with respect to n e (θ 2 )) to make the nematicon bounce back in the halfplane y > 0 (figure 13), i.e. undergoing total internal reflection, with α reaching ≈ −15 • for A direct comparison of these data with the theoretical results in figures 5 and 6 shows a discrepancy in angles. In particular, comparing figure 13(f) with figures 6(b) and (d), the best fit of the measured angles with the theoretical results can be achieved for a tilt of 5 • in the input wave vector. Such apparent disagreement can be explained by considering that the comb electrodes do not start in z = 0 due to technological limitations: at the beginning of the sample, spurious components of the low-frequency electric field are present alongx andẑ and induce a modulation of the index along y able to tilt the wave vector. Consistently, nematicons propagate straight in region 1 when V 1 = 0 V (see figure 12(e)), at an angle equal to the calculated walk-off when the wave vector is at θ B = 10 • with the director; the latter (θ B = θ R ) is a confirmation that the input wave vector is along the z-axis.
Due to the high anisotropy and the large electro-optic response of NLC, these experimental results encompass the most important features of light reflection at the interface between anisotropic dielectrics, emphasizing the differences with the case of isotropic media. The quasispecular reflection observed for the beam implies a marked non-specularity for the wave vectors, as discussed in section 2.5. Moreover, owing to the dependence of the walk-off on V 1 and the quasi-specular reflection of the Poynting vector, anisotropic (uniaxial) media allow changes in TIR even without any tilt of the input wave vector.
Unfortunately, large voltages destabilize beam propagation due to director rotations out of the plane yz, occurring when the x-components of the low-frequency electric field(s) are no longer negligible. Nematicon stability with bias is also affected by the alignment of the comb electrodes (at 90 • with respect to previously used geometries [56]) despite the smallness of the periodic modulation of the director distribution (thus of the extraordinary refractive index) in the cell mid-plane. In fact, the periodic perturbation in our geometry is perpendicular to the input wave vector k; hence, the effects on light propagation are stronger than those associated with modulation along k [56]. The sample used in the experiments, in fact, corresponds to an array of shallow one-dimensional waveguides, their guiding properties getting stronger with the applied bias. For this reason, in these cells we did not observe nematicon self-steering [52] as the interaction with the periodic index perturbation (superposed to the index barrier) introduces energy leakage toward the bias-defined guides, leading to instabilities for excitations above 5 mW.
Similar dynamics, including nematicon refraction and TIR, is observed when both V 1 and V 2 are non-zero, with a repelling or attracting linear force F lin + F δ depending on the sign of V .  depending on both V 1 and V , in particular decreasing as V 1 increases. The condition V 1 > V 2 is not sufficient to ensure TIR: the output slope becomes negative only for V 2 = 1.3 V, whereas for V 2 = 1.5 and 1.8 V TIR is never achieved in the range of available biases V 1 . Noteworthy, for V 2 = 1.3 V the output slope is non-monotonic as δ 1 decreases for V 1 = 2 V, thus giving rise to a smaller angle of reflection, as discussed in section 2.5. Finally, the output position for V 2 = 1.5 and 1.8 V increases on going from the first to the second value of V 1 (see figure 14(a)) due to the dominant increase in δ 1 with respect to the reduced index gradient. Figure 14(d) plots the output slope versus V 2 . The biases are such that the nematicon always remains in region 1 (the output position is always negative, see figure 14(c)), even if the output slope is positive because the beam does not interact with the interface within the observation window. For the smallest V 1 = 1.2 V the output angle undergoes a continuous transition, whereas for V 1 = 1.5 and 2.3 V the transition is sharp, in agreement with figure 6. In summary, in the range of reported powers and voltages, nematicon robustness allows conserving self-confinement despite the interaction with a refractive potential barrier. Hence, by exploiting the light-guiding property of these spatial solitons, effective electrically controlled signal addressing/routing can be obtained with transverse displacements as large as 1.5 mm after propagation for 4 mm, i.e. an angular span of about 40 • (see figures 12 and 13).

Conclusions
Using a theoretical model, numerically tested and experimentally validated in several cases, we investigated in detail the interaction of spatial solitons in anisotropic media with a linearly graded dielectric interface. With specific reference to positively uniaxial and non-local NLC, we studied the dependence of the soliton trajectory on its own power by applying Ehrenfest's theorem. We showed that two main contributions produce nonlinear changes in the soliton path: the gradient in refractive index and the longitudinal variations of walk-off. For the first contribution, when the interface barrier is wider than the soliton waist, asymmetries in the nonlinear index well negligibly affect the beam path. In the same limit, power dependence in the force generated by a transverse index gradient is mainly due to perturbation of the linear index gradient via light-induced reorientation (the latter effect can be observed whenever (c) Partial sums of (A.1) versus number of terms considered. In all plots we took into account that the odd terms in (9) are null. a nonlinearity links the refractive index to the physical parameter modified by light). By full numerical simulations based on BPM, we studied nonlinear beam deflection through nonlinear walk-off. For fixed walk-off, we verified the dependence of index gradient on power, demonstrating that walk-off variations dominate the changes in the soliton path.
We employed Ehrenfest's theorem to understand the role of walk-off in a beam trajectory, preparing a simple particle-like model to describe the effects of anisotropy. We modeled known phenomena such as (wave vector) non-specular reflection at the interface between two uniaxial media, where Poynting vectors essentially undergo specular reflection (within a good approximation).
Experimentally, we prepared planar cells filled up with the liquid crystal E7, two couples of interdigitated comb-electrodes being deposited on each glass interface. The electrodes allow one to induce a (quasi-) linear interface across two NLC regions with optic axes oriented by external voltages. The measurements confirmed that nematicons are robust enough to survive the interaction with the graded interface, with exit beam angles in agreement with the theoretical predictions throughout the accessible ranges of voltages and excitations.
These results represent a step forward in the investigation of highly nonlinear effects in the interaction of spatial solitons with linear potential barriers; they show that NLC are an ideal platform for soliton physics, thanks to their large all-optical and electro-optic responses as well as their non-locality. We foresee potential applications in the realization of a new class of alloptical circuits encompassing both optical and electrical control.

Appendix. Series expansion for F ψ
The dependence of force F ψ on beam waist w can be cast in a closed form using the series (9). The derivatives of ψ odd to every order can be easily computed from (14) (their lengthy expressions are not reproduced here). Figure A.1 graphs the behavior of F ψ versus w for m up to 8 (corresponding to the 9th derivative of ψ odd ); the sign of the series terms alternates with m, with partial sums oscillating around the exact value. The convergence of the series can be more easily investigated in the highly non-local limit. Since µ m = 2 −m (m − 1)!!w m and ∂ m ψ odd /∂ y m | y s =C √ 2/π2 m+1 w 1−m (−1) k δ m,2k+1 (m − 2)!! (C = γ √ π κ 0 8 √ 2 D y dn e dθ | y=y s ρ 1 cos[2(θ in + ρ 0 − δ b )]P 2D ), in the limit w → 0 we find According to Leibniz's criterion, equation (A.1) is a convergent alternating series.