Introduction

The research on accelerating beams has been growing rapidly since it was introduced into the domain of optics, by Demetri Christodoulides, in 2007 [1,2]. The best known accelerating beam is the paraxial Airy beam, which propagates along a parabolic trajectory while preserving its amplitude structure indefinitely. The effect is caused by interference: the waves emitted from all points on the Airy profile maintain a propagation-invariant structure, which shifts laterally along a parabola. This beautiful phenomenon has led to many intriguing ideas such as guiding particles along a curve [3], light-induced curved plasma channels [4], and to recent studies on shape-preserving accelerating beams in nonlinear optics [5–12]. The idea of nonlinear accelerating beams goes back to Giannini and Joseph [13] who found temporal self-accelerating solitons in Kerr media. However, that pioneering idea has been largely forgotten. Instead, scientists attempted to launch Airy beams into nonlinear media, and have found, in experiments and in simulations, that the beams do not propagate as an accelerating shape-preserving entity [9–11]. Except for cases where the nonlinear effects on the beam were weak [4], the beams deform and even breakup [9–11]. Thus, the prevailing opinion around 2010 considered strong nonlinearity as an obstacle for having accelerating beams. In spite of this general opinion, we have shown [5] that a specific design of the beam profile can make it accelerate in a shape-preserving manner while propagating in the nonlinear Kerr [5,6] nonlocal [7], and quadratic [5,8] media. Moreover, a recent study [12] has shown, by analyzing the spectrum, how an initial linear Airy beam launched into the nonlinear medium is evolving into our nonlinear solution [5,6]. In fact [12], proves that our nonlinear accelerating solution is an attractor for the nonlinear dynamics. All of these studies, however, were in the paraxial domain.

In fact, until recently, all studies on accelerating beams were on paraxial beams, namely, they have all dealt with solutions of the paraxial wave equation, where they are analogous to the accelerating wave-packets that are exact solutions of the Schrödinger equation [14]. Therefore, the propagation of such Airy beams had to stay paraxial, i.e., the trajectory was fundamentally limited to small (paraxial) angles. This restriction is a serious limitation, because the paraxial accelerating beam is moving on a curve which bends ever faster, and is eventually bound to break its own domain of existence. Indeed, using Airy beams beyond the paraxial domain leads to fast breakup – whenever the fine structure (oscillations) on the Airy waveform become on the scale of a few wavelength or smaller [15–17]. Interestingly, though, a different method to find accelerating wavepackets has done better beyond paraxiality: the caustics method. Using the caustics method one can design beams that accelerate along any convex trajectory [18]. A recent study [19] has used the caustics method to design beams that could bend to large non-paraxial angles, but such method inherently cannot really find shape-preserving solutions. Hence, the caustics-based non-paraxial accelerating beams in [19] did not display propagation invariance. Altogether, until recently it was believed that accelerating wavepackets are a unique phenomenon associated with Schrödinger-type equations, rather than being a general feature of waves. In spite of this prevailing opinion, several months ago we presented [20] non-paraxial shape-preserving accelerating beams that are solutions to the full Maxwell equations for time-harmonic fields. Finding these solutions for the Helmholtz equation has made them general and applicable for a very large class of wave phenomena in nature: ranging from acoustics and ultrasonics to waves on two dimensional membranes, and more. In the past few months, these beams were observed experimentally by both Dudley's group [21] and our group [22]. Among other things, we have found that non-paraxial accelerating beams can be designed to have features on the scale of a single wavelength. This stretches the spatial spectrum of the beam to the maximum possible, without crossing the line to the evanescent regime. But even so, the fundamental limit remains: the features of the beams cannot be smaller than a single wavelength without being composed of evanescent waves. Otherwise, the evanescent part would inevitably decay, causing the beam to change its shape dramatically. This raise an intriguing question: are we restricted by the same limitation in nonlinear optics? Is it possible to design a nonlinear accelerating beam with sub-wavelength features that would not decay? Another interesting question is related to the maximum possible bending angle: when the beam is launched from single planar aperture, bending up to 90° is the theoretical bound on linear beams [20–22]. What is the limit of bending of nonlinear self-accelerating beams? Can we design a self-trapped beam that would turn itself backward?

Here we present the first nonlinear self-trapped accelerating beams of the full time harmonic Maxwell equations [23]. These beams accelerate along a circular trajectory in a shape-preserving manner. We show that such nonlinear non-paraxial accelerating beams involve coupling between forward- and backward-propagating beams, giving rise to unique effects: Specifically, part of the beam is reflected due to the nonlinear change of the permittivity. However, the reflections are exactly compensated by the backward moving part of the wavepacket. We study the behavior of the beam in the absence of the backward part, showing that in the presence of a strong nonlinearity the trajectory of the shape-preserving beam is altered from being exactly circular acceleration. Moreover, at a very high nonlinearity, such nonlinear coupling between the forward and the backscattered waves leads to regime of instability, where unexpected highly nonlinear, possibly chaotic, dynamics arises. In addition, accelerating shape-preserving propagation occurs even when the beam is constructed of subwavelength features, such that a significant part of its power resides in the evanescent regime. In this context, this work proves that the nonlinearity takes a unique role in the propagation, preventing the evanescent spectrum from decaying.

Derivation of the solution

We begin from Maxwell's equations with nonlinearity of the form n²(E) = n₀²(1 + ε(|E|)) for a general saturable or Kerr-like nonlinearity that depends only on the intensity I = |E|². In general, the nonlinear permittivity ε is a tensor [24] that also depends on the electric field E, but we analyze here the case of a TE-polarized beam which does not mix the polarization. For a y-polarized time-harmonic field propagating in the x-z plane, Ē = E_Y(x,z,t)ŷ, Maxwell's equation yield the nonlinear Helmholtz equation:

This last issue is very important, namely, making sure there are physical parameters that support a nonlinear Maxwell accelerating beam. This is not so simple because, as we explained above, taking higher nonparaxial bending makes the nonlinear effects smaller in comparison, and vice versa. Hence, higher beam bending inevitably involves larger nonlinearities. To explain the tradeoff, consider an accelerating beam of high acceleration. As known from the paraxial [1,2] and the non-paraxial [20] regimes, this requires narrow lobes, which therefore experience faster diffraction broadening. This implies that in order to observe significant nonlinear effects – the nonlinear index change must be very strong. On the other hand, for a beam with wider lobes diffraction effects are weaker, hence weaker nonlinearities will affect the propagation dynamics. To give a qualitative explanation, we first find a relation between the width of the lobes to the bending of the beam, for the linear case: compare the acceleration of the linear non-paraxial beam, k/α ([20]), to the acceleration of an Airy beam, (2d³k²)⁻¹, where d is the scale of the Airy beam and is typically the width of the main lobe. From such comparison we get that α = 2(dk)³ = 16π³(d/λ)³~500 × (the number of wavelengths in the main lobe)³. To proceed, we obtain that needed number by considering the nonlinear effect on a single lobe. That is, from intuitive soliton considerations (the existence curve) [26], for a nonlinear index change of ~10⁻⁴, the typical width of a beam which would be significantly affected by the nonlinearity is around 20 wavelengths. Substituting this estimate forces us to take α to be at least 4 × 10⁶. From α, it is easy to find the needed transverse width of the beam: to achieve significant bending, the initial aperture should be at least half the radius of the circular trajectory ([20]), 0.5α/k, which is more than 4 meters, for the above parameters in the visible wavelength regime. The important consequence of this calculation is that one needs a strong nonlinearity to observe the nonlinear self-accelerating Maxwell beam. In particular, for the Kerr nonlinearity, the dimensionless nonlinear refractive index change is n₂I = 2(dk)⁻². Hence we can write α directly using the refractive index change: α = 2^(5/2)(n₂I)^(−3/2). A strong nonlinearity of n₂I~10⁻², which is nevertheless accessible in organic materials, yields α ~5,500, which would require an initial beam width of ~1cm to observe significant non-paraxial nonlinear bending.

Figures 1(a) and 1(b) show the beam profile and the corresponding plane-wave spectrum of a nonlinear accelerating beam, compared with those of the linear accelerating (Airy) beam. Both beams are truncated at the same finite aperture with its edges outside of Fig. 1(a). The aperture creates the ripples and the decaying tail found in Fig. 1(b). In all cases, the result is a beam which is shape-preserving along any circular curve (as is shown in Fig. 2(a) which is the solution of Eq. (2) transformed back to the x-z plane). One might be tempted to assume that this beam can propagate in circles endlessly, anti-clockwise, with its phase velocity α. However, the boundary conditions here are still at plane z = 0: the incident beam is launched as E(x,z = 0) and is propagating in the forward + z direction only. Therefore, the solution of (2) should be used only as the initial condition for x>0 in the z = 0 plane, where α>0 yields a beam bending anti-clockwise, or for x<0 in the z = 0 plane where α<0 generates a clockwise-bending beam. This selection of initial conditions breaks the full symmetry of the polar coordinates, and limits the acceleration of the beam to one quarter of the circle, after which the trajectory should (according to the solution of Eq. (2)) bend back toward -z. However, only a backward propagating beam can actually go back along the second quarter of the circle, and such a beam does not agree with the physical boundary conditions [20]. We add the black arrows in Fig. 2(a) to clarify this issue: the beam in Fig. 2(a) is composed of both forward and backward propagating parts, making it a physical solution only if one allows two boundaries, below and above, initiating two beams that propagate forward and backward respectively. Interestingly, nonlinearity creates an exception; since a backscattered wave coming from the nonlinear change in the refractive index can be back reflected. Hence, waves that propagate back exist even when boundary conditions are set such that only forward propagating waves are launched. Consequently, the exact solution of Eq. (2) is not the actual physical solution for a beam launched at z = 0, but instead it approximates the desired beam with increasing accuracy for larger α. This peculiarity is discussed in more details in [20,27], although an important difference should be noted: in linear media, the forward and backward propagating waves are completely decoupled, hence launching a forward propagating beam in homogeneous linear media will never cause any part of the beam to propagate back. On the other hand, adding nonlinearity to the same medium causes the forward and backward propagating waves to interact through the nonlinearity, hence they cannot be separated. In any case, to understand the actual physical dynamics we need to simulate the nonparaxial nonlinear dynamics of our solution. Such a simulation is challenging even with simple Gaussian-like beams, and has received much attention from past researchers [28–33]. Unfortunately, none of the previous techniques worked for our accelerating beams, which is why we had to develop two specific methods described below.

 

Fig. 1 Profile of the electric field (a) and amplitude of the Fourier transform (b) of a nonlinear (linear) accelerating beam in blue (red dashed). The edges of the evanescent regime are marked by a black dashed rectangle in (b), emphasizing in the enlarged window that the peak of the Fourier amplitude is at k_x>k₀. These beams have angular momentum of α = 60 for the nonlinear (blue) and α = 45 for the red (linear). The nonlinear change in the permittivity is ε_mean = 0.05, on the background of ε_max = 1.3; with λ = 1μm.

Download Full Size | PDF

 

Fig. 2 (a) The full self-accelerating self-trapped beam composed of both the forward and backward wavepackets, shows perfect shape-preserving bending of the nonlinear beam. (b) The self-accelerating nonlinear beam for physical boundary conditions at z = 0, showing an approximately shape-preserving bending. (c) The forward propagating part. (d) The backward propagating part. The arrows indicate the direction of propagation. Both panels have α = 160 and ε_max = 0.1. All scales are in microns; λ = 1μm.

Download Full Size | PDF

Decoupling the forward and backward waves

The coupling of the forward and backward waves is crucial for finding the correct physical solution. Had the system been linear, we could have used the projection in Fourier plane (x→k_x) which takes the form of $-\frac{i}{2}{\left(k_0^2-k_x^2\right)}^{-\frac{1}{2}}{\partial }_z+\frac{1}{2}$to exactly eliminates the backward propagating waves and gets the solution which satisfies the physical boundary conditions (i.e. beam is launched at z = 0). However, in a nonlinear wave system, the task of separating the wavepakcet into the forward and backward propagating parts is more complicated, but possible; see Fig. 2(c,d) for the forward and backward beams, calculated by the method explained in this section. We suggest a general technique for this separation, that can actually be utilized for any wave system which is second-order in the propagation parameter z. For each plane of constant z₀, we substitute the field into the nonlinear term to define a linear operator M = ∂_x² + k₀² + k₀²ε(I); Eq. (1) is now E_zz = -ME. Using the eigenvalues λ_i and the eigenfunctions u_i = u_i(x) of M, one can express any solution (close enough to the plane z = z₀) as a sum of the forward and backward propagating parts (first and second sums respectively)

Unfortunately, however, splitting the solution of Eq. (2) only at z = 0 into the forward and the backward parts is not sufficient for identifying the physical solution arising from a beam launched at plane z = 0 and propagating in the forward ( + z) direction only. This is because the problem is nonlinear, hence even after the division; one cannot eliminate the backward waves, since they are part of the physical solution. Mathematically, this means that the correct boundary condition (e.g., at z = 0) has an unknown back-propagating part. Physically, any incoming accelerating beam is expected to backscatter part of it, since it has a wide k-space spectrum. This creates a wave that propagates back which includes contributions from waves back-reflected from many planes, and might even reflect part of it forward again. In any case, the backward part of the beam at z = 0 cannot be set to zero. Rather, the backward propagating part at z = 0 depends on the incoming beam in a complicated manner. This issue will be revisited below, and as we show there, one can resolve it numerically via “reverse simulation”.

Coupled mode theory

Before proceeding to the full numerical solution, we suggest a simplified scenario. Assume that two beams are launched at opposite azimuthal directions. The first is launched counter-clockwise from θ = 0 and the second is launched clockwise from some other angle θ₁<90°. An exact solution of (2), such as Fig. 2(a), would be achieved only under three terms: First, if the forward propagating part of the beam (Fig. 2(c)) is chosen as the boundary at θ = 0. Second, the clockwise propagating part of the beam is chosen as the boundary condition at θ₁ (due to azimuthal symmetry we can use the values of the backward part of Fig. 2(d) on any θ). Third, the two beams must be mutually coherent. We should highlight the fascinating underlying physical insight in this: the backscattered beam of the counterclockwise part destructively interferes with part of the clockwise propagating beam, and vise versa. All in all, the two nonlinear wavepackets (Fig. 2(c,d)) add up into the exact solutions presented in Fig. 2(a). Consequently, the nonlinear coupling of the backward (or backscattered) and to the forward waves is vital in understanding of the full nonlinear dynamics. Yet, simulating the nonlinear dynamics is essential for checking the stability of the solution which is by no means guaranteed, and the numerical simulation is the strongest method to verify stability. However, simulating this coupled dynamics of the two beams is very complicated. Such simulation can be carried out in an iterative manner similar to that of [34], which studied the coupled nonlinear interaction between solitons propagating in opposite direction. Instead, we prefer to make use here of another approach, reverse simulation, which numerically simulates the propagation arising from a single beam launched at z = 0, with no pre-assumptions.

The reverse simulation

At an infinite distance z from the launch plane z = 0, we expect the beam to have only a forward propagating part, since no backward propagating beam would be back reflected from infinity, where the nonlinear index change is negligible so no reflections would occur. This means that at some large enough z₁ the desired solution should consist of only a forward propagating part. Our simulation utilize this fact: First, we take the solution of Eq. (2) (Fig. 2(a)) at z = z₁ and use the nonlinear projection of it onto the forward propagating wavepacket (Fig. 2(c)). Second, we propagate this beam in reverse until z = 0, by combining an explicit Runge-Kutta method with nonlinear beam propagation terms. In practice, of course we cannot take an infinite z value, but still our numerical simulation yields precise accelerating beams (Fig. 2(b) presents the beam calculated by our reverse simulation method). Note that the solution is close to the full (forward + backward) solution of Eq. (2) (Fig. 2(a)). The reason is that large α give a strongly anti-clockwise beam which has almost all the power in the forward part. Therefore, the initial beam can be chosen directly from the result of Fig. 2(a) with almost no change. The beam bends to an angle of 45°, while maintaining an almost shape-preserving profile. However, note that backscattering does play an important role – the trajectory is no longer an exact circle, and is instead stretched due to the absence of the balancing backward wave. The comparison to the exact circle is marked by the black dashed ark. For clarification, we emphasize that the beam of Fig. 2(b) is the physical solution of a beam launched at z = 0, while the beam of Fig. 2(a) requires two counter propagating beams, and is otherwise only a mathematical solution of Eq. (2).

The upper boundary of Fig. 2(c) is our upper boundary condition z₁ for the reverse simulation. Note that, since we removed the backward propagating beams at z = z₁ instead of removing it at infinite z, what we obtain is an approximate solution, which is actually the exact solution of the related problem of a finite-size nonlinear medium. This medium can be thought as being coated by some anti-reflecting layer which prevents reflections from the surface. Anyway, this “finite size medium approach” can be taken as large z-values as we need, and it is found to be a good (asymptotically exact) approximation for the true propagation dynamics inside an infinite nonlinear material.

Our method of nonlinear reverse propagation naturally assumes that the propagation dynamics of the nonlinear Helmholtz equation is reversible in z. Is it also correct in the nonlinear case? At first glance, it appears that the Eq. (2) is symmetric to the exchange z↔-z so we expect it to be reversible. However, the answer is actually that the simulation is reversible only if also the noise is unchanged. That is, in order for the nonlinear dynamics to be completely reversible, all the boundary conditions must be the same for the forward and the backward propagation. This means that one should use the exact field amplitude at the output plane (including the weak noise, amplitude and phase). In practice, this means that our numerical method of reverse simulation could lead to errors in the presence of very strong nonlinearities, which might causes the dynamics to be unstable. That is, under very strong nonlinearities small values of noise could play a major role, hence the reverse simulation might be extremely sensitive to numerical errors. However, in practice, we find that as long as the maximum nonlinear index change is smaller than ~1, the reverse simulation method works well.

Boomerang beam?

Figure 1(b) displays the Fourier transform of the beams of Fig. 1(a). The dashed line indicates the edge of the evanescent regime, highlighting that a considerable part of the spatial spectrum of the nonlinear beam is evanescent when the nonlinearity is absent. As such, this accelerating nonlinear beam heavily relies on the nonlinearity to propagate without decaying, otherwise an essential part of its spectrum would decay and the accelerating propagation would stop. This finding is manifested by the thin sub-wavelength lobes shown in Fig. 1(a), where even the main lobe is about half a wavelength wide. As a consequence, when the beam completes near 90° bending, it breaks up, causing the nonlinear effect to drop quickly, causing most of the beam to decay. For stronger nonlinearities, we expect an even more abrupt decay, causing a larger portion of the power to backscatter. This breaks a fundamental rule of linear optics: linear propagation prohibits the backward radiation to be created from an incident beam in a homogenous medium; the essential reason is simple: forward propagating waves cannot change to become backward propagating waves, so nothing that was launched forward can turn backwards. Indeed, linear wave dynamics do not couple the forward and backward waves. And yet, this rule breaks down in the presence of nonlinear effects, because nonlinearity couples all waves, in particular – in this case it couples the forward and the backward propagating waves. Nonlinear optics allows the nonlinear beam to reflect itself (or part of it) – which is expressed in an effective bend of more than 90°. This beautiful phenomenon is at its best in the nonlinear dynamics of our accelerating beams, since the back-reflection changes from zero - along the propagation, to its maximum - right after the nonlinearity drops. This dynamics suggests the exciting possibility of making a beam that would somehow fully reflect on itself. Such a notion is known in nonlinear optics in the context of phase conjugation; however, here the beams are accelerating (bending), hence the trajectory of such a beam would more resemble a boomerang. Such an envisioned beam will have a 'U' trajectory: it will start propagating along one branch of the 'U' and then bend to propagate back along the other branch of the 'U'. It seems that such a beam is indeed in principle physically possible. However, very high nonlinearities are required to cause enough power to backscatter so that an actual “boomerang soliton” would form. Further study of this phenomenon is left for future research.

Conclusions

We presented shape-preserving nonlinear self-accelerating beams of Maxwell's equations, and analyzed their physical features.

Acknowledgments

This work was funded by an Advanced Grant from the European Research Council, by the Israel Science Foundation, and by the Binational US-Israel Science Foundation.