Open Access
Add to BibSonomyAbstract
One difficulty of angular spectrum representation method in studying optical propagation inside anisotropic crystals is to calculate the double integrals containing highly oscillating functions. In this paper, we introduce an accuracy and numerically cheap method based on asymptotic expansion theory to overcome this difficulty, which therefore allows to compute optical fields with arbitrary incident beam and is not restricted to the paraxial limit. This numerical method is benchmarked against the analytical solutions in uniaxial crystals and excellent agreements between them are obtained. As an application, we adopt it to investigate the propagation of a Gaussian vortex-beam in a biaxial crystal. The general features of anisotropic dynamics and power conversion between field components are revealed. The numerical results is interpreted by making appropriate analytical approximation to the wave equations.
© 2013 Optical Society of America
1. Introduction
Optical vortex is of significance in modern optics due to unique topological characters and its carrying angular momentum [1, 2]. It is a topological defect in the wave structure at which the field amplitude vanishes while the phase becomes uncertain. One important aspect of research on optical vortex is its propagating behaviors, which have already been revealed to exhibit many interesting features [1, 3–11]. For the free space or homogeneous isotropic medium, the propagation of light beam carrying optical vortices is ruled by the conservation law of total topological charge and angular momentum [3]. In uniaxial crystals, the displacement and splitting of optical vortex, as well as the coupling between spin and orbital angular momentum, are predicted and experimentally observed [4–12]. The studies on the propagating behaviors of vortex-beam in biaxial crystals are on the beginning which are based on simplified approximation to the wave equations [13, 14].
As far as the propagation of a general laser beam in uniaxially anisotropic media is concerned, various approaches have been proposed and developed in the past several decades (see for instance Refs. [15–22] and references therein). In particular, based on a plane-wave angular-spectrum representation, Ciattoni et al. develop a general scheme to describe the propagation of paraxial beam in uniaxial crystal [20]. The basic idea of this method is to express the optical field in uniaxial crystal as an superposition of an ordinary and extraordinary plane waves, whose propagation independently satisfy two decoupled parabolic equations. Two limits are discussed in details by Ciattoni et al.: the propagation parallel to and orthogonal to the optical axis [5, 20, 23], and many propagating properties are revealed [4, 24–26].
One crucial step of the angular spectrum representation method is to calculate double integrals containing highly oscillating functions [20, 27]. Nevertheless, the integration of highly oscillating function is generally difficult in practice, which also is an important ongoing issue in the community of applied Mathematics [28–32]. In the numerical respect, the traditional quadrature methods, which evaluate the integrand at a finite set of points and use a weighted sum of these values to approximate the integral, such as Newton-Cotes quadrature, Gauss-Legendre quadrature, etc., are inefficient for highly oscillatory integrals. The number of integration points (or function evaluations) in traditional quadrature methods grows as the oscillation of integrand increases, resulting in slow convergence speed and large absolute error [28,29]. Hence, the angular spectrum approach had thought to be not suit for numerical computations [27]. In [33], Lu et al. adopt the quadpack routines (Clenshaw-Curtis method) [34] to compute the oscillating integral and obtain excellent accuracy at small propagating distance, where the oscillation of integrand is not rapidly. However, the calculation becomes slower and also the accuracy decreases as the propagating distance increases. This situation motives us to look for better method to approximate the highly oscillatory integrals in the angular spectrum method, especially at large propagating distance.
On the other hand, previous studies on biaxial crystals mainly focus on the condition of conical diffraction (see, e.g., Refs. [15, 35–43]), or adopt assumptions such as slowly varying approximation [13, 44, 45]. In this paper, by extending the angular spectrum approach to the biaxial crystals [33], we introduce a numerical method based on asymptotic theory to accurately compute the optical field inside biaxial crystals with arbitrary incident beams along the crystallographic principal axes. The paper is organized as follows: In section 2, we review the general expressions for electric field inside biaxial crystals. In section 3, we present the numerical method and benchmark it against the analytical solutions in uniaxial crystals. In section 4, the propagation of a Gaussian vortex beam in biaxial crystals are investigated. Finally, we give conclusions in section 5.
2. The expression for electric field inside anisotropic crystals
The propagation of monochromatic light in an anisotropic crystal is described by the following equation [15, 20],
Here, E(r) is the complex amplitude of electric field E(r, t) with a relation E(r, t) = Re[E(r)exp(−iωt)], k0 = 2π/λ with λ being the vacuum wave length of light, and ε is the relative dielectric tensor which has form with nx, ny, and nz being the refractive index along crystalline x, y, and z axis, respectively. For the uniaxial crystals with nx = ny ≠ nz or nx ≠ ny = nz, the light propagating behaviors have been fully discussed [5, 20, 23, 46, 47]. What we will focus on in this paper is the general case of biaxial crystals with nx ≠ ny ≠ nz. We assume that the light is propagating along the z direction and the crystal fills the whole z > 0 space. In order to achieve the optical field inside the crystals, we first expand it using two-dimensional Fourier transformation [20, 23], where r⊥ = xêx + yêy and k⊥ = kxêx + kyêy denote the vectors in the transverse planes of real and momentum spaces. By substituting Eq. (3) into the wave equation of Eq. (1), it is straightforward to obtain the exact expression for Ẽ(k⊥, z), The coefficients λ1, λ2, c1, c2, c3, c4, c5, and c6 are given by [33] in which the 2 × 2 matrices P and Q are and the vector Ẽ⊥(0) is the two-dimensional Fourier transformation of the transverse incident field E⊥ (r⊥, 0) on the plane z = 0, In the above equations, the parameters α, β, γ, p, K, L, M are defined to beThe equations (3)–(7) consist of a complete description of optical propagation along z direction in an anisotropic crystal, whose forms are complex. For special input field such as paraxial Gaussian beam, one can simplify these equations to obtain analytical solutions for uniaxial crystal [5,20]. The analytical solution strongly depends on the form of incident beam and thus is complicated to find for biaxial crystals and non-paraxial limit. A numerical solution is needed for a general case. The main difficult of numerical method lies in computing the oscillating Fourier integral of Eq. (3), where the oscillation of integrand increases as the propagating distance z is increasing.
3. Numerical method
In this section, we present an effective scheme to calculate the highly oscillating Fourier integral in Eq. (3), based on the asymptotic expansion theory. For simplicity, we will present it in the framework of uniaxial crystal and benchmark it against the exact analytical solutions. For the uniaxial crystal with beam propagating along the optical axis, one can denote nx = ny = n0 and nz = ne. The expressions for the coefficients λ1 (λ2) and the matrices P (Q) in uniaxial crystals can be simplified as [20],
As a benchmark, We consider the propagation of a left-handed circularly polarized Gaussian beam in the uniaxial crystal. The incident field at z = 0 is where s is the waist of Gaussian beam and , which as well as its partner forms the basis of circular polarizations. Generally, the transverse part of any electric field E⊥ can be decomposed by using either the Cartesian basis or the circular basis, i.e., E⊥ = Exêx + Eyêy = E+ê+ + E−ê−. The relation between different components is The initial field in k space is a Fourier transformation of Eq. (17), which is given by The analytical solutions for the propagation of Gaussian beam of Eq. (17) in uniaxial crystal is derived by Ciattoni et al. in the paraxial limit [5, 24].3.1. Asymptotic approximation
The asymptotic expansion is a commonly used method to approximate the highly oscillatory double integral [15, 48],
with λ being a large positive parameter and D a bounded domain. In contrast to the conventional quadrature methods, the accuracy of asymptotic expansion improves as the oscillation (or λ) increases. The key point of asymptotic approximation is that the main contributions to the asymptotic expansion of I(λ) come only from regions in the vicinity of certain critical points [15, 48]. There are three types of critical point. The first type is the stationary points of f(x, y) within the domain D at which . The other two types of critical points are associated with the boundary of integral domain D. As far as the Eq. (3) is concerned, the boundary of integration goes to infinity and for an inputting Gaussian beam the kernel of integration vanishes on large boundary. So that only the first kind of critical points are important for the integral of Eq. (3) and we do not need to consider the other kinds of critical points.In the following, we present the leading term of asymptotic expansion to the integral of Eq. (3). As an example, we explicitly write out the expression for Ex,
where For large distance z, Ex consists of two typical highly oscillatory integral, which can be approximated by asymptotic expansion. It’s easy to show that both functions f1 and f2 have only one stationary point, which in the case of uniaxial crystals are at where ro = (x2 + y2 + z2)1/2, . For the dominated contributions to the integral in Ex come from the critical points of first kind, one can expand the functions f1 and f2 on the critical points (kxo, kyo) and (kxe, kye), respectively. For example, function f1 can be expanded as where the coefficients are , , and . In the same way, we also can obtain the partial derivatives for f2, , , and . After truncated the Taylor expansion of f1 and f2 to second order and substituting them into Eq. (21), we can derive the asymptotic approximation to Ex, Note that since the leading contributions to the integral Eq. (21) come from the stationary points we have approximated the functions c1(kx, ky) and c2(kx, ky) in the integral by constants c1(kxo, kyo) and c2(kxe, kye), which can be evaluated by substituting (kxo, kyo) and (kxe, kye) into Eq. (6). However, this approximation is valid only when the functions c1 and c2 are differentiable [48]. We will discuss this point in detail latter.The asymptotic expansion for the other two components, Ey and Ez, can be obtained in a similar way, by replacing c1, c2 with c3, c4 and c5, c6, respectively. After some straightforward algebra, the electric field E, in the case of uniaxial crystal, can be simplified as
This is one of the main results of present paper.After computing the values of function ci (i = 1, ⋯ ,6) at respective stationary points with a given incident field Ẽ⊥ (0), it’s straightforward to obtain the asymptotic expansion for E(r⊥, z) by using Eq. (27). In the following, we will qualify this approximation of asymptotic expansion by comparing it with the analytical solutions in uniaxial crystal [5,24]. The values of parameters in the calculation are as following: the refractive index no = 1.656, ne = 1.458, the wave length λ = 0.633μm, and the waist of initial Gaussian beam s = 4.59μm[33,49]. As mentioned above, the more oscillatory the integral, the more accuracy the asymptotic approxiamtion. Therefore we choose a large propagating distance: z = 100000μm. In Figs. 1(a) and 1(b), both asymptotic and analytical results for the moduli of E+ and E−, the circularly polarized components of propagating beam, are shown as a function of transverse radius r⊥. The results of asymptotic expansion agree very well with the analytical solutions, except for a small region close to the r⊥ = 0 point. For clarity, we show in Figs. 1(c) and 1(d) the absolute errors of asymptotic approximation with respect to the analytical solutions. One can see that the errors are smaller than 1 × 10−5 for E+ and most region of E− (e.g., r⊥ > 600μm). In the following, we will discuss why the errors of E− are large in the region close to r⊥ = 0.
Fig. 1 A comparison of asymptotic expansion and analytical solutions. (a) and (b) show the moduli of E+ and E−, the circularly polarized components of Gaussian beam propagating in uniaxial crystal, as a function of transverse radius r⊥. (c) and (d) show the absolute errors of asymptotic expansion for |E+| and |E−| with respect to the analytical results. The inset inside (d) is the errors of |E−| in the region of small radius r⊥. The values of parameters used in calculations: z = 100000μm, no = 1.656, ne = 1.458, λ = 0.633μm, and s = 4.59μm.
The asymptotic approximation of Eq. (27) is obtained by assuming that the functions ci(kx, ky) (i = 1, ⋯,6) are smooth and differentiable [48]. For circular components E±, one should look at the behavior of functions u± = c1 ± ic3 and v± = c2 ± ic4, respectively. It is interesting to observe that although the functions u+ and v+ corresponding to E+ are smooth everywhere, the functions u− and v− for E− are singular in the momentum space. As an example, we plot the real parts of functions u+ and u− in Fig. 2, from which we can see a singularity at the point (0, 0) for the function u−. For small values of r⊥ and large distance z, the stationary points (kxo, kyo) and (kxe, kye) are very close to the singular point (0, 0) in u− and v−, also see Eqs. (23) and (24). Therefore, the asymptotic approximation of Eq. (27), which does not consider the effect of singularity, is not good for E− in this case. But the effect of singularity is negligible when the radius r⊥ is large, for the stationary points are far from the singular point. This results in excellent accuracy for E− in large r⊥ region, as can be seen in Fig. 1(b).
Fig. 2 The real parts of functions u+ = c1 − ic3 (a) and u− = c1 + ic3 (b) in momentum space. The initial field is Gaussian beam of Eq. (17).
Therefore, we conclude that as long as the functions ci are smooth, the accuracy of asymptotic approximation is great for large propagating distance z. Note that the functions ci depend on the incident beams, which often contains singularities. In the following, we propose a method to take the effect of singularity into account.
3.2. Method of integrating near important points
As discussed above, the highly oscillatory integrals in the optical field are dominated by the contributions from the stationary points of functions f1 and f2, and also the singular points in functions ci. Thus, a straightforward idea of approximation is to numerically integrate over small regions near such important points and throw away the contributions from the other more oscillatory and smooth regions [29]. The part of the integral which is thrown away decays exponentially fast as the oscillation increases [29,48]. Therefore, we can deduce that, the larger the propagating distance z is, the smaller region is needed to be integrated over in order to obtain a fixed accuracy. The integrand close to the stationary points does not oscillate rapidly, therefore the integration over a nearby region can be performed by using the classic quadrature methods such as Clenshaw-Curtis and Gauss-Kronrod [33, 34], which makes this scheme easy to implement.
To verify this method, we evaluated the optical field by integrating over small areas centered at the stationary points (kxo, kyo) and (kxe, kye). How to choose the domain of integration is a little bit tricky. It is not a good idea to choose a circular domain around the stationary point, for the value of integral will oscillate as the radius of domain is increasing. It turns out that a square domain is a good choice (see the inset in Fig. 3(b) for the shape and position). The oscillating part of the integrand on the four corners of square can partially cancel each other, making the result of integration much more smooth.
Fig. 3 The absolute errors of |E+| (a) and |E−| (b) as a function of the width dk of the integration domain for various transverse radii r⊥. The inset in subfigure (b) is a sketch of the square integration domain, in which S denotes the position of the stationary point and the concentric circles are the oscillating periods of the integrand. The values of parameters used in the calculations are the same as those in Fig. 1.
Figure 3 presents the results for the Gaussian incident beam of Eq. (17), in which the errors of both E+ and E− are plotted as a function of the width dk of the square domain. All parameters used in the calculations are the same as those in Fig. 1. One can see that the errors of both E+ and E− decrease as the width dk is increasing. For a small domain with dk = 0.5, we already obtain an excellent accuracy, i.e., the errors for all values of transverse radius r⊥ are smaller than 2.0×10−5 for E+ and smaller than 4.0 × 10−5 for E−. The relative large errors of E− compared with E+ may come from the singularities in functions u− and v−, which make the numerical integration difficult. We emphasize that in contrast with the asymptotic expansion in the above section, the errors of E− in small r⊥ region are very small even for a quite small integration width dk; see Fig. 3(b) for details. This is basically due to the fact that the singularities in u− and v− have been included into the square domain of integration.
The accuracy of the above method can be further improved by a coarse-grained treatment. As the value of integral varies as the width of domain dk is increasing, we can do a coarse-grained average in one oscillating period of the integrand. For example, we can choose ten values of dk in a period, e.g., between the two squares in the inset of Fig. 4(b), calculate the integral for each dk, and treat the average as the final result. In this way, the highly oscillating part far away from the stationary point is canceled, leaving the most important contribution from the small area close to the stationary point. The calculated results for a Gaussian incident beam are shown in Fig. 4. One can see that the accuracy is improved greatly comparing with the Fig. 3, e.g., for a small width dk = 0.4, the errors for both E+ and E− are smaller than 1.0×10−5 for all radii r⊥. Note that for even smaller width such as dk = 0.2 the accuracy is still good.
Fig. 4 The absolute errors of |E+| (a) and |E−| (b) as a function of the width dk of the integration domain for various transverse radii r⊥. A coarse-grained average is performed between two square domains of integration; see the inset in subfigure (b) for a sketch. The values of parameters used in the calculations are the same as those in Fig. 1.
Therefore, we conclude that the method of integrating near the important points is very efficient and numerically cheap. It allows one to evaluate any kinds of wave propagation along the crystal axis, provided that the Fourier form of the incident beam is known. Moreover, this method is not restricted to the paraxial approximation and can be easily extended to compute the wave propagation in biaxial crystals.
4. The anisotropic propagating dynamics in biaxial crystals
For biaxial crystals, the expressions for functions f1 and f2 in Eq. (22) are complicated. Nevertheless, one can numerically show that f1 and f2 still have only one stationary point, whose position can also be obtained numerically. With these points, the method of integrating near important points can be straightforward applied to compute the optical field in biaxial crystals. The only difference from uniaxial crystals is that the domain of integration is no longer square, but is rectangle due to the breaking symmetry in the transverse plane. In the following, we adopt this method to investigate the propagation of vortex beam in the biaxial crystals.
We consider the case that an elliptically polarized Gaussian beam containing a single vortex is incident at the z = 0 plane,
where the two-dimensional complex vector d̂in = êx cos(θ)+êy sin(θ)eiϕ representing the state of polarization. The Fourier transformation of incident field Eq. (28) is In order to study the properties of propagation in biaxial crystals, we fix the values nx = 1.656 and nz = 1.458, while let ny vary in the range between nx and nz. Correspondingly, we introduce parameters Δn = nx − ny and η = Δn/(nx − nz) to describe the degree of deviation from the uniaxial crystals. The cases of η = 0 and η = 1 correspond to two limits of uniaxial crystal with propagation direction parallel to and orthogonal to the optical axis respectively (nx = ny ≠ nz and nx ≠ ny = nz). In the past, the propagating behaviors of vortex beam of Eq. (28) in these two uniaxial limits have been extensively studied [4–11]. What we are interested in this paper is the physics happens between these two limits. We firstly consider the case when the incident vortex beam is circularly polarized, e.g., d̂in = ê+. The beam parameters are still chosen as: λ = 0.633μm and s = 4.59μm, so that the paraxial condition k0s ≫ 1 is satisfied. We investigate the evolutions of Cartesian components Ex and Ey with different values of deviations η at a fixed propagating distance z = 5000μm. The results of numerical simulations are presented in Fig. 5.Fig. 5 The profiles of |Ex| (first row), |Ey| (second row), and intensity I = |Ex|2 + |Ey|2 (third row) of a circularly polarized Gaussian vortex beam propagating in biaxial crystals at a distance z = 5000μm. The deviations for each subfigures are: (a), (f), and (k): Δn = 0, η = 0; (b), (g), and (l): Δn = 0.0001, η = 0.0005; (c), (h), and (m): Δn = 0.001, η = 0.0051; (d), (i), and (n): Δn = 0.05, η = 0.2525; (e), (j), and (o): Δn = 0.198, η = 1. The x and y coordinates are in unit of μm. The colorbars for each rows are shown on the right sides.
One evident property of beam propagation in Fig. 5 is that both the Ex and Ey components for biaxial crystals exhibit anisotropic propagating behaviors due to the broken symmetry in the transverse plane. However, the anisotropic characters of the Ex and Ey components are quite different. The profile of Ex is extended in the x direction, nevertheless the profile of Ey is extended in the y direction. In the region where the deviation η is very small (e.g., the panels (b) and (g), (c) and (h)), the regular pattern formed in the uniaxial limit with η = 0 [5, 49] is destroyed quickly as the value of η is increasing. The propagation and angular momentum dynamics in this region are studied in detail before [33], where a formation of vortex pairs in the circular polarized component are observed. For larger value of η, the propagating behaviors of Ex and Ey components can be well understood by an analytical approximation.
For a paraxial incident field, it’s transverse width s is much greater than the vacuum wave length λ. Thus, the incident field Ẽ⊥ (0) in the Fourier space is nonvanishing only in a very small region with |k⊥| ≪ k0. For biaxial crystals with strong birefringence, e.g. η > 0.05 as in our calculations, we expect the following relation is satisfied for a paraxial beam,
In this limit, the complicated optical field inside biaxial crystals can be greatly simplified. We first look at the expression of L in Eq. (13). Since the dominate term of L is (α − β)2 when the value of η is not very small, it’s reasonable to change the form of terms in it. In order to simplify the expression, we therefore add to L a term proportional to , . Note that this approximate treatment is exact in the uniaxial limit of η = 1, for the added term is zero when β = γ. After the treatment, the expression of L in Eq. (13) can be written as which is actually the square of M in Eq. (14), L = M2. Then the eigenvalues λ1, λ2 and the matrix P, Q can be simplified as For the case nx ≠ ny = nz (i.e., β = γ and η = 1), the matrix P and Q reduce to be which is exactly the same as those in Ref. [23]. Following the procedures in Refs. [20, 23], we can expand the phase and amplitude of optical fields up to the second and the zeroth order respectively, in terms of the small parameters kx/k0 and ky/k0. Finally, we get The optical field in Eqs. (35) and (36) is a generalization of the results in Ref. [23], where the uniaxial limit of η = 1 is discussed. Many features of the beam propagation in Fig. 5 can be explained by employing these equations.An important property of the optical field described by Eqs. (35) and (36) is that two Cartesian components Ex and Ey are uncoupled, that is, the Ex and Ey depend on the x and y components of the incident field respectively, and the polarization conversion between them is forbidden during propagation.
For a circularly polarized incident beam, the two Cartesian components of it carry equal optical power. One can see from Fig. 5 that, for the figures with η > 0.05, the total power of Ex and Ey is almost equal and do not change as propagating, which is in agreement with the predictions from Eqs. (35) and (36).
The general features of anisotropy in Fig. 5 can also be understood. The Eq. (35) shows that the anisotropy of Ex component comes from the exponential factor inside the integrand. Note that the incident components Ẽx(0) and Ẽy(0) are isotropic in the transverse plane. As in our calculations, the value of nx = 1.656 is larger than that of nz = 1.458, so that the profile of Ex will extend in the x direction. Another interesting feature of Ex in Eq. (35) is that it does not depend on the value of ny, which implies that its profile will not change as the value of η changes. This feature can be seen from the panels (d) and (e) in Fig. 5, where the value of η is large enough. The component Ey exhibits a different anisotropic behavior from that of Ex. The anisotropy is caused by a different factor in the integrand, and therefore the profile is extended in y direction in case when the value ny > nz. Moreover, the anisotropy in the profile of Ey loses gradually as the value of ny approaches to that of nz (η approaches to 1). At the uniaxial limit of η = 1, Ey turns out to be isotropic in the transverse plane, see the panel (j) in Fig. 5.
To further illustrate the power exchange between two Cartesian components, we show in Fig. 6 the propagation of a linearly polarized Gaussian vortex beam at a distance z = 5000μm, i.e., d̂in in Eq. (28) is chosen to be êx. For η = 0, a significant amount of power is transferred from Ex to Ey; see Fig. 6(a) and (f). The power conversion in this limit has been well studied both experimentally and theoretically in Refs. [5, 25, 49, 50] for the case of Gaussian beam. As the value of η is increasing, the amount of transferred power decreases sharply, and for large η the intensity of Ey is much smaller than that of Ex, implying that the components Ex and Ey are almost decoupled. This agrees with our previous analysis on Eqs. (35) and (36). Note that the colorbar on the right side of the second row is only for subfigures (i) and (j). The colorbar for subfigures (f), (g), and (h) is the same as that of subfigures in the first row, which is on the right side of it. For completeness, the total intensity of the transverse fields I = |Ex|2 + |Ey|2 are shown in the third rows in Fig. 5 and Fig. 6 for circularly and linearly polarized beam, respectively. Regular pattern forms when η = 0, which is destroyed quickly as the value of η is increasing. For large η, the intensity I basically forms a elliptical profile and does not change a lot as η varies.
Fig. 6 The profiles of |Ex| (first row), |Ey| (second row), and intensity I = |Ex|2 + |Ey|2 (third row) of a linearly polarized Gaussian vortex beam propagating in biaxial crystals at a distance z = 5000μm. The deviations for each subfigures are: (a), (f), and (k): Δn = 0, η = 0; (b), (g), and (l): Δn = 0.0001, η = 0.0005; (c), (h), and (m): Δn = 0.001, η = 0.0051; (d), (i), and (n): Δn = 0.05, η = 0.2525; (e), (j), and (o): Δn = 0.198, η = 1. The x and y coordinates are in unit of μm. The colorbars for the first and third rows are shown on the right sides of them. In the second row, the colorbar of subfigures (f), (g), and (h) is the same as that of the first row, and the colorbar on the right side is only for subfigures (i) and (j).
5. Conclutions
We have proposed a numerical method to calculate the optical fields inside biaxial crystals in the plane-wave angular spectrum representation, which is based on the asymptotic expansion theory. The asymptotic approximation to the optical fields is examined and its advantages and disadvantages are discussed. We showed that the method of asymptotic approximation is not good when the highly oscillating integrand includes singular points. In order to overcome this shortcoming, we proposed the method of integrating near important points. By benchmarking the method against the exact solutions in uniaxial limits, we verified that it is very efficient and numerically easy to implement. Another asset is that it allows to compute optical fields inside a biaxial crystal with arbitrary incident beam.
In the second part of the paper, we adopted this method of integrating near important points to investigate the propagation of a Gaussian vortex beam in biaxial crystals. We focused on the anisotropic dynamics between two uniaxial limits with η = 0 and η = 1. We found that for the paraxial limit and a large value of η (the condition is satisfied), the two Cartesian components Ex and Ey are uncoupled and the power exchange between them is very small. The general features of anisotropic propagation are identified and explained by an analytical approximation.
It is noted that our analysis focuses on the case when the incident beam propagates along the crystallographic principal axes. However, the angular spectrum representation as well as our method in principle can be extended to the case of beam propagation along arbitrary directions [44]. To this end, we need to make a proper rotation from the crystalline coordinates to the experimental one in advance. Assuming an arbitrary propagation direction is specified by (θ, ϕ) in the crystalline coordinates, where θ is the polar angle and ϕ the azimuthal angle. Then the transformation matrix can be given as,
In this perspective, the diagonal dielectric tensor of Eq. (2) will become ε′ = TεTt in the experimental frame. It thus allows one to make our numerical description to the well-known phenomena of conical diffraction when the light propagates along the optic axis or in vicinity of “the diabolic points” [37, 43]. We would like to leave this possibility to our future studies.Acknowledgments
The authors acknowledge the financial support from the National Natural Science Foundation of China (NSFC) (Grant No. 11004164 and 11104233), the Doctoral Fund of Ministry of Education of China (Grant No. 2011012112003), the Fundamental Research Funds for the Central Universities (Grant No. 2011121043 and 2012121015) and the Natural Science Foundation of Fujian Province of China (Grant No. 2011J05010 and 2012J01285).
References and links
1. J. F. Nye, Natural Focusing and Fine Structure of Light (Institute of Physics Publishing, Bristol, 1999).
2. L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes,” Phys. Rev. A 45, 8185–8189 (1992) [CrossRef] [PubMed] .
3. M. S. Soskin, V. N. Gorshkov, M. V. Vasnetsov, J. T. Malos, and N. R. Heckenberg, “Topological charge and angular momentum of light beams carrying optical vortices,” Phys. Rev. A 56, 4064–4075 (1997) [CrossRef] .
4. G. Cincotti, A. Ciattoni, and C. Palma, “Laguerre-Gauss and Bessel-Gauss beams in uniaxial crystals,” J. Opt. Soc. Am. A 19, 1680–1688 (2002) [CrossRef] .
5. A. Ciattoni, G. Cincotti, and C. Palma, “Circularly polarized beams and vortex generation in uniaxial media,” J. Opt. Soc. Am. A 20, 163–171 (2003) [CrossRef] .
6. F. Flossmann, U. T. Schwarz, M. Maier, and M. R. Dennis, “Polarization singularities from unfolding an optical vortex through a birefringent crystal,” Phys. Rev. Lett. 95, 253901 (2005) [CrossRef] [PubMed] .
7. T. Fadeyeva, A. Rubass, Y. Egorov, A. Volyar, and J. Grover Swartzlander, “Quadrefringence of optical vortices in a uniaxial crystal,” J. Opt. Soc. Am. A 25, 1634–1641 (2008) [CrossRef] .
8. T. A. Fadeyeva, A. F. Rubass, and A. V. Volyar, “Transverse shift of a high-order paraxial vortex-beam induced by a homogeneous anisotropic medium,” Phys. Rev. A 79, 053815 (2009) [CrossRef] .
9. T. Fadeyeva, C. Alexeyev, B. Sokolenko, M. Kudryavtseva, and A. Volyar, “Non-canonical propagation of high-order elliptic vortex beams in a uniaxially anisotropic medium,” Ukr. J. Phys. Opt. 12, 62–82 (2011) [CrossRef] .
10. T. Fadeyeva, C. Alexeyev, A. Rubass, A. Zinov’ev, V. Konovalenko, and A. Volyar, “Subwave spikes of the orbital angular momentum of the vortex beams in a uniaxial crystal,” Opt. Lett. 36, 4215–4217 (2011) [CrossRef] [PubMed] .
11. J. Li, Y. Xin, Y. Chen, S. Xu, Y. Wang, M. Zhou, Q. Zhao, and F. Chen, “Diffraction of Gaussian vortex beam in uniaxial crystals orthogonal to the optical axis,” Eur. Phys. J. Appl. Phys. 53, 20701 (2011) [CrossRef] .
12. C. N. Alexeyev, “Circular array of anisotropic fibers: A discrete analog of a q plate,” Phys. Rev. A 86, 063830 (2012) [CrossRef] .
13. T. A. Fadeyeva, C. N. Alexeyev, P. M. Anischenko, and A. V. Volyar, “Engineering of the space-variant linear polarization of vortex-beams in biaxially induced crystals,” Appl. Opt. 51, C224–C230 (2012) [CrossRef] [PubMed] .
14. T. Fadeyeva, C. Alexeyev, P. Anischenko, and A. Volyar, “The fine structure of the vortex-beams in the biaxial and biaxially-induced birefringent media caused by the conical diffraction,” arXiv:1107.5775 (2011).
15. M. Born and E. Wolf, Principles of Optics (Pergamon, Oxford, 1999).
16. H. C. Chen, Theory of Electromagnetic Waves (McGraw-Hill, New York, 1983).
17. J. J. Stamnes and G. C. Sherman, “Radiation of electromagnetic fields in uniaxially anisotropic media,” J. Opt. Soc. Am. 66, 780–788 (1976) [CrossRef] .
18. J. J. A. Fleck and M. D. Feit, “Beam propagation in uniaxial anisotropic media,” J. Opt. Soc. Am. 73, 920–926 (1983) [CrossRef] .
19. M. Trippenbach and Y. B. Band, “Propagation of light pulses in nonisotropic media,” J. Opt. Soc. Am. B 13, 1403–1411 (1996) [CrossRef] .
20. A. Ciattoni, B. Crosignani, and P. D. Porto, “Vectorial theory of propagation in uniaxially anisotropic media,” J. Opt. Soc. Am. A 18, 1656–1661 (2001) [CrossRef] .
21. A. Volyar and T. Fadeeva, “Generation of singular beams in uniaxial crystals,” Optics and Spectroscopy 94, 235–244 (2003) [CrossRef] .
22. A. Volyar and T. Fadeeva, “Laguerre-Gaussian beams with complex and real arguments in a uniaxial crystal,” Opt. Spectros. 101, 450–457 (2006) [CrossRef] .
23. A. Ciattoni and C. Palma, “Optical propagation in uniaxial crystals orthogonal to the optical axis: paraxial theory and beyond,” J. Opt. Soc. Am. A 20, 2163–2171 (2003) [CrossRef] .
24. A. Ciattoni, G. Cincotti, and C. Palma, “Angular momentum dynamics of a paraxial beam in a uniaxial crystal,” Phys. Rev. E 67, 036618 (2003) [CrossRef] .
25. A. Ciattoni, G. Cincotti, C. Palma, and H. Weber, “Energy exchange between the Cartesian components of a paraxial beam in a uniaxial crystal,” J. Opt. Soc. Am. A 19, 1894–1900 (2002) [CrossRef] .
26. A. Ciattoni and C. Palma, “Anisotropic beam spreading in uniaxial crystals,” Opt. Commun. 231, 79–92 (2004) [CrossRef] .
27. A. Ciattoni, G. Cincotti, D. Provenziani, and C. Palma, “Paraxial propagation along the optical axis of a uniaxial medium,” Phys. Rev. E 66, 036614 (2002) [CrossRef] .
28. S. Olver, Numerical Approximation of Highly Oscillatory Integrals (PhD thesis, 2008).
29. D. Huybrechs and S. Olver, “Highly oscillatory quadrature,” in Highly Oscillatory Problems (Cambridge University, 2009), pp. 25–50 [CrossRef] .
30. A. Iserles and S. P. Nørsett, “Efficient quadrature of highly oscillatory integrals using derivatives,” Proc. R. Soc. A 461, 1383–1399 (2005) [CrossRef] .
31. S. Olver, “Moment-free numerical integration of highly oscillatory functions,” IMA J. Numer. Anal. 26, 213–227 (2006) [CrossRef] .
32. D. Huybrechs and S. Vandewalle, “On the evaluation of highly oscillatory integrals by analytic continuation,” SIAM J. Numer. Anal. 44, 1026–1048 (2006) [CrossRef] .
33. X. Lu and L. Chen, “Spin-orbit interactions of a Gaussian light propagating in biaxial crystals,” Opt. Express 20, 11753–11766 (2012) [CrossRef] [PubMed] .
34. R. Piessens, E. deDoncker Kapenga, C. Uberhuber, and D. Kahaner, Quadpack: A Subroutine Package for Automatic Integration (Springer Verlag, 1983).
35. N. S. Kazak, N. A. Khilo, and A. A. Ryzhevich, “Generation of Bessel light beams under the conditions of internal conical refraction,” Quantum Electron. 29, 1020 (1999) [CrossRef] .
36. A. Belsky and M. Stepanov, “Internal conical refraction of coherent light beams,” Opt. Commun. 167, 1–5 (1999) [CrossRef] .
37. V. N. Belyi, T. A. King, N. S. Kazak, N. A. Khilo, E. G. Katranji, and A. A. Ryzhevich, “Methods of formation and nonlinear conversion of Bessel optical vortices,” Proc. SPIE 4403, 229–240 (2001) [CrossRef] .
38. M. V. Berry, “Conical diffraction asymptotics: fine structure of Poggendorff rings and axial spike,” J. Opt. A, Pure Appl. Opt. 6, 289 (2004) [CrossRef] .
39. M. V. Berry, M. R. Jeffrey, and M. Mansuripur, “Orbital and spin angular momentum in conical diffraction,” J. Opt. A: Pure Appl. Opt. 7, 685 (2005) [CrossRef] .
40. M. V. Berry and M. R. Jeffrey, “Conical diffraction: Hamilton’s diabolical point at the heart of crystal optics,” in Progress in Optics, E. Wolf, ed. (Elsevier, 2007), pp. 13–50 [CrossRef] .
41. M. V. Berry, “Conical diffraction from an N-crystal cascade,” J. Opt. 12, 075704 (2010) [CrossRef] .
42. D. P. O’Dwyer, C. F. Phelan, Y. P. Rakovich, P. R. Eastham, J. G. Lunney, and J. F. Donegan, “The creation and annihilation of optical vortices using cascade conical diffraction,” Opt. Express 19, 2580–2588 (2011) [CrossRef] .
43. N. A. Khilo, “Conical diffraction and transformation of Bessel beams in biaxial crystals,” Opt. Commun. 286, 1–5 (2013) [CrossRef] .
44. M. A. Dreger, “Optical beam propagation in biaxial crystals,” J. Opt. A: Pure Appl. Opt. 1, 601 (1999) [CrossRef] .
45. Z. Chen and Q. Guo, “Rotation of elliptic optical beams in anisotropic media,” Opt. Commun. 284, 3183–3191 (2011) [CrossRef] .
46. T. A. Fadeyeva and A. V. Volyar, “Extreme spin-orbit coupling in crystal-traveling paraxial beams,” J. Opt. Soc. Am. A 27, 381–389 (2010) [CrossRef] .
47. T. A. Fadeyeva, V. G. Shvedov, Y. V. Izdebskaya, A. V. Volyar, E. Brasselet, D. N. Neshev, W. Krolikowski, and Y. S. Kivshar, “Spatially engineered polarization states and optical vortices in uniaxial crystals,” Opt. Express 18, 10848–10863 (2010) [CrossRef] [PubMed] .
48. R. Wong, Asymptotic Approximations of Integrals (SIAM, Philadelphia, 2001) [CrossRef] .
49. E. Brasselet, Y. Izdebskaya, V. Shvedov, A. S. Desyatnikov, W. Krolikowski, and Y. S. Kivshar, “Dynamics of optical spin-orbit coupling in uniaxial crystals,” Opt. Lett. 34, 1021–1023 (2009) [CrossRef] [PubMed] .
50. Y. Izdebskaya, E. Brasselet, V. Shvedov, A. Desyatnikov, W. Krolikowski, and Y. Kivshar, “Dynamics of linear polarization conversion in uniaxial crystals,” Opt. Express 17, 18196–18208 (2009) [CrossRef] [PubMed] .
