Parabasal field decomposition and its application to non-paraxial propagation

Daniel Asoubar; Site Zhang; Frank Wyrowski; Michael Kuhn

doi:10.1364/OE.20.023502

1. Introduction

Nowadays field tracing enables optical modeling and design taking the wave nature of light into consideration [1]. This approach provides the tracing of electromagnetic fields through optical systems. An essential part of this simulation technique is the propagation of harmonic fields through homogeneous media. Common propagation techniques are based on the Fast Fourier Transformation (FFT) algorithm, due to its short computational time. A rigorous propagation technique is the Spectrum of Plane Waves (SPW) operator [2]. Here each harmonic field component V_ℓ(ρ, 0), given on a plane boundary ρ = (x, y)^T orthogonal to the propagation direction, is decomposed into a set of plane waves by the Fourier transformation

A_{ℓ} (κ, 0) = ℱ [V_{ℓ} (ρ, 0)] = \frac{1}{2 π} \iint_{- \infty}^{\infty} V_{ℓ} (ρ, 0) e^{- i κ \cdot ρ} d ρ .

Here ℓ = 1...6 denotes the 3 electric and the 3 magnetic field components of the harmonic field and κ = (k_x, k_y)^T the corresponding spatial frequency vector. Then the plane waves are propagated over the distance z and superimposed using an inverse Fourier transformation

V_{ℓ} (ρ, z) = ℱ^{- 1} [A_{ℓ} (κ, 0) e^{i k_{z} (κ) z}]

with the spherical phase function

k_{z} (κ) = \sqrt{k^{2} - k_{x}^{2} - k_{y}^{2}}

, leading to the propagated harmonic field components V_ℓ(ρ, z) in the parallel target plane. Please note that we have used the wavenumber in the homogeneous medium k = k₀n with refractive index n and vacuum wavenumber k₀. In practice just two harmonic field components have to be propagated and the remaining components can be calculated by Maxwell’s Equations on demand [1]. From the physical point of view the SPW operator is valid for all propagation distances z and all maximum spatial frequencies |κ_max|. Nevertheless for long propagation distances the numerical effort to sample Eq. (2) is too high. For non-paraxial fields, which are containing high spatial frequencies, the numerical effort is getting even worse which is mainly due to the high oscillating spherical phase k_z(κ). In Fig. 1 this limited scope of the SPW operator due to fast growing numerical effort is schematically sketched. In the past several approaches [3, 4, 5, 6, 7] have been tried to overcome this limited range. Nevertheless a general solution of the problem is not known.

Fig. 1 Schematic illustration of the scope for common FFT-based propagation operators in field tracing. Colored areas indicate parameter spaces for which the techniques are numerically feasible. The rigorous SPW operator suffers from high numerical effort for high propagation distances or spatial frequencies. The far field and Fresnel propagation operators have a limited range of validity due to their nature of physical approximations.

Download Full Size | PDF

Maybe the most famous approaches in literature are two approximations of Eq. (2). The Fresnel operator [8] is using a Taylor expansion of the spatial frequency component k_z(κ). In this way the spherical phase function is replaced by a parabolic one, leading to a semi-analytical solution of the inverse Fourier transformation integral of Eq. (2). As illustrated in Fig. 1 this concept is just valid for paraxial fields, which only have small spatial frequencies and sufficiently large propagation distances.

For long propagation distances z the far field operator [2] can be applied as illustrated in Fig. 1. Here the inverse Fourier transformation integral of Eq. (2) can be solved semi-analytically using the idea of stationary phase.

In the following we will develop a rigorous FFT based propagation operator to fill the remaining white areas in Fig. 1, using a combination of a parabasal decomposition technique (PDT) and a semi-analytical SPW propagation operator. For this purpose we will describe the fundamentals of parabasal fields first. In Chapter 3, we will derive a semi-analytical expression of Eq. (2) for parabasal fields. We will generalize this concept using a PDT to propagate non-paraxial fields rigorously in Chapter 4. All simulations were done with the optics software VirtualLab^TM [9].

2. Fundamentals of parabasal fields

An electromagnetic field is called parabasal one when the angular spectrum of each of its harmonic field components is located around a central spatial frequency vector κ₀. It means that a single parabasal harmonic field component possesses a low divergence and propagates with mean direction along

k = k {\hat{s}}_{0} = k {(s_{0 x}, s_{0 y}, s_{0 z})}^{T} = k {(sin θ_{0} cos ϕ_{0}, sin θ_{0} sin ϕ_{0}, cos θ_{0})}^{T}

defined in spherical coordinates θ and ϕ with unity radius. If we restrict to non-lossy dielectrics, the direction vector defines the central spatial frequency vector

κ_{0} = {(k_{0 x}, k_{0 y})}^{T} = k {(s_{0 x}, s_{0 y})}^{T}

with the wavenumber in the homogeneous medium k = k₀n. In the special case of κ₀ = (0, 0)^T the parabasal field propagates along the optical axis z and is called paraxial. In principle the following theory also can be applied on complex refractive indices. Then the central spatial frequency is given by the projection of 𝔕k on the k_x−k_y-plane. For simplification we will just discuss non-absorbing dielectrics in this paper.

Figure 2 illustrates the lower limit of a parabasal harmonic field with zero divergence, which is a plane wave pointing into the direction ŝ₀. For the sake of simplicity the two dimensional case where s_0y = 0 is chosen. In general a parabasal field can be written as

V_{ℓ} (ρ, 0) = V_{ℓ}^{'} (ρ, 0) e^{i κ_{0, ℓ} \cdot ρ}

which is meaning that it is always possible to split a parabasal harmonic field component into a residual field V′_ℓ(ρ, 0) and a linear phase term. Please note that in general the central spatial frequency κ_0,ℓ can be different for different components. For a given parabasal field V_ℓ(ρ, 0) we get the corresponding residual field using the shift theorem of the Fourier transformation

A_{ℓ} (κ, 0) = A_{ℓ}^{'} (κ - κ_{0, ℓ}, 0)

V_{ℓ}^{'} (ρ, 0) = ℱ^{- 1} [A_{ℓ}^{'} (κ, 0)] .

Fig. 2 Illustration of fundamental properties of a plane wave: Projecting the wavefront (dashed line in the left figure), which is pointing into the direction ŝ₀, on the xy-plane leads to a linear phase term (central figure). According to the shift theorem of the Fourier transformation this linear phase results in a spectrum (right figure) which is located around the spatial frequency κ₀.

Download Full Size | PDF

3. Semi-analytical SPW propagation operator

In the previous section it was shown that a parabasal field is located around a certain spatial frequency κ_0,ℓ. Using this property we can decompose k_z(κ) independently for each component ℓ around κ_0,ℓ by

k_{z, ℓ} (κ) = γ_{0, ℓ} + γ_{1, ℓ} \cdot κ + γ_{ℓ} (κ - κ_{0, ℓ})

with the constant term

γ_{0, ℓ} = k_{0 z, ℓ} + \frac{1}{k_{0 z, ℓ}} (k_{0 x, ℓ}^{2} + k_{0 y, ℓ}^{2}) - \frac{k_{0 x, ℓ}^{2} k_{0 y, ℓ}^{2}}{k_{0 z, ℓ}^{3}},

and the linear term

γ_{1, ℓ} = {[\frac{k_{0 x, ℓ}}{k_{0 z, ℓ}} (\frac{k_{0 y, ℓ}^{2}}{k_{0 z, ℓ}^{2}} - 1), \frac{k_{0 y, ℓ}}{k_{0 z, ℓ}} (\frac{k_{0 x, ℓ}^{2}}{k_{0 z, ℓ}^{2}} - 1)]}^{T}

with

k_{0 z, ℓ} = \sqrt{k^{2} - k_{0 x, ℓ}^{2} - k_{0 y, ℓ}^{2}} .

Here the residual term γ_ℓ(κ − κ_0,ℓ) denotes a function of (κ − κ_0,ℓ). We can calculate this residual term by rearranging Eq. (8) and performing a coordinate transformation into a new system κ′ = (k′_x, k′_y)^T with center at κ_0,ℓ. That results in

γ_{ℓ} (κ^{'}) = k_{z, ℓ} (κ^{'} + κ_{0, ℓ}) - γ_{0, ℓ} - γ_{1, ℓ} \cdot (κ^{'} + κ_{0, ℓ})

with

k_{z, ℓ} (κ^{'} + κ_{0, ℓ}) = \sqrt{k^{2} - {(k_{x}^{'} + k_{0 x, ℓ})}^{2} - {(k_{y}^{'} + k_{0 y, ℓ})}^{2}} .

Please note that the decomposition of k_z(κ) in Eq. (8) represents a Taylor expansion around κ₀ with the inclusion of all higher order terms in γ_ℓ(κ − κ₀). Plugging Eq. (8) into Eq. (2) of the SPW operator and applying the shift theorem of the Fourier transformation leads to

V_{ℓ} (ρ, z) = V_{ℓ}^{'} (ρ + z γ_{1, ℓ}, z) exp [i (γ_{0, ℓ} z + γ_{1, ℓ} \cdot κ_{0, ℓ} z + ρ \cdot κ_{0, ℓ})]

with the propagated shifted residual field

V_{ℓ}^{'} (ρ, z) = ℱ^{- 1} [A_{ℓ}^{'} (κ, 0) e^{i γ_{ℓ} (κ) z}]

where the initial parabasal field is written in the separated form of Eq. (5).

Please note that for parabasal fields around κ_0,ℓ there are two reasons lowering the sampling effort of Eq. (15) compared to Eq. (2). At first the phase of the residual field A′_ℓ(κ, 0) in spatial frequency domain has a flatter slope than the complete angular spectrum A_ℓ(κ, 0) due to analytical handling of the initial linear phase e^{(iκ_0,ℓρ)}. Secondly in Eq. (15) just the higher order terms γ_ℓ(κ) have to be sampled instead of the whole spherical phase k_z(κ) in Eq. (2). As shown in Fig. 3 the slope of γ_ℓ(κ − κ_0,ℓ) for κ around κ_0,ℓ is much flatter than the slope of k_z(κ). According to the sampling theorem [10] a flatter slope of the phase function will dramatically reduce the required sampling effort. We like to emphasize that the 2D case in Fig. 3 is just chosen for a convenient display without any restriction.

Fig. 3 2D line evaluation of the spherical phase k_z(κ) (blue line) and its decomposition (magenta and red line) according to Eq. (8) along the diagonal direction ${(k_{x}, k_{y})}^{T} = \frac{k_{ρ}}{\sqrt{2}} {(1, 1)}^{T}$ . Left side: Normalized spherical phase function k_z(κ) and linear phase function γ_0,ℓ + γ_1,ℓ · κ around $\frac{k_{0 ρ}}{k} = 0.424$ . Right side: Higher order phase term γ_ℓ(κ − κ_0,ℓ) for $\frac{k_{0 ρ}}{k} = 0.424$ calculated by Eq. (12).

Download Full Size | PDF

Figure 4 shows as an example a parabasal super Gaussian beam with an additional linear phase of 10° in x-direction and 4° in y-direction at a wavelength of 532 nm. The waist diameter is 50 μm and the beam is linearly polarized in x-direction.

Fig. 4 Example for a parabasal field: Amplitude (left) and phase (right) of a super Gaussian beam with 50 μm waist diameter at a wavelength of 532 nm. An additional linear phase term of 10° in x-direction and 4° in y-direction replaces the paraxial properties of the super Gaussian beam by a parabasal field behavior. The positions of the 1D cross section profiles are given by the black arrows.

Download Full Size | PDF

Table 1 shows the simulation accuracy and the numerical effort for the propagation of the parabasal super Gaussian beam over a distance z = 10 mm in vacuum using different propagation operators. The accuracy of the propagation which leads to V_ℓ_,operator(x, y, z) is given by its deviation

d : = \frac{\sum_{x, y} {| V_{ℓ, operator} (x, y, z) - V_{ℓ, reference} (x, y, z) |}^{2}}{\sum_{x, y} {| V_{ℓ, reference} (x, y, z) |}^{2}}

where the propagated reference field V_ℓ_,reference(x, y, z) is obtained by the rigorous SPW propagation operator. The numerical effort of the different propagation techniques is given by the number of sampling points which is needed to store the propagated field V_ℓ_,operator(x, y, z). We can see that the SPW operator requires high numerical effort. The Fresnel and the far field operators cannot be applied due to the high deviation caused by the approximative models. In comparison to the other propagation operators the rigorous semi-analytical SPW operator possesses a low deviation as well as a reduced numerical effort.

Table 1. Accuracy and numerical effort for different propagation techniques to propagate the parabasal super Gaussian beam by 10 mm:

View Table | View all tables in this article

The numerical effort in terms of sampling points only allows the estimation of required computer memory. There is no statement about the computational speed. We do not make a comparison of CPU times here because we have not made any effort to optimize the computer codes in this respect, yet. That is why the numerical complexity of the semi-analytical SPW propagation operator algorithm is investigated roughly only in what follows. We define a complexity factor as

η : = \frac{Complexity of SPW algorithm}{Complexity of semi-analytical SPW algorithm}

where we assume that the complexity of both algorithms can be estimated by the complexity of the 2D-FFT algorithm (O(N² log₂ N)) and the complexity of the multiplication (O(N²)) of phase terms with a harmonic field with N × N sampling points. Applying Eq. (17) on the previous propagation of the parabasal super Gaussian leads to a complexity factor η of 73.6 for the semi-analytical SPW operator.

4. Parabasal field decomposition technique

In the previous section the semi-analytical SPW operator was derived and applied to parabasal fields. In principle the same algorithm can be applied to general non-parabasal harmonic fields because of the rigorous treatment of higher order phase terms in the Taylor expansion of Eq. (8). However as shown in Fig. 3 for non-parabasal fields the sampling effort for the higher order phase terms in Eq. (15) will increase exorbitantly, resulting in a rapid increase of the numerical effort and a reduction of the complexity factor. To overcome this problem a parabasal field decomposition technique (PDT) will be introduced in the following. Figure 5 illustrates the combination of the PDT with the semi-analytical SPW operator to propagate non-paraxial fields rigorously and with reduced computational effort.

Fig. 5 Flowchart for the efficient propagation of non-paraxial fields using a combination of a parabasal decomposition technique (PDT) and the semi-analytical SPW operator: A general harmonic field is decomposed into a set of parabasal subfields using a PDT technique. After that each parabasal harmonic subfield is propagated by the semi-analytical SPW operator. Finally all propagated parabasal subfields are superimposed in the target plane leading to the propagated non-paraxial field.

Download Full Size | PDF

Figure 6 shows the decomposition of a Laguerre Gaussian beam of first order in the spatial frequency domain into 2 × 2 parabasal subfields.

Fig. 6 Decomposition of a first order Laguerre Gaussian beam into 2 × 2 parabasal subfields in the spatial frequency domain.

Download Full Size | PDF

Mathematically the decomposition into N subfields in the spatial frequency domain can be written as

A_{ℓ} (κ, 0) = \sum_{i, j = 1}^{N} Φ_{i, j} (κ) A_{ℓ} (κ, 0)

where Φ_i_,j(κ) is the basis function with local carrier and Φ_i,j(κ) := Φ(κ − κ_i,j) with a plenty of supporting points κ_i,j = (k_i,x, k_j,y)^T located e.g. on an equidistant rectangular grid. Also the property

\sum_{i, j = 1}^{N} Φ_{i, j} (κ) = 1

should be fulfilled for all points κ inside the definition area of the field A_ℓ(κ, 0). An obvious choice for Φ(κ) is the rect-function with Φ(κ) = 1 for max(k_x, k_y) < 1 and Φ(κ) = 0 anywhere else. As shown in Fig. 7 such basis function would result in a decomposition with sharp edges of the parabasal subfields. Then an inverse Fourier transformation will result in an infinite extended subfield in the space domain.

Fig. 7 Single parabasal subfield created by a decomposition of a Laguerre Gaussian beam of first order using a rect-basis function: The sharp edges of the rect-function in the spatial frequency domain (left) lead to an infinite extended subfield in the space domain (right).

Download Full Size | PDF

To overcome such numerically inconvenient sharp edges we define Φ(κ) for the one dimensional case in the following way:

\begin{array}{l} Φ (k_{x}) = 0.5 [sin (\frac{π}{a} (k_{x} - \frac{a}{2})) + 1] & for 0 < k_{x} \leq a, \\ Φ (k_{x}) = 1 & for a < k_{x} < b - a, \\ Φ (k_{x}) = 1 - 0.5 [sin (\frac{π}{a} (k_{x} - \frac{2 b - a}{2})) + 1] & for b - a \leq k_{x} < b, \\ and Φ (k_{x}) = 0 & anywhere else, \end{array}

with a, b ∈ ℝ and 0 < a < b. This definition can be generalized to the two-dimensional case using a separation Ansatz

Φ_{i, j} (κ) = Φ (k_{x} - k_{i, x}) Φ (k_{y} - k_{j, y})

and a desired scaling of the definition range b and edge width a can be implemented. The first and third line of the Eq. (19) represent the smooth edges of the basis function. In what follows these edges are called boundaries of the parabasal subfields. Figure 8 shows the same parabasal subfield as in Fig. 7 using the new basis function with smooth edges of Eq. (19). Please note that the basis function given in Eq. (19) is not the only possible candidate to truncate the parabasal subfields. However it is beyond the scope of this paper to compare different basis functions. We have chosen the particular function due to its simple implementation and well performance in practice.

Fig. 8 Single parabasal subfield created by a decomposition of a Laguerre Gaussian beam of first order using the basis function of Eq. (19): The smooth edges of the basis function in the spatial frequency domain (left) leads to a finitely extended subfield in the space domain (right).

Download Full Size | PDF

From the physical point of view all fields can be decomposed into parabasal fields in the spatial frequency domain due to the definition of a parabasal field in the spatial frequency domain. However this is not always convenient from a numerical point of view. It is useful to distinguish between two basic cases of non-paraxial fields:

The field is very divergent because of small features in the field function but it can be sampled without problems in the space domain. In this case the FFT algorithm can transform the field into the spatial frequency domain, where the PDT described above can be applied. A Gaussian beam with small waist or a strongly scattered field are examples of such fields of first kind.
The field possesses a smooth but strong phase function, which does not allow its sampling in space domain. Here a FFT algorithm cannot be applied, meaning that the spectrum in the spatial frequency domain is not accessible. In this case the PDT must be performed in the space domain. Such fields, e.g. a spherical or cylindrical wave usually are given in an analytical form. In what follows we will refer to such fields as fields of second kind.

Figure 9 illustrates the fundamental differences between the two basic cases of non-paraxial fields.

Fig. 9 Illustration of two basic cases of non-paraxial fields. On the left-hand side a wavefront is illustrated, which possesses rough wavefront features. Then only a PDT in spatial frequency domain is reasonable. On the right hand-side a wavefront is shown, which is smooth but has a strong phase function. Adjacent positions of the wavefront contribute to almost the same spatial frequencies which enable a PDT in space domain.

Download Full Size | PDF

So far the shift theorem of the Fourier transformation has been used to extract linear phase terms from sampled parabasal fields (Eq. (6) and (7)). However for fields of second kind, which are given analytically, there is another way to bring them into the form of Eq. (5). In principle any analytically stored harmonic field component can be written in the form

V_{ℓ} (ρ, 0) = {\hat{V}}_{ℓ} (ρ, 0) e^{i φ_{ℓ} (ρ, 0)}

where V_ℓ(ρ, 0) is representing the amplitude and φ_ℓ(ρ, 0) the phase function. Using the Taylor expansion around the parabasal field center ρ_0,ℓ in space domain leads to

{φ (ρ, 0) |}_{(ρ \approx ρ_{0, ℓ})} = ε_{ℓ} + κ_{0, ℓ} \cdot ρ + Δ_{ℓ} (ρ, 0)

with

ε_{ℓ} = φ_{ℓ} (ρ_{0, ℓ}, 0) + {\frac{\partial^{2} φ_{ℓ}}{\partial x \partial y} |}_{ρ_{0, ℓ}} \cdot (x_{0, ℓ} y_{0, ℓ}) - {\frac{\partial φ_{ℓ}}{\partial x} |}_{ρ_{0, ℓ}} \cdot x_{0, ℓ} - {\frac{\partial φ_{ℓ}}{\partial y} |}_{ρ_{0, ℓ}} \cdot y_{0, ℓ}

and

κ_{0, ℓ} = {(k_{0 x, ℓ}, k_{0 y, ℓ})}^{T} = {({\frac{\partial φ_{ℓ}}{\partial x} |}_{ρ_{0, ℓ}} - {\frac{\partial^{2} φ_{ℓ}}{\partial x \partial y} |}_{ρ_{0, ℓ}} \cdot y_{0, ℓ}, {\frac{\partial φ_{ℓ}}{\partial y} |}_{ρ_{0, ℓ}} - {\frac{\partial^{2} φ_{ℓ}}{\partial x \partial y} |}_{ρ_{0, ℓ}} \cdot x_{0, ℓ})}^{T} .

Please note that Δ_ℓ(ρ, 0) again contains all higher order terms, so that the Taylor expansion of the phase function in Eq. (22) is exact. Finally we get the same form with the extracted linear phase term like in Eq. (5):

V_{ℓ} (ρ, 0) = {\hat{V}}_{ℓ} (ρ, 0) exp [i (ε_{ℓ} + Δ_{ℓ} (ρ - ρ_{0, ℓ}, 0)] exp [i (ρ \cdot κ_{0, ℓ})]

For the decomposition in space domain the same basis function with smooth edges (see Eq. (19)) like in the case of the parabasal decomposition in spatial frequency domain is used.

5. Examples for non-paraxial field propagation

In the previous sections an efficient and rigorous semi-analytical SPW propagation operator for parabasal fields was derived. It was also shown how to decompose a general non-paraxial harmonic field into parabasal subfields using a PDT. Now the combination of both techniques is applied to propagate non-paraxial fields. In the following two examples are discussed.

5.1. Example I: convergent spherical wave

At first a non-paraxial convergent spherical wave with a spherical wave radius of R = −4 mm is propagated in vacuum by z = 3.8 mm, meaning that the target plane is slightly defocused. The diameter of the x-direction linearly polarized field is 1.28 mm at a wavelength of 532 nm. Figure 10 shows the initial spherical phase. A PDT in space domain is applied due to the field of second kind behavior of the spherical wave. The initial field is decomposed in 20 × 20 subfields to get parabasal subfields with a sufficiently small divergence. Figure 11 shows the residual phase of a single parabasal subfield after extracting the local linear phase. Next each single parabasal subfield is propagated 3.8 mm by the semi-analytical SPW operator. Figure 12 displays the residual amplitude and phase of the propagated subfield. Finally all propagated parabasal subfields are superimposed, including the extracted linear phase factors, to obtain the propagated non-paraxial defocused spherical wave. The result is given in Fig. 13.

Fig. 10 Initial phase of a convergent spherical wave which will be propagated by z = 3.8 mm. Please note that the artefacts in the phase distribution is a Moiré pattern due to the low resolution of the picture in comparison with the very fine sampling which is required for the phase.

Download Full Size | PDF

Fig. 11 Residual amplitude (left) and phase (right) of a parabasal subfield created by a PDT in space domain of a spherical wave into 20 × 20 subfields and extraction of the linear phase term.

Download Full Size | PDF

Fig. 12 Residual amplitude (left) and phase (right) of a propagated parabasal subfield calculated by the semi-analytical SPW operator.

Download Full Size | PDF

Fig. 13 Amplitude (left) and phase (right) of the defocused non-paraxial spherical wave propagated by a combination of PDT and semi-analytical SPW operator.

Download Full Size | PDF

Table 2 compares the deviation and the numerical effort for the propagation of the defocused spherical wave using different propagation techniques. Again the conventional SPW operator suffers from high numerical effort and the Fresnel and far field approximations result in large simulation error. However the combination of the PDT and the rigorous semi-analytical SPW operator leads to dramatically reduced numerical effort in terms of required RAM access as well as high accuracy. Also the complexity factor evaluated by Eq. (17) roughly yields an eight times reduced complexity of the algorithm compared to the conventional SPW operator.

Table 2. Accuracy and numerical effort for different propagation techniques to propagate the defocused spherical wave by 3.8 mm. In case of the PDT and semi-analytical SPW Operator the numerical effort in brackets gives the number of used parabasal subfields.

View Table | View all tables in this article

5.2. Example II: astigmatic Gaussian beam

To illustrate the PDT in spatial frequency domain a linearly polarized astigmatic Gaussian beam is propagated in vacuum by z = 0.5 mm. The waist radius in x-direction is 7 μm and in y-direction 10 μm. The wavelength is 532 nm. Additionally an astigmatic aberration is added by the corresponding Zernike coefficient of 60 in x and 120 in y-direction. Figure 14 shows the initial phase of the astigmatic Gaussian beam. Due to small wavefront features of the astigmatic Gaussian beam a PDT in spatial frequency in 15 × 15 subfields is used. Figure 15 shows the residual phase of a parabasal subfield after decomposition and extraction of the linear phase term. After that the semi-analytical SPW operator yields the propagated subfield of Fig. 16. Finally, after adding of all analytically stored linear phase terms this leads to the propagated non-paraxial astigmatic Gaussian beam of Fig. 17.

Fig. 14 Amplitude (left) and phase (right) of an astigmatic Gaussian beam with waist radius of 7 μm in x-direction and 10 μm in y-direction. Additional astigmatism is added by Zernike coefficients.

Download Full Size | PDF

Fig. 15 Residual amplitude (left) and phase (right) of a parabasal subfield created by a PDT in spatial frequency domain of an astigmatic Gaussian beam into 15 × 15 subfields and extraction of the local linear phase terms.

Download Full Size | PDF

Fig. 16 Residual amplitude (left) and phase (right) of a propagated parabasal subfield calculated by the semi-analytical SPW operator.

Download Full Size | PDF

Fig. 17 Amplitude of the non-paraxial astigmatic beam propagated by a combination of PDT and semi-analytical SPW operator.

Download Full Size | PDF

Table 3 compares the deviation and the numerical effort for the propagation of the astigmatic Gaussian beam using different propagation techniques. Again the conventional SPW operator suffers from high numerical effort and the Fresnel and far field approximations result in large simulation error. However the combination of the PDT and the rigorous semi-analytical SPW operator leads to dramatically reduced numerical effort as well as high accuracy. Also the complexity factor evaluated by Eq. 17 roughly yields to a reduced complexity of the algorithm by 10% compared to the conventional SPW operator.

Table 3. Accuracy and numerical effort for different propagation techniques to propagate the astigmatic Gaussian beam by 0.5 mm. In case of the PDT and semi-analytical SPW Operator the numerical effort in brackets gives the number of used parabasal subfields.

View Table | View all tables in this article

6. Conclusion

In this paper we have derived a rigorous semi-analytical SPW propagation operator. Basically the idea is an extraction of linear phase terms from the conventional SPW propagation integral kernel and its analytical handling due to the shift theorem of the Fourier transformation. The analytical handling of linear phase terms results in case of parabasal fields in a significant reduction of the required sampling points. It was shown by use of an example that this semi-analytical concept allows a rigorous propagation of parabasal fields with drastically reduced numerical effort.

Also it was proved that a combination of a parabasal decomposition technique with the semi-analytical SPW operator enables the rigorous propagation of general non-paraxial fields with reduced numerical effort. There are two different kinds of parabasal decomposition techniques:

in spatial frequency domain for fields with a rough wavefront
in space domain for fields with a smooth but strongly curved wavefront

Each decomposition technique was presented with the aid of an example.

A rough estimate of the complexity factor also showed the potential of the new semi-analytical propagation algorithm to reduce the required computational time solving the SPW-integral. In further work this advantage in computational speed has to be proven by performing a program runtime test.

References and links

1. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58(5–6), 449–466 (2011). [CrossRef]

2. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University Press, 1995).

3. A. Wuttig, M. Kanka, H. J. Kreuzer, and R. Riesenberg, “Packed domain Rayleigh-Sommerfeld wavefield propagation for large targets,” Opt. Express 18(26), 27036–27047 (2010). [CrossRef]

4. J. A. C. Veerman, J. J. Rusch, and H. P. Urbach, “Calculation of the Rayleigh–Sommerfeld diffraction integral by exact integration of the fast oscillating factor,” J. Opt. Soc. Am. A 22(4), 636–646 (2005). [CrossRef]

5. M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6(5), 786–805 (1989). [CrossRef]

6. J. Braat, P. Dirksen, and A. J. E. M. Janssen, “Assessment of an extended Nijboer–Zernike approach for the computation of optical point-spread functions,” J. Opt. Soc. Am. A 19(5), 858–870 (2002). [CrossRef]

7. P. Valtr and P. Pechac, “Domain decomposition algorithm for complex boundary modeling using the Fourier split-step parabolic equation,” IEEE Trans. Antennas Propag. 6, 152–155 (2007).

8. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).

9. LightTrans GmbH, LightTrans VirtualLab Advanced^TM, www.lighttrans.com (2012).

10. E. O. Brigham, The Fast Fourier Transform and its Applications (Prentice Hall, 1988).

Operator	Deviation vs. SPW	numerical effort (sampling points)

SPW	0%	9156 × 9156
Fresnel	107%	266 × 266
Far field	46%	533 × 533
semi-analytical SPW	0.02%	279 × 279

Operator	Deviation vs. SPW	numerical effort (sampling points)

SPW	0%	8585 × 8585
Fresnel	159%	4658 × 4658
Far field	200%	8419 × 8419
PDT + semi-analytical SPW	0.02%	98 × 98 × (20 × 20)

Operator	Deviation vs. SPW	numerical effort (sampling points)

SPW	0%	10869 × 10869
Fresnel	124%	3843 × 3843
Far field	25.3%	3843 × 3843
PDT + semi-analytical SPW	0.08%	240 × 240 × (15 × 15)

Operator	Deviation vs. SPW	numerical effort (sampling points)

SPW	0%	9156 × 9156
Fresnel	107%	266 × 266
Far field	46%	533 × 533
semi-analytical SPW	0.02%	279 × 279

Operator	Deviation vs. SPW	numerical effort (sampling points)

SPW	0%	8585 × 8585
Fresnel	159%	4658 × 4658
Far field	200%	8419 × 8419
PDT + semi-analytical SPW	0.02%	98 × 98 × (20 × 20)

Parabasal field decomposition and its application to non-paraxial propagation

Abstract

1. Introduction

2. Fundamentals of parabasal fields

3. Semi-analytical SPW propagation operator

4. Parabasal field decomposition technique

5. Examples for non-paraxial field propagation

5.1. Example I: convergent spherical wave

5.2. Example II: astigmatic Gaussian beam

6. Conclusion

References and links

Cited By

Figures (17)

Tables (3)

Equations (25)

Optics Express