Generalized revival and splitting of an arbitrary optical field in GRIN media

Assuming a non-paraxial propagation operator, we study the propagation of an electromagnetic field with an arbitrary initial condition in a quadratic GRIN medium. We show that at certain specific periodic distances, the propagated field is given by the fractional Fourier transform of a superposition of the initial field and of a reflected version of it. We also prove that for particular wavelengths, there is a revival and a splitting of the initial field. We apply this results, first to an initial field given by a Bessel function and show that it splits into two generalized Bessel functions, and second, to an Airy function. In both cases our results are compared with the exact numerical ones.


Introduction
Graded index (GRIN) media are mainly used in image formation applications [1]. On the other hand, it has been established that a GRIN medium can form self-images of periodic fields [2,3] and can support invariant propagation modes, either in the paraxial [4] and the non-paraxial domains [5]. In a different context, light propagation in a quadratic GRIN medium can be employed as a form of optical emulation of quantum phenomena. An example is the mimicking of quantum mechanical invariants by the propagation of light through the interface of two coupled GRIN devices [6]. Because the Schrödinger equation and the paraxial wave equation in classical optics are formally equivalent, cross applications between quantum mechanics and classical optics are common. One can extend the application in order to consider not only the paraxial regime but also the non-paraxial one, i.e. the complete Helmholtz equation. For instance, supersymmetric methods, common to quantum mechanics, have been proposed in classical optics [7,8].
An interesting conceptual and mathematical result is that the paraxial propagated field in a quadratic GRIN medium can be expressed as the Fractional Fourier transform (FrFT) of the incoming beam [4]. Indeed, the FrFT operator has been considered for discussing additional similitudes between classical optics and quantum formalism. For instance, Agarwal and Simon [9] have shown that Fresnel diffraction leads to the FrFT by noting that, when constructing the (quantum) harmonic oscillator evolution operator, it contains a term proportional to paraxial free propagation. In another report, by Fan and Chen [10], it is shown that the quantum-mechanical position-momentum mutual transformation operator is the core element for constructing the integration kernel of FrFT.
In the context of GRIN media, with quadratic refractive index dependence in the radial coordinate, knowledge from harmonic oscillator-type Hamiltonians can be used for the solution beyond the paraxial regime. As a significant result in this context, we established previously the long period revival and splitting of a Gaussian beam, transmitted by a quadratic GRIN medium [11,12]. Such effects are theoretically predicted expressing the Taylor series for the propagation operator up to the second order, i. e. including an additional term beyond the paraxial approximation.
In the present paper, we describe a significant generalization of such long period effects, assuming again a second order form of the propagator. We establish propagation distances that exhibit the revival and the splitting of an arbitrary input field. We also establish propagation distances for which any input field E(x) splits into two fields given by the FrFT of E(x) + E(−x). We apply this result to an initial field given by a Bessel function, and as a solution, we obtain the superposition of two so-called Generalized Bessel functions [13][14][15][16][17][18].

Helmholtz equation for GRIN media
The Helmholtz equation in two dimensions for a GRIN medium is where κ is the wave number and n(X) is the variable refraction index. For a quadratic medium, the refraction index can be written as where g is the gradient index in the X direction. So, for a quadratic dependence in the index of refraction, the Helmholtz equation is expressed as where we have introduced the dimensionless variables x = √ ηX and z = √ ηZ, and where we have defined η = n 0 gκ and k = n 0 κ/ √ η.

The paraxial approximation.
As a background for our main result in next section, we present here an alternate derivation of the paraxial propagation in GRIN media, in terms of the FrFT [10,21]. The paraxial approximation is obtained when the square root in the exponential of Eq. (6) is expanded to first order, that is On the other hand, it has been established that the fractional Fourier transform of a well behaved function f (x) can be obtained in terms of the number operatorn, as where alpha is the transform order. Considering this result, (7) can be rewritten as Thus, the paraxial propagation to a distance z is proportional to the fractional Fourier transform of order z π of the initial condition E(x, 0) [10, 21].

Beyond the paraxial approximation
We now allow ourselves to go one step further than the paraxial approximation. In Eq. (5), we again expand the square root in Taylor series, but we hold terms to second order instead, to obtain where, for simplicity, we have defined We develop the initial condition in terms of the Gauss-Hermite functions where H m (x) are the Hermite polynomials, to obtain Now, for z = lπk 3 with l any non-negative integer, we have Next, considering the identities and we obtain But, as we already said, Hence, at these periodic distances the field is the fractional Fourier transform of a superposition of the initial field and its specular image. It is clear that if the initial condition is symmetric, E(−x, 0) = E(x, 0), then at those periodic distances, we will have just the fractional Fourier transform of it. In particular, when l is congruent with 0 modulo 4, the field is the fractional Fourier transform of the initial condition times a phase factor. If l is congruent with 1 or with 3 modulo 4, we get the fractional Fourier transform of a superposition of the initial field and its specular image. In the case of l congruent with 2 modulo 4, we obtain a phase factor times the fractional Fourier transform of the specular image of the initial condition; of course, in this last case, if the initial condition is symmetric we will have just the fractional Fourier transform of the initial condition.
An interesting case is obtained when the dimensionless wave vector k is chosen such that lπγ 2 = 2mπ, where now m is another positive integer. As F 2πm is the identity operator, we get where k c = (2m/l − 1/2) 1/2 and γ c = 3/8 + m (4m − 3l) /l 2 . Thus, for those periodic distances and values of k, we can obtain the revival of the initial condition (when l is congruent with 0 mod 4), a superposition of the initial condition and its specular image (when l is congruent with 1 or with 3 mod 4) and the specular image of the initial condition (when l is congruent with 2 mod 4). A similar result occurs when lπγ 2 = (2m + 1) π, where again m is a positive integer, for which the fractional Fourier transform becomes the parity operator; but in this case E(x, 0), in Eq. (18) is replaced by E(−x, 0), and vice versa. The easiest situation is when we pick l = 1, and then, Below, we will study the Bessel functions and the Airy function as initial conditions, and this particular case will be exemplified.

A Bessel function as initial condition
In the particular case of a Bessel function as initial condition, E(x, 0) = J ν (bx + a), where ν is non-negative integer, we know from the Appendix A its fractional Fourier transform, Eq (35), thus where J ν (ξ, ζ; i) is the second order generalized Bessel function, defined in Eq. (32) of Appendix A. In Figure 1, we show the field intensity at the first splitting distance z = πk 3 , when the initial condition is a Bessel In

An Airy function as initial condition
We take now as initial condition the Airy function Ai(bx + a). Considering the fractional Fourier transform of this initial condition, Eq. (45) in Appendix B, the field propagated to the periodic distances z = lπk 3 is given by x tan (lπγ 2 ) sec (lπγ 2 ) Ai bx sec (lπγ 2 ) + a − b 4 4 tan 2 (lπγ 2 ) In Figure 3, we show the field intensity at the first splitting distance z = πk 3 , when the initial condition is an Airy In Fig. 4, we plot (22) (black continuous line) and the exact numerical solution (red dotted line). The GRIN media parameters are the same than in Fig. 3, and we took l = 1 and m = 4 × 10 5 , so k = 894.43. The result of the propagation is just a linear combination of the initial condition with its specular image.

Conclusions
We have shown that the propagation of an initial field distribution in a quadratic GRIN media, beyond the paraxial wave approximation, produces the revival and the splitting of the input field, at specific propagation distances. It is also proved that the field at another propagation distances is given by the fractional Fourier transform of the superposition of the initial field with a reflected version of it. We have applied the results to an initial Bessel function and found that the propagated field is given by a superposition of Generalized Bessel functions [13,14]. Also, the example when the initial field is an Airy function is examined. In both concrete cases, our predictions are checked against the exact numerical solution and good agreement has been established.

A The fractional Fourier transform of a Bessel function of the first kind of integer order
It is known [21] that the fractional Fourier transform F α {f } of an arbitrary function f is given by, wheren is the number operator of the quantum harmonic oscillator. To find the fractional Fourier transform of an integer order Bessel function of the first kind, we use the Jacobi-Angers integral representation [24,25] J n (bx + a) = 1 2π Then we have, We now write the number operator asn =p 2 2 + x 2 2 − 1 2 , and so exp (iαn) = exp −i α 2 exp i α 2 p 2 + x 2 . But we know that [9] exp iζ where f (ζ) = 1 2 tan(2ζ) and g(ζ) = 1 2 ln [cos(2ζ)]. Thus, then, we have Also then the fractional Fourier transform of the Bessel functions of the first kind can be written as Writing sin 2 (τ ) = 1 2 [1 − cos(2τ )] and changing the integration variable from τ to −τ in the last formula, we arrive to Introducing now the Generalized Bessel Functions, defined as [13] J (m) and its integral representation we can cast (31) as Substituting the functions f (ζ) = 1 2 tan(2ζ) and g(ζ) = 1 2 ln [cos(2ζ)], we finally obtain

B The fractional Fourier transform of an Airy function
We use again the known fact [21] that the fractional Fourier transform F α {f } of an arbitrary function f is given by, F α {f } = e iαn {f } wheren is the number operator of the quantum harmonic oscillator. To find the fractional Fourier transform of an Airy function, we use the integral representation [24,25] Ai(bx + a) = 1 2π Then we have, We now write the number operator asn =p 2 2 + x 2 2 − 1 2 , and so exp (iαn) = exp −i α 2 exp i α 2 p 2 + x 2 . But we know that [9] exp iζ p where f (ζ) = 1 2 tan(2ζ) and g(ζ) = 1 2 ln [cos(2ζ)]. Thus, then, we have then the fractional Fourier transform of the Airy function can be written as Completing a cube binomial and changing the integration variable, we get Finally, recalling the integral representation of the Airy function (36) and the definitions of the functions f and g, we obtain