## Abstract

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 OSA

## 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],

**E**(

**r**) is the complex amplitude of electric field

**E**(

**r**,

*t*) with a relation

**E**(

**r**,

*t*) =

*Re*[

**E**(

**r**)exp(−

*iωt*)],

*k*

_{0}= 2

*π*/

*λ*with

*λ*being the vacuum wave length of light, and

*ε*is the relative dielectric tensor which has form

*n*,

_{x}*n*, and

_{y}*n*being the refractive index along crystalline

_{z}*x*,

*y*, and

*z*axis, respectively. For the uniaxial crystals with

*n*=

_{x}*n*≠

_{y}*n*or

_{z}*n*≠

_{x}*n*=

_{y}*n*, 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

_{z}*n*≠

_{x}*n*≠

_{y}*n*. We assume that the light is propagating along the

_{z}*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],

**r**

_{⊥}=

*x*

**ê**

*+*

_{x}*y*

**ê**

*and*

_{y}**k**

_{⊥}=

*k*

_{x}**ê**

*+*

_{x}*k*

_{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*

_{y}**Ẽ**(

**k**

_{⊥},

*z*),

*λ*

_{1},

*λ*

_{2},

*c*

_{1},

*c*

_{2},

*c*

_{3},

*c*

_{4},

*c*

_{5}, and

*c*

_{6}are given by [33]

**P**and

**Q**are

**Ẽ**

_{⊥}(0) is the two-dimensional Fourier transformation of the transverse incident field

**E**

_{⊥}(

**r**

_{⊥}, 0) on the plane

*z*= 0,

*α*,

*β*,

*γ*,

*p*,

*K*,

*L*,

*M*are defined to be

The 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 *n _{x}* =

*n*=

_{y}*n*

_{0}and

*n*=

_{z}*n*. The expressions for the coefficients

_{e}*λ*

_{1}(

*λ*

_{2}) and the matrices

**P**(

**Q**) in uniaxial crystals can be simplified as [20],

*z*= 0 is

*s*is the waist of Gaussian beam and ${\widehat{\mathbf{e}}}_{+}=\frac{1}{\sqrt{2}}\left({\widehat{\mathbf{e}}}_{x}+i{\widehat{\mathbf{e}}}_{y}\right)$, which as well as its partner ${\widehat{\mathbf{e}}}_{-}=\frac{1}{\sqrt{2}}\left({\widehat{\mathbf{e}}}_{x}-i{\widehat{\mathbf{e}}}_{y}\right)$ 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**

_{⊥}=

*E*

_{x}**ê**

*+*

_{x}*E*

_{y}**ê**

*=*

_{y}*E*

_{+}

**ê**

_{+}+

*E*

_{−}

**ê**

_{−}. The relation between different components is

*k*space is a Fourier transformation of Eq. (17), which is given by

*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],

*λ*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 $\frac{\partial f}{\partial x}=\frac{\partial f}{\partial y}=0$. 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 *E _{x}*,

*z*,

*E*consists of two typical highly oscillatory integral, which can be approximated by asymptotic expansion. It’s easy to show that both functions

_{x}*f*

_{1}and

*f*

_{2}have only one stationary point, which in the case of uniaxial crystals are at

*r*= (

_{o}*x*

^{2}+

*y*

^{2}+

*z*

^{2})

^{1/2}, ${r}_{e}={\left[{x}^{2}+{y}^{2}+{\left(\frac{{n}_{o}}{{n}_{e}}\right)}^{2}{z}^{2}\right]}^{1/2}$. For the dominated contributions to the integral in

*E*come from the critical points of first kind, one can expand the functions

_{x}*f*

_{1}and

*f*

_{2}on the critical points (

*k*,

_{xo}*k*) and (

_{yo}*k*,

_{xe}*k*), respectively. For example, function

_{ye}*f*

_{1}can be expanded as

*f*

_{2}, ${a}_{2}={\frac{{\partial}^{2}{f}_{2}}{\partial {k}_{x}^{2}}|}_{\left({k}_{xe},{k}_{ye}\right)}$, ${b}_{2}={\frac{{\partial}^{2}{f}_{2}}{\partial {k}_{y}^{2}}|}_{\left({k}_{xe},{k}_{ye}\right)}$, and ${c}_{2}={\frac{{\partial}^{2}{f}_{2}}{\partial {k}_{x}\partial {k}_{y}}|}_{\left({k}_{xe},{k}_{ye}\right)}$. After truncated the Taylor expansion of

*f*

_{1}and

*f*

_{2}to second order and substituting them into Eq. (21), we can derive the asymptotic approximation to

*E*,

_{x}*c*

_{1}(

*k*,

_{x}*k*) and

_{y}*c*

_{2}(

*k*,

_{x}*k*) in the integral by constants

_{y}*c*

_{1}(

*k*,

_{xo}*k*) and

_{yo}*c*

_{2}(

*k*,

_{xe}*k*), which can be evaluated by substituting (

_{ye}*k*,

_{xo}*k*) and (

_{yo}*k*,

_{xe}*k*) into Eq. (6). However, this approximation is valid only when the functions

_{ye}*c*

_{1}and

*c*

_{2}are differentiable [48]. We will discuss this point in detail latter.

The asymptotic expansion for the other two components, *E _{y}* and

*E*, can be obtained in a similar way, by replacing

_{z}*c*

_{1},

*c*

_{2}with

*c*

_{3},

*c*

_{4}and

*c*

_{5},

*c*

_{6}, respectively. After some straightforward algebra, the electric field

**E**, in the case of uniaxial crystal, can be simplified as

After computing the values of function *c _{i}* (

*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

*n*= 1.656,

_{o}*n*= 1.458, the wave length

_{e}*λ*= 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.

The asymptotic approximation of Eq. (27) is obtained by assuming that the functions *c _{i}*(

*k*,

_{x}*k*) (

_{y}*i*= 1, ⋯,6) are smooth and differentiable [48]. For circular components

*E*

_{±}, one should look at the behavior of functions

*u*

_{±}=

*c*

_{1}±

*ic*

_{3}and

*v*

_{±}=

*c*

_{2}±

*ic*

_{4}, 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 (

*k*,

_{xo}*k*) and (

_{yo}*k*,

_{xe}*k*) are very close to the singular point (0, 0) in

_{ye}*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).

Therefore, we conclude that as long as the functions *c _{i}* are smooth, the accuracy of asymptotic approximation is great for large propagating distance

*z*. Note that the functions

*c*depend on the incident beams, which often contains singularities. In the following, we propose a method to take the effect of singularity into account.

_{i}#### 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 *f*_{1} and *f*_{2}, and also the singular points in functions *c _{i}*. 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 (*k _{xo}*,

*k*) and (

_{yo}*k*,

_{xe}*k*). 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.

_{ye}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 *d _{k}* 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

*d*is increasing. For a small domain with

_{k}*d*= 0.5, we already obtain an excellent accuracy, i.e., the errors for all values of transverse radius

_{k}*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

*d*; see Fig. 3(b) for details. This is basically due to the fact that the singularities in

_{k}*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 *d _{k}* is increasing, we can do a coarse-grained average in one oscillating period of the integrand. For example, we can choose ten values of

*d*in a period, e.g., between the two squares in the inset of Fig. 4(b), calculate the integral for each

_{k}*d*, 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

_{k}*d*= 0.4, the errors for both

_{k}*E*

_{+}and

*E*

_{−}are smaller than 1.0×10

^{−5}for all radii

*r*

_{⊥}. Note that for even smaller width such as

*d*= 0.2 the accuracy is still good.

_{k}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 *f*_{1} and *f*_{2} in Eq. (22) are complicated. Nevertheless, one can numerically show that *f*_{1} and *f*_{2} 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,

**d̂**

*=*

_{in}**ê**

*cos(*

_{x}*θ*)+

**ê**

*sin(*

_{y}*θ*)

*e*representing the state of polarization. The Fourier transformation of incident field Eq. (28) is

^{iϕ}*n*= 1.656 and

_{x}*n*= 1.458, while let

_{z}*n*vary in the range between

_{y}*n*and

_{x}*n*. Correspondingly, we introduce parameters Δ

_{z}*n*=

*n*−

_{x}*n*and

_{y}*η*= Δ

*n*/(

*n*−

_{x}*n*) to describe the degree of deviation from the uniaxial crystals. The cases of

_{z}*η*= 0 and

*η*= 1 correspond to two limits of uniaxial crystal with propagation direction parallel to and orthogonal to the optical axis respectively (

*n*=

_{x}*n*≠

_{y}*n*and

_{z}*n*≠

_{x}*n*=

_{y}*n*). 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.,

_{z}**d̂**

*=*

_{in}**ê**

_{+}. The beam parameters are still chosen as:

*λ*= 0.633

*μm*and

*s*= 4.59

*μm*, so that the paraxial condition

*k*

_{0}

*s*≫ 1 is satisfied. We investigate the evolutions of Cartesian components

*E*and

_{x}*E*with different values of deviations

_{y}*η*at a fixed propagating distance

*z*= 5000

*μm*. The results of numerical simulations are presented in Fig. 5.

One evident property of beam propagation in Fig. 5 is that both the *E _{x}* and

*E*components for biaxial crystals exhibit anisotropic propagating behaviors due to the broken symmetry in the transverse plane. However, the anisotropic characters of the

_{y}*E*and

_{x}*E*components are quite different. The profile of

_{y}*E*is extended in the

_{x}*x*direction, nevertheless the profile of

*E*is extended in the

_{y}*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

*E*and

_{x}*E*components can be well understood by an analytical approximation.

_{y}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**_{⊥}| ≪ *k*_{0}. 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,

*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 ${\mathbf{k}}_{\perp}^{4}$ terms in it. In order to simplify the expression, we therefore add to

*L*a term proportional to ${\mathbf{k}}_{\perp}^{4}$, $4\left(1-\frac{\alpha}{\gamma}\right)\left(1-\frac{\beta}{\gamma}\right){k}_{x}^{2}{k}_{y}^{2}$. 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

*M*in Eq. (14),

*L*=

*M*

^{2}. Then the eigenvalues

*λ*

_{1},

*λ*

_{2}and the matrix

**P**,

**Q**can be simplified as

*n*≠

_{x}*n*=

_{y}*n*(i.e.,

_{z}*β*=

*γ*and

*η*= 1), the matrix

**P**and

**Q**reduce to be

*k*/

_{x}*k*

_{0}and

*k*/

_{y}*k*

_{0}. Finally, we get

*η*= 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 *E _{x}* and

*E*are uncoupled, that is, the

_{y}*E*and

_{x}*E*depend on the

_{y}*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 *E _{x}* and

*E*is almost equal and do not change as propagating, which is in agreement with the predictions from Eqs. (35) and (36).

_{y}The general features of anisotropy in Fig. 5 can also be understood. The Eq. (35) shows that the anisotropy of *E _{x}* component comes from the exponential factor
${n}_{x}^{2}{k}_{x}^{2}+{n}_{z}^{2}{k}_{y}^{2}$ inside the integrand. Note that the incident components

*Ẽ*(0) and

_{x}*Ẽ*(0) are isotropic in the transverse plane. As in our calculations, the value of

_{y}*n*= 1.656 is larger than that of

_{x}*n*= 1.458, so that the profile of

_{z}*E*will extend in the

_{x}*x*direction. Another interesting feature of

*E*in Eq. (35) is that it does not depend on the value of

_{x}*n*, which implies that its profile will not change as the value of

_{y}*η*changes. This feature can be seen from the panels (d) and (e) in Fig. 5, where the value of

*η*is large enough. The component

*E*exhibits a different anisotropic behavior from that of

_{y}*E*. The anisotropy is caused by a different factor ${n}_{z}^{2}{k}_{x}^{2}+{n}_{y}^{2}{k}_{y}^{2}$ in the integrand, and therefore the profile is extended in

_{x}*y*direction in case when the value

*n*>

_{y}*n*. Moreover, the anisotropy in the profile of

_{z}*E*loses gradually as the value of

_{y}*n*approaches to that of

_{y}*n*(

_{z}*η*approaches to 1). At the uniaxial limit of

*η*= 1,

*E*turns out to be isotropic in the transverse plane, see the panel (j) in Fig. 5.

_{y}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

**ê**

*. For*

_{x}*η*= 0, a significant amount of power is transferred from

*E*to

_{x}*E*; 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

_{y}*η*is increasing, the amount of transferred power decreases sharply, and for large

*η*the intensity of

*E*is much smaller than that of

_{y}*E*, implying that the components

_{x}*E*and

_{x}*E*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

_{y}*I*= |

*E*|

_{x}^{2}+ |

*E*|

_{y}^{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.

## 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
$2{n}_{x}\left({n}_{x}-{n}_{z}\right)\eta {k}_{0}^{2}\gg {\mathbf{k}}_{\perp}^{2}$ is satisfied), the two Cartesian components *E _{x}* and

*E*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.

_{y}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,

*ε*′ =

*TεT*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.

^{t}## 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] .