## Abstract

The fast and accurate propagation of general optical fields in free space is still a challenging task. Most of the standard algorithms are either fast or accurate. In the paper we introduce a new algorithm for the fast propagation of non-paraxial vectorial optical fields without further physical approximations. The method is based on decomposing highly divergent (non-paraxial) fields into subfields with small divergence. These subfields can then be propagated by a new semi-analytical spectrum of plane waves (SPW) operator using fast Fourier transformations. In the target plane, all propagated subfields are added coherently. Compared to the standard SPW operator, the numerical effort is reduced drastically due to the analytical handling of linear phase terms arising after the decomposition of the fields. Numerical results are presented for two examples demonstrating the efficiency and the accuracy of the new method.

© 2012 OSA

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

*ℓ*= 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*) in the parallel target plane. Please note that we have used the wavenumber in the homogeneous medium

*k*=

*k*

_{0}

*n*with refractive index

*n*and vacuum wavenumber

*k*

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

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

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

*κ*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 *κ*_{0}. It means that a single parabasal harmonic field component possesses a low divergence and propagates with mean direction along

*θ*and

*ϕ*with unity radius. If we restrict to non-lossy dielectrics, the direction vector defines the central spatial frequency vector

*k*=

*k*

_{0}

*n*. In the special case of

*κ*_{0}= (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

**on the**

*k**k*−

_{x}*k*-plane. For simplification we will just discuss non-absorbing dielectrics in this paper.

_{y}Figure 2 illustrates the lower limit of a parabasal harmonic field with zero divergence, which is a plane wave pointing into the direction *ŝ*_{0}. 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) 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**

*ρ*## 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

*γ*(

_{ℓ}**−**

*κ*

*κ*_{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}**) in Eq. (8) represents a Taylor expansion around**

*κ*

*κ*_{0}with the inclusion of all higher order terms in

*γ*(

_{ℓ}**−**

*κ*

*κ*_{0}). Plugging Eq. (8) into Eq. (2) of the SPW operator and applying the shift theorem of the Fourier transformation leads to

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

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

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

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

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

*N*

^{2}log

_{2}

*N*)) and the complexity of the multiplication (O(

*N*

^{2})) 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.

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

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

_{i}_{,j}(

**) is the basis function with local carrier and Φ**

*κ**(*

_{i,j}**) := Φ(**

*κ***−**

*κ**) with a plenty of supporting points*

**κ**_{i,j}*= (*

**κ**_{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}{\mathrm{\Phi}}_{i,j}(\mathit{\kappa})=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*) < 1 and Φ(

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

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

*a*,

*b*∈ ℝ and 0 <

*a*<

*b*. This definition can be generalized to the two-dimensional case using a separation Ansatz

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

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.

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) 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) 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):**

*ρ*## 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.

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.

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

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.

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

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