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

H. M. Moya-Cessa; F. Soto-Eguibar; V. Arrizon; A. Zúñiga-Segundo

doi:10.1364/OE.24.010445

1. 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. In previous works [11, 12], we reported the revival and splitting of Gaussian beams, propagating in quadratic GRIN media. Such effects occur at propagation distances which are relatively large in comparison to the revival length associated to a GRIN medium, obtained elsewhere using either geometric optics or the first order expansion of the wave optics propagation operator. The indicated results 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 develop a significant generalization of such long periodic effects, assuming again a second order form of the propagator. We establish propagation distances that exhibit the revival and the splitting of any arbitrary input field, not as in our previous works [11,12], where the results were valid only for an initial Gaussian beam. 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). Both conclusions are derived analytically for any well behaved initial condition. As a simple example, we apply these results 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–20]. We also present, in Section 4, the case when the initial condition is an Airy function.

2. Helmholtz equation for GRIN media

The Helmholtz equation in two dimensions for a GRIN medium is

[\frac{\partial^{2}}{\partial X^{2}} + \frac{\partial^{2}}{\partial Z^{2}} + κ^{2} n^{2} (X)] E (X, Z) = 0,

where κ is the wave number and n(X) is the variable refraction index. For a quadratic medium, the refraction index can be written as

n^{2} (X) = n_{0}^{2} (1 - g^{2} X^{2}),

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

(\frac{\partial^{2}}{\partial z^{2}} + \frac{\partial^{2}}{\partial x^{2}} + k^{2} - x^{2}) E (x, z) = 0,

where we have introduced the dimensionless variables

x = \sqrt{η} X

and

z = \sqrt{η} Z

, and where we have defined η = n₀gκ and

k = n_{0} κ / \sqrt{η}

.

Introducing the operator $\hat{p} = - i \frac{d}{d x}$ , we can cast the last expression as [21]

\frac{\partial^{2} E (x, z)}{\partial z^{2}} = - (k^{2} - {\hat{p}}^{2} - x^{2}) E (x, z),

whose formal solution is

E (x, z) = exp [- i z \sqrt{k^{2} - {\hat{p}}^{2} - x^{2}}] E (x, 0) = exp [- i z k \sqrt{1 - \frac{{\hat{p}}^{2} + x^{2}}{k^{2}}}] E (x, 0),

where E(x, 0) is the boundary condition for z = 0.

We define the lowering ladder operator, â = (1/2)^1/2 (x + ip̂), its adjoint, the raising ladder operator, â^† = (1/2)^1/2 (x − ip̂) and the number operator n̂ = â^†â, and we write Eq. (5) as

E (x, z) = exp [- i z k \sqrt{1 - \frac{2 \hat{n} + 1}{k^{2}}}] E (x, 0),

2.1. The paraxial approximation

As a background for our main result in the next section, we present here an alternate derivation of the paraxial propagation in GRIN media, in terms of the FrFT [10, 22, 23]. The paraxial approximation is obtained when the square root in the exponential of Eq. (6) is expanded to first order, that is

E (x, z) = exp [- i z k (1 - \frac{1}{2 k^{2}})] exp (i \frac{z}{k} \hat{n}) E (x, 0) .

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 operator n̂, as 𝔉_α{f(x)} = exp(iαn̂){f(x)}, where alpha is the transform order. Considering this result, (7) can be rewritten as

E (x, z) = exp [- i z (k - \frac{1}{2 k})] 𝔉_{\frac{z}{k}} {E (x, 0)} .

Thus, the paraxial propagation to a distance z is proportional to the fractional Fourier transform of order

\frac{z}{π}

of the initial condition E(x, 0) [10, 22, 23].

2.2. 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

E (x, z) = exp (- i \frac{γ_{1}}{k^{3}} z) exp (i \frac{γ_{2}}{k^{3}} z \hat{n}) exp [i \frac{1}{2 k^{3}} z {\hat{n}}^{2}] E (x, 0),

where, for simplicity, we have defined

γ_{1} \equiv k^{4} - \frac{k^{2}}{2} - \frac{1}{8}, γ_{2} \equiv k^{2} + \frac{1}{2} .

We develop the initial condition in terms of the Gauss-Hermite functions

φ_{m} (x) = {(\frac{1}{π})}^{1 / 4} \frac{1}{\sqrt{2^{m} m!}} exp (- \frac{1}{2} x^{2}) H_{m} (x),

where H_m (x) are the Hermite polynomials, to obtain

E (x, z) = exp (- i \frac{γ_{1}}{k^{3}} z) exp (i \frac{γ_{2}}{k^{3}} z \hat{n}) \sum_{j = 0}^{\infty} c_{j} exp [i \frac{1}{2 k^{3}} z j^{2}] φ_{j} (x) .

Now, for z = lπk³ with l any non-negative integer, we have

E (x, z = l π k^{3}) = exp (- i l γ_{1} π) exp (i l γ_{2} π \hat{n}) [\sum_{j = 0}^{\infty} c_{2 j} φ_{2 j} (x) + i^{l} \sum_{j = 0}^{\infty} c_{2 j + 1} φ_{2 j + 1} (x)] .

Next, considering the identities

2 \sum_{j = 0}^{\infty} c_{2 j} φ_{2 j} (x) = \sum_{j = 0}^{\infty} c_{j} φ_{j} (x) + {(- 1)}^{j} \sum_{j = 0}^{\infty} c_{j} φ_{j} (x) = E (x, 0) + E (- x, 0)

and

2 \sum_{j = 0}^{\infty} c_{2 j + 1} φ_{2 j + 1} (x) = \sum_{j = 0}^{\infty} c_{j} φ_{j} (x) - {(- 1)}^{j} \sum_{j = 0}^{\infty} c_{j} φ_{j} (x) = E (x, 0) - E (- x, 0);

we obtain

E (x, z = l π k^{3}) = exp (- i l γ_{1} π) [\frac{1 + i^{l}}{2} exp (i l γ_{2} π \hat{n}) E (x, 0) + \frac{1 - i^{l}}{2} exp (i l γ_{2} π \hat{n}) E (- x, 0)] .

But, as we already said, 𝔉_α{f(x)} = exp(iαn̂){f(x)} [10], thus

E (x, z = l π k^{3}) = exp (- i l γ_{1} π) [𝔉_{l γ_{2} π} {\frac{1 + i^{l}}{2} E (x, 0) + \frac{1 - i^{l}}{2} E (- x, 0)}] .

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πγ₂ = 2mπ, where now m is another positive integer. As 𝔉_2πm is the identity operator, we get

E (x, z = l π k_{c}^{3}) = exp (- i l γ_{c} π) [\frac{1 + i^{l}}{2} E (x, 0) + \frac{1 - i^{l}}{2} E (- x, 0)],

where k_c = (2m/l − 1/2)^1/2 and γ_c = 3/8 + m (4m − 3l)/l². 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πγ₂ = (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,

z = π {(2 m - \frac{1}{2})}^{3 / 2}

,

k = \sqrt{2 m - \frac{1}{2}}

and

γ_{1} = m (4 m - 3) + \frac{3}{8}

. Below, we will study the Bessel functions and the Airy function as initial conditions, and this particular case will be exemplified.

In what follows, we present two simple examples of the behavior pointed above. However, we want to remark that our main results are analytical, our principal findings are Eqs. (17) and (18); the following examples are presented in order to clarify the ideas and to make them more concrete.

3. 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

\begin{array}{l} E (x, z = l π k^{3}) = & exp {i [1 - (γ_{1} + \frac{1}{2} γ_{2}) l π + \frac{1}{2} (x^{2} + \frac{b^{2}}{2}) tan (l π γ_{2})]} \sqrt{sec (l π γ_{2})} \times \\ {\frac{1 + i^{l}}{2} J_{ν}^{(2)} [x b sec (l π γ_{2}) + a, \frac{b^{2}}{4} tan (l π γ_{2}); - i] \\ + \frac{1 - i^{l}}{2} J_{ν}^{(2)} [- x b sec (l π γ_{2}) + a, \frac{b^{2}}{4} tan (l π γ_{2}); i]} . \end{array}

where

J_{ν}^{(2)} (ξ, ζ; i)

is the second order generalized Bessel function, defined in Eq. (32) of Appendix A.

In Figs. 1 and 2, we show the field intensity at the first splitting distance z = πk³, when the initial conditions are the Bessel functions J₀(x + 3) and J₃(x + 3), respectively. The parameters of the quadratic GRIN medium are n₀ = 1.5 and g = 10 mm⁻¹, in both cases. The dimensionless wave number is k = 1099.7, which corresponds to a wave number of κ = 8062.27 m⁻¹. In the left side of both figures, we show a 3D graphic of the propagation around z_c, which is z_c = 4.18 × 10⁹, and in the original variable Z_c = 3.80 × 10⁸ m. The right side of both figures show a cut of the 3D graphic at z = z_c, in order to display the agreement between the fractional Fourier transform, given by Eq. (19), the black line, and the numerical solution, the red dotted line; what we mean by numerical solution is explained in Appendix C.

Fig. 1 The square of the propagated electric field at the first critical distance when the initial condition is J₀ (x + 3). The quadratic GRIN media parameters are n₀ = 1.5 and g = 10 mm⁻¹. The beam has k = 1099.7, which corresponds to a wave number κ = 8062.27 m⁻¹. a). In the Z direction the 3D graphic goes from z_c − 200 to z_c + 200, where z_c = 4.18 × 10⁹, or Z_c = 3.80 × 10⁸ m. b). The black continuous line is the graphic of expression (19) and the red dotted line is the numerical solution Eq. (49).

Download Full Size | PDF

Fig. 2 The square of the propagated electric field at the first critical distance when the initial condition is J₃ (x + 3). The quadratic GRIN media parameters are n₀ = 1.5 and g = 10 mm⁻¹. The beam has k = 1099.7, which corresponds to a wave number κ = 8062.27 m⁻¹. a). In the Z direction the 3D graphic goes from z_c − 200 to z_c + 200, where z_c = 4.18 × 10⁹, or Z_c = 3.80 × 10⁸ m. b). The black continuous line is the graphic of expression (19) and the red dotted line is the numerical solution Eq. (49).

Download Full Size | PDF

For the special case, when the wavenumber is so that the fractional Fourier transform is just the identity, indicated in Eq. (18), we have

E (x, z = l π k_{c}^{3}) = exp (- i l γ_{c} π) [\frac{1 + i^{l}}{2} J_{ν} (b x + a) + \frac{1 - i^{l}}{2} J_{ν} (- b x + a)] .

In Figs. 3 and 4, we show the behavior of the propagating field around and in the critical Z distance, that in this circumstance is z_c = 6.36 × 10⁹ (Z_c = 5.03 × 10⁸ m). The 3D graphics, 3a and 4a, present the field propagation from z_c − 1000 to z_c + 1000. In the graphics, 3b and 4b, we plot (20) (black continuous line) and the numerical solution (red dotted line), Eq. (49). The GRIN media parameters are the same as in Figs. 1 and 2, n₀ = 1.5 and g = 10 mm⁻¹, and we took l = 1 and m = 8 × 10⁵, so k_c = 1264.91, that corresponds to a κ_c = 10666.6 m⁻¹. The result of the propagation it is just a linear combination of the initial condition with its specular image.

Fig. 3 The square of the propagated electric field at the first critical distance when the initial condition is J₀ (x + 3). The quadratic GRIN media has n₀ = 1.5 and g = 10mm⁻¹. We have taken l = 1, m = 8 × 10⁵ and thus k_c = 1264.91 (κ= 10666.6 m⁻¹). a). In the Z–direction the graphics goes from z_c − 1000 to z_c + 1000. b). The black continuous line is the graphic of expression (20) and the red dotted line is the numerical solution, Eq (49).

Download Full Size | PDF

Fig. 4 The square of the propagated electric field at the first critical distance when the initial condition is J₃ (x + 3). The quadratic GRIN media has n₀ = 1.5 and g = 10mm⁻¹. We have taken l = 1, m = 8 × 10⁵ and thus k_c = 1264.91 (κ = 10666.6 m⁻¹). a). In the Z–direction the graphics goes from z_c − 1000 to z_c + 1000. b). The black continuous line is the graphic of expression (20) and the red dotted line is the numerical solution Eq. (49).

Download Full Size | PDF

4. 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³ is given by

\begin{array}{l} E (x, z = l π k^{3}) = \sqrt{sec (l π γ_{2})} exp {i [- \frac{1}{2} l π γ_{2} + \frac{1}{2} (x^{2} - b^{2} a) tan (l π γ_{2}) + \frac{b^{6}}{12} {tan}^{3} (l π γ_{2})]} \\ {\frac{1 + i^{l}}{2} exp [- i \frac{b^{3}}{2} x tan (l π γ_{2}) sec (l π γ_{2})] Ai [b x sec (l π γ_{2}) + a - \frac{b^{4}}{4} {tan}^{2} (l π γ_{2})] \\ + \frac{1 - i^{l}}{2} exp [i \frac{b^{3}}{2} x tan (l π γ_{2}) sec (l π γ_{2})] Ai [- b x sec (l π γ_{2}) + a - \frac{b^{4}}{4} {tan}^{2} (l π γ_{2})]} . \end{array}

In Fig. 5, we show the field intensity at the first splitting distance z = πk³, when the initial condition is an Airy function Ai(x). The parameters of the quadratic GRIN medium are n₀ = 1.5 and g = 10mm⁻¹, and we have taken k = 1000.5 (κ = 6673.33 m⁻¹). The critical distance is z_c = 3.15 × 10⁹, or Z_c = 3.14 × 10⁸ m. The 3D graphic exhibit the comportment of the field around the critical distance z_c. In Fig. 5(b), the black continuous line is the fractional Fourier transform given in Eq. (21) and the red dotted line is the numerical solution, Eq. (49).

Fig. 5 The square of the electric field at the critical distance when the initial condition is Ai(x). The quadratic GRIN media parameters are n₀ = 1.5 and g = 10 mm⁻¹, and the wave number is k = 1000.5 (κ = 6673.33 m⁻¹). a). Behavior of the propagation around the critical distance z_c. The graphics shows from z_c − 200 to z_c + 200. b). The black continuous line is the graphic of expression (21) and the red dotted line is the numerical solution.

Download Full Size | PDF

For the special case indicated in Eq. (18), we have

E (x, z = l π k_{c}^{3}) = exp (- i l γ_{c} π) [\frac{1 + i^{l}}{2} Ai (x + 1) + \frac{1 - i^{l}}{2} Ai (- x + 1)] .

In Fig. 6, we show the behavior of the propagating field around and in the critical Z distance, that in this case is z_c = 2.25 × 10⁹, Z_c = 2.51 × 10⁸ m. The 3D graphic 6a, presents the field propagation from z_c − 200 to z_c + 200. In the graphic 6b, we plot (20) (black continuous line) and the numerical solution (red dotted line), Eq. (49). The parameters of the quadratic GRIN medium are n₀ = 1.5 and g = 10 mm⁻¹, and we took l = 1 and m = 8 × 10⁵, so k_c = 894.43, that corresponds to κ_c = 5333.37 m⁻¹. The result of the propagation it is just a linear combination of the initial condition with its specular image.

Fig. 6 The square of the propagated electric field at the critical distance when the initial condition is Ai(x). The quadratic GRIN media has n₀ = 1.5 and g = 10 mm⁻¹. We have taken l = 1, m = 4 × 10⁵ and then the wave number is k = 894.43 (κ = 5333.37 m⁻¹). a). Behavior of the propagation around the critical distance z_c. The graphics shows from z_c − 200 to z_c + 200. b). The black continuous line is the graphic of expression (22) and the red dotted line is the numerical solution.

Download Full Size | PDF

5. Conclusions

We have shown that for specific wave numbers, the propagation of any initial field distribution in a quadratic GRIN media, beyond the paraxial wave approximation, produces the revival and the splitting of the input field, at particular propagation distances. It is also proved that at some distances, the field is given by the fractional Fourier transform of the superposition of the initial field with a reflected version of it. These results are analytically and independent of the initial condition, as long as this last one is well behaved. We have applied these results, just as an example, 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 numerical solution and good agreement has been established.

Appendices

In Appendices A and B, we calculate the fractional Fourier transform of a Bessel function of the first kind of integer order and the fractional Fourier transform of an Airy function; we use the definition 𝔉_α{f(x)} = exp(iαn̂){f(x)} [10]. We want to point out, that there are other methods to obtain this fractional transforms [19, 20].

A. The fractional Fourier transform of a Bessel function of the first kind of integer order

It is known [23] that the fractional Fourier transform 𝔉_α{f} of an arbitrary function f is given by, 𝔉_α{f} = e^iαn̂{f}, where n̂ 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} (b x + a) = \frac{1}{2 π} \int_{- π}^{π} d τ exp {i [n τ - (b x + a) sin (τ)]} .

Then we have,

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} & = e^{i α \hat{n}} \frac{1}{2 π} \int_{- π}^{π} d τ exp {i [n τ - (b x + a) sin (τ)]} \\ = \frac{1}{2 π} \int_{- π}^{π} d τ exp (i n τ) exp (i α \hat{n}) exp [- i (b x + a) sin (τ)] . \end{array}

We now write the number operator as

\hat{n} = \frac{{\hat{p}}^{2}}{2} + \frac{x^{2}}{2} - \frac{1}{2}

and so

exp (i α \hat{n}) = exp (- i \frac{α}{2}) exp [i \frac{α}{2} ({\hat{p}}^{2} + x^{2})]

. But we know that [9]

exp [i ζ ({\hat{p}}^{2} + x^{2})] = exp [i f (ζ) x^{2}] exp [- i g (ζ) (x \hat{p} + \hat{p} x)] exp [i f (ζ) {\hat{p}}^{2}],

where

f (ζ) = \frac{1}{2} tan (2 ζ)

and

g (ζ) = \frac{1}{2} ln [cos (2 ζ)]

. Thus,

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} & = \frac{exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}]}{2 π} \int_{- π}^{π} d τ exp (i n τ) exp [- i g (\frac{α}{2}) (x \hat{p} + \hat{p} x)] \\ exp [i f (\frac{α}{2}) {\hat{p}}^{2}] exp [- i (b x + a) sin (τ)] . \end{array}

But,

exp [i f (\frac{α}{2}) {\hat{p}}^{2}] exp [- i (b x + a) sin (τ)] = exp [i f (\frac{α}{2}) b^{2} {sin}^{2} (τ)] exp [- i (b x + a) sin (τ)]

then, we have

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} = & \frac{exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}]}{2 π} \int_{- π}^{π} d τ exp (i n τ) exp [i f (\frac{α}{2}) b^{2} {sin}^{2} (τ)] \\ exp [- i g (\frac{α}{2}) (x \hat{p} + \hat{p} x)] exp [- i (b x + a) sin (τ)] . \end{array}

Also

\begin{array}{l} exp [- i g (\frac{α}{2}) (x \hat{p} + \hat{p} x)] exp [- i (b x + a) sin (τ)] = \\ = exp [- g (\frac{α}{2})] exp {- i (exp [- 2 g (\frac{α}{2})] b x + a) sin (τ)}, \end{array}

then the fractional Fourier transform of the Bessel functions of the first kind can be written as

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} = exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}] exp [- g (\frac{α}{2})] \\ \times \frac{1}{2 π} \int_{- π}^{π} d τ exp (i n τ) exp [i f (\frac{α}{2}) b^{2} {sin}^{2} (τ)] exp {- i (exp [- 2 g (\frac{α}{2})] b x + a) sin (τ)} . \end{array}

Writing

{sin}^{2} (τ) = \frac{1}{2} [1 - cos (2 τ)]

and changing the integration variable from τ to −τ in the last formula, we arrive to

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} = exp [- g (\frac{α}{2})] exp [- i \frac{α}{2} + i f (\frac{α}{2}) x^{2} + \frac{i}{2} b^{2} f (\frac{α}{2})] \\ \frac{1}{2 π} \int_{- π}^{π} d τ exp (- i n τ) exp [- \frac{i}{2} f (\frac{α}{2}) b^{2} cos (2 τ)] exp {i (exp [- 2 g (\frac{α}{2})] b x + a) sin (τ)} . \end{array}

Introducing now the Generalized Bessel Functions, defined as [13]

J_{n}^{(m)} (x, y; c) = \sum_{l = - \infty}^{\infty} c^{l} J_{n - m l} (x) J_{l} (y),

and its integral representation

J_{n}^{(m)} (x, y; e^{i θ}) = \frac{1}{2 π} \int_{- \infty}^{\infty} d φ exp [i x sin φ + i y sin (m φ + θ) - i n φ],

we can cast (31) as

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} = & exp [- g (\frac{α}{2}) - i \frac{α}{2} + i f (\frac{α}{2}) (x^{2} + \frac{b^{2}}{2})] \\ \times J_{n}^{(2)} [e^{- 2 g (\frac{α}{2})} b x + a, \frac{b^{2}}{2} f (\frac{α}{2}); - i] . \end{array}

Substituting the functions

f (ζ) = \frac{1}{2} tan (2 ζ)

and

g (ζ) = \frac{1}{2} ln [cos (2 ζ)]

, we finally obtain

\begin{array}{l} 𝔉_{α} {J_{n} (b x + a)} = & \sqrt{sec α} exp {\frac{i}{2} [- α + (x^{2} + \frac{b^{2}}{2}) tan α]} \\ \times J_{n}^{2} (b x sec α + a, \frac{b^{2}}{4} tan α; - i) . \end{array}

B. The fractional Fourier transform of an Airy function

We use again the known fact [23] that the fractional Fourier transform 𝔉_α{f} of an arbitrary function f is given by, 𝔉_α{f} = e^iαn̂{f} where n̂ 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 (b x + a) = \frac{1}{2 π} \int_{- \infty}^{\infty} d τ exp [i τ (b x + a)] exp (i \frac{τ^{3}}{3}) .

Then we have,

\begin{array}{l} 𝔉_{α} {Ai (b x + a)} = & e^{i α \hat{n}} \frac{1}{2 π} \int_{- \infty}^{\infty} d τ exp (i τ b x) exp (i τ a) exp (i \frac{τ^{3}}{3}) \\ = \frac{1}{2 π} \int_{- \infty}^{\infty} d τ exp (i \frac{τ^{3}}{3}) exp (i τ a) e^{i α \hat{n}} exp (i τ b x) . \end{array}

We now write the number operator as

\hat{n} = \frac{{\hat{p}}^{2}}{2} + \frac{x^{2}}{2} - \frac{1}{2}

, and so

exp (i α \hat{n}) = exp (- i \frac{α}{2}) exp [i \frac{α}{2} ({\hat{p}}^{2} + x^{2})]

. But we know that [9]

exp [i ζ ({\hat{p}}^{2} + x^{2})] = exp [i f (ζ) x^{2}] exp [- i g (ζ) (x \hat{p} + \hat{p} x)] exp [i f (ζ) {\hat{p}}^{2}],

where

f (ζ) = \frac{1}{2} tan (2 ζ)

and

g (ζ) = \frac{1}{2} ln [cos (2 ζ)]

. Thus,

\begin{array}{l} 𝔉_{α} {Ai (b x + a)} = & \frac{exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}]}{2 π} \int_{- \infty}^{\infty} d τ exp (i \frac{τ^{3}}{3}) exp [- i g (\frac{α}{2}) (x \hat{p} + p x)] \\ \times exp (i τ a) exp [i f (\frac{α}{2}) {\hat{p}}^{2}] exp (i τ b x) . \end{array}

But,

exp [i f (\frac{α}{2}) {\hat{p}}^{2}] exp (i τ b x) = exp [i f (\frac{α}{2}) b^{2} τ^{2}] exp (i τ b x),

then, we have

\begin{array}{l} 𝔉_{α} {Ai (b x + a)} = & \frac{exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}]}{2 π} \int_{- \infty}^{\infty} d τ exp (i \frac{τ^{3}}{3}) exp [i f (\frac{α}{2}) b^{2} τ^{2}] \\ \times exp (i τ a) exp [- i g (\frac{α}{2}) (x \hat{p} + \hat{p} x)] exp (i b x τ) . \end{array}

Also

exp [- i g (\frac{α}{2}) (x \hat{p} + \hat{p} x)] exp (i b x τ) = exp [- g (\frac{α}{2})] exp {i exp [- 2 g (\frac{α}{2})] b x τ},

then the fractional Fourier transform of the Airy function can be written as

\begin{array}{l} 𝔉_{α} {Ai (b x + a)} = exp (- i \frac{α}{2}) exp [i f (\frac{α}{2}) x^{2}] exp [- g (\frac{α}{2})] \\ \times \frac{1}{2 π} \int_{- \infty}^{\infty} d τ exp (i \frac{τ^{3}}{3}) exp (i τ a) exp [i f (\frac{α}{2}) b^{2} τ^{2}] exp {i exp [- 2 g (\frac{α}{2})] b x τ} . \end{array}

Completing a cube binomial and changing the integration variable, we get

\begin{matrix} 𝔉_{α} {Ai (b x + a)} = exp {- i \frac{α}{2} + i f (\frac{α}{2}) x^{2} + i \frac{2}{3} b^{6} f^{3} (\frac{α}{2}) - i b^{3} f (\frac{α}{2}) exp [- 2 g (\frac{α}{2})] x \\ - g (\frac{α}{2}) - i b^{2} f (\frac{α}{2}) a} \\ \times \frac{1}{2 π} \int_{- \infty}^{\infty} d τ exp (i \frac{τ^{3}}{3}) exp (i τ a) exp (i τ {exp [- 2 g (\frac{α}{2})] b x - b^{4} f^{2} (\frac{α}{2})}) . \end{matrix}

Finally, recalling the integral representation of the Airy function (36) and the definitions of the functions f and g, we obtain

\begin{array}{l} 𝔉_{α} {Ai (b x + a)} = & \sqrt{sec α} exp {i [- \frac{α}{2} + \frac{1}{2} tan α (x^{2} - b^{2} a - x b^{3} sec α + \frac{b^{6}}{6} {tan}^{2} α)]} \\ \times Ai (x b sec α + a - \frac{b^{4}}{4} {tan}^{2} α) . \end{array}

C. The numerical solution

In this appendix, we explain what we meant by numerical solution. Let us write the completely arbitrary initial condition as

E (x, 0) = \sum_{j = 0}^{\infty} c_{j} φ_{j} (x),

where φ_j(x) are the Gauss-Hermite functions (11), and

c_{j} = \int_{- \infty}^{\infty} E (x, 0) φ_{j} (x) d x .

Substituting in Eq. (6) and using the fact that n̂φ_j(x) = jφ_j(x), we obtain

E (x, z) = \sum_{j = 0}^{\infty} c_{j} exp [i z \sqrt{k^{2} - 1 - 2 j}] φ_{j} (x) .

What we call numerical solution is the finite version of the previous infinite sum (48). That is, the numerical solution is given by

E (x, z) = \sum_{j = 0}^{N_{\max}} c_{j} exp [i z \sqrt{k^{2} - 1 - 2 j}] φ_{j} (x) .

In the examples analyzed; i.e., an initial Bessel function and an initial Airy function, the coefficients (47) can be calculated analytically, and then evaluated numerically. We have used two criteria to determine the number of coefficients (47) and the number of terms in the sum (48); firstly, we have evaluated enough coefficients to reproduce precisely the initial conditions in a large x range; secondly, we have added terms to the sum until it does not change anymore. The number of terms needed, following the exposed criteria, was 300. All the graphics and the numerical calculations have been done in Mathematica 10.3.

References and links

1. C. Gomez-Reino, M. V. Perez, and C. Bao, Gradient-Index Optics (Springer-Verlag, 2002). [CrossRef]

2. E. Silvestre-Mora, P. Andres, and J. Ojeda-Castaneda, “Self-imaging in GRIN media,” Proc. SPIE 2730, Second Iberoamerican Meeting on Optics, 468 (February 5, 1996).

3. M. T. Flores-Arias, C. Bao, M. V. Prez, and C. Gmez-Reino, “Talbot effect in a tapered gradient-index medium for nonuniform and uniform illumination,” J. Opt. Soc. Am. A 16, 2439–2446 (1999). [CrossRef]

4. H. M. Ozaktas and D. Mendlovic, “Fractional Fourier transforms and their optical implementation. II,” J. Opt. Soc. Am. A 10, 2522–2531 (1993). [CrossRef]

5. J. Ojeda-Castañeda and P. Szwaykowski, “Novel modes in α-power GRIN,” Proc. SPIE 1500, 246 (1991). [CrossRef]

6. H. M. Moya-Cessa, M. Fernández Guasti, V. M. Arrizon, and S. Chávez-Cerda, “Optical realization of quantum-mechanical invariants,” Opt. Lett. 34, 1459–1461 (2009). [CrossRef] [PubMed]

7. S. M. Chumakov and K. B. Wolf, “Supersymmetry in Helmholtz optics,” Phys. Lett. A 193, 51–53 (1994). [CrossRef]

8. A. Zuñiga-Segundo, B. M. Rodríguez-Lara, D. J. Fernández, and H. M. Moya-Cessa, “Jacobi photonic lattices and their SUSY partners,” Opt. Express 22, 987–994 (2014). [CrossRef] [PubMed]

9. G. S. Agarwal and R. Simon, “A simple realization of fractional Fourier transform and relation to harmonic oscillator Green’s function,” Opt. Commun. 110, 23–26 (1994). [CrossRef]

10. H. Y. Fan and J. H. Chen, “On the core of the fractional Fourier transform and its role in composing complex fractional Fourier transformations and Fresnel transformations,” Front. Phys. 10, 1–6 (2015). [CrossRef]

11. F. Soto-Eguibar, V. Arrizon, A. Zuñiga-Segundo, and H. M. Moya-Cessa, “Optical realization of quantum Kerr medium dynamics,” Opt. Lett. 39, 6158–6161 (2014). [CrossRef] [PubMed]

12. V. Arrizon, F. Soto-Eguibar, A. Zuñiga-Segundo, and H. M. Moya-Cessa, “Revival and splitting of a Gaussian beam in gradient index media,” J. Opt. Soc. of Am. A 32, 1140–1145 (2015). [CrossRef]

13. G. Dattoli, L. Giannessi, L. Mezi, and A. Torre, “Theory of generalized Bessel functions,” Il Nuovo Cimento B 105(3), 327–348 (1990). [CrossRef]

14. G. Dattoli, A. Torre, S. Lorenzutta, and G. Maino, “A note on the theory of n-variable generalized Bessel functions,” Il Nuovo Cimento B 106(10), 1159–1166 (1991). [CrossRef]

15. G. Dattoli and A. Torre, “Note on discrete-like diffraction dynamics in free space: highlighting the variety of solving procedures,” J. Opt. Soc. of Am. B 31, 2214–2220 (2014). [CrossRef]

16. A. Perez-Leija, F. Soto-Eguibar, S. Chavez-Cerda, A. Szameit, H. M. Moya-Cessa, and D. N. Christodoulides, “Discrete-like diffraction dynamics in free space,” Opt. Express 21, 17951–17960 (2013). [CrossRef] [PubMed]

17. T. Eichelkraut, C. Vetter, A. Perez-Leija, H. M. Moya-Cessa, D. N. Christodoulides, and A. Szameit, “Coherent random walks in free space,” Optica 1, 268–271 (2014). [CrossRef]

18. A. Torre, “Propagating Airy wavelet-related patterns,” J. Opt. 17, 075604 (2015). [CrossRef]

19. M. A. Bandres and J. C. Gutíerrez-Vega, “Airy-Gauss beams and their transformation by paraxial optical systems,” Opt. Expr. 15, 16719–16728 (2007). [CrossRef]

20. A. Torre, “Gaussian modulated Ai- and Bi-based solutions of the 2D PWE: a comparison,” Appl. Phys. B99, 775–799 (2010). [CrossRef]

21. K. B. Wolf, Geometric Optics in Phase Space (Springer-Verlag, 2004).

22. V. Namias, “The Fractional order Fourier transform and its application to quantum mechanics,” J. Inst. Maths Applics. 25, 241–265 (1980). [CrossRef]

23. Pierre Pellat-Finet, Optique de Fourier. Théorie Métaxiale et Fractionnaire (Springer-Verlag, France2009). [CrossRef]

24. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables (Dover, 1972).

25. F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (NIST, 2010)

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

Abstract

1. Introduction

2. Helmholtz equation for GRIN media

2.1. The paraxial approximation

2.2. Beyond the paraxial approximation

3. A Bessel function as initial condition

4. An Airy function as initial condition

5. Conclusions

Appendices

A. The fractional Fourier transform of a Bessel function of the first kind of integer order

B. The fractional Fourier transform of an Airy function

C. The numerical solution

References and links

Cited By

Figures (6)

Equations (49)

Optics Express