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

## 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)]=12π∬−∞∞Vℓ(ρ,0)e−iκ⋅ρdρ.$
Here = 1...6 denotes the 3 electric and the 3 magnetic field components of the harmonic field and κ = (kx, ky)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)eikz(κ)z]$
with the spherical phase function $kz(κ)=k2−kx2−ky2$, 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 = k0n with refractive index n and vacuum wavenumber k0. 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 kz(κ). 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.

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 kz(κ). 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 VirtualLabTM [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

$k=ks^0=k(s0x,s0y,s0z)T=k(sinθ0cosϕ0,sinθ0sinϕ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=(k0x,k0y)T=k(s0x,s0y)T$
with the wavenumber in the homogeneous medium k = k0n. 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 𝔕k on the kxky-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 ŝ0. For the sake of simplicity the two dimensional case where s0y = 0 is chosen. In general a parabasal field can be written as

$Vℓ(ρ,0)=Vℓ′(ρ,0)eiκ0,ℓ⋅ρ$
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 ŝ0, 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 κ0.

## 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 kz(κ) independently for each component around κ0, by

$kz,ℓ(κ)=γ0,ℓ+γ1,ℓ⋅κ+γℓ(κ−κ0,ℓ)$
with the constant term
$γ0,ℓ=k0z,ℓ+1k0z,ℓ(k0x,ℓ2+k0y,ℓ2)−k0x,ℓ2k0y,ℓ2k0z,ℓ3,$
and the linear term
$γ1,ℓ=[k0x,ℓk0z,ℓ(k0y,ℓ2k0z,ℓ2−1),k0y,ℓk0z,ℓ(k0x,ℓ2k0z,ℓ2−1)]T$
with
$k0z,ℓ=k2−k0x,ℓ2−k0y,ℓ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 κ′ = (kx, ky)T with center at κ0,. That results in
$γℓ(κ′)=kz,ℓ(κ′+κ0,ℓ)−γ0,ℓ−γ1,ℓ⋅(κ′+κ0,ℓ)$
with
$kz,ℓ(κ′+κ0,ℓ)=k2−(kx′+k0x,ℓ)2−(ky′+k0y,ℓ)2.$
Please note that the decomposition of kz(κ) 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
$Vℓ(ρ,z)=Vℓ′(ρ+zγ1,ℓ,z)exp[i(γ0,ℓz+γ1,ℓ⋅κ0,ℓz+ρ⋅κ0,ℓ)]$
with the propagated shifted residual field
$Vℓ′(ρ,z)=ℱ−1[Aℓ′(κ,0)eiγℓ(κ)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 kz(κ) in Eq. (2). As shown in Fig. 3 the slope of γ(κκ0,) for κ around κ0, is much flatter than the slope of kz(κ). 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 kz(κ) (blue line) and its decomposition (magenta and red line) according to Eq. (8) along the diagonal direction $(kx,ky)T=kρ2(1,1)T$. Left side: Normalized spherical phase function kz(κ) and linear phase function γ0, + γ1, · κ around $k0ρk=0.424$. Right side: Higher order phase term γ(κκ0,) for $k0ρk=0.424$ calculated by Eq. (12).

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.

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:=∑x,y|Vℓ,operator(x,y,z)−Vℓ,reference(x,y,z)|2∑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:

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

$η:=Complexity of SPW algorithmComplexity 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(N2 log2 N)) and the complexity of the multiplication (O(N2)) 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.

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.

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

$Aℓ(κ,0)=∑i,j=1NΦ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 = (ki,x, kj,y)T located e.g. on an equidistant rectangular grid. Also the property $∑i,j=1NΦ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(kx, ky) < 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).

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

$Φ(kx)=0.5[sin(πa(kx−a2))+1]for 0
with a, b ∈ ℝ and 0 < a < b. This definition can be generalized to the two-dimensional case using a separation Ansatz
$Φi,j(κ)=Φ(kx−ki,x)Φ(ky−kj,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).

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:

1. 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.
2. 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.

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)=V^ℓ(ρ,0)eiφℓ(ρ,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)|(ρ≈ρ0,ℓ)=εℓ+κ0,ℓ⋅ρ+Δℓ(ρ,0)$
with
$εℓ=φℓ(ρ0,ℓ,0)+∂2φℓ∂x∂y|ρ0,ℓ⋅(x0,ℓy0,ℓ)−∂φℓ∂x|ρ0,ℓ⋅x0,ℓ−∂φℓ∂y|ρ0,ℓ⋅y0,ℓ$
and
$κ0,ℓ=(k0x,ℓ,k0y,ℓ)T=(∂φℓ∂x|ρ0,ℓ−∂2φℓ∂x∂y|ρ0,ℓ⋅y0,ℓ,∂φℓ∂y|ρ0,ℓ−∂2φℓ∂x∂y|ρ0,ℓ⋅x0,ℓ)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)=V^ℓ(ρ,0)exp[i(εℓ+Δℓ(ρ−ρ0,ℓ,0)]exp[i(ρ⋅κ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.

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.

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

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.

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.

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

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.

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

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

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.

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

1. in spatial frequency domain for fields with a rough wavefront
2. 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 AdvancedTM, www.lighttrans.com (2012).

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

### References

• View by:
• |
• |
• |

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 AdvancedTM, www.lighttrans.com (2012).
10. E. O. Brigham, The Fast Fourier Transform and its Applications (Prentice Hall, 1988).

#### 2011 (1)

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

#### 2007 (1)

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

#### Brigham, E. O.

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

#### Goodman, W.

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

#### Kuhn, M.

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

#### Mandel, L.

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

#### Pechac, P.

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

#### Valtr, P.

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

#### Wolf, E.

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

#### Wyrowski, F.

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

#### IEEE Trans. Antennas Propag. (1)

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

#### J. Mod. Opt. (1)

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

#### Other (4)

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

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

LightTrans GmbH, LightTrans VirtualLab AdvancedTM, www.lighttrans.com (2012).

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

### Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

### Figures (17)

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.

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 ŝ0, 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 κ0.

Fig. 3

2D line evaluation of the spherical phase kz(κ) (blue line) and its decomposition (magenta and red line) according to Eq. (8) along the diagonal direction $( k x , k y ) T = k ρ 2 ( 1 , 1 ) T$. Left side: Normalized spherical phase function kz(κ) and linear phase function γ0, + γ1, · κ around $k 0 ρ k = 0.424$. Right side: Higher order phase term γ(κκ0,) for $k 0 ρ k = 0.424$ calculated by Eq. (12).

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.

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.

Fig. 6

Decomposition of a first order Laguerre Gaussian beam into 2 × 2 parabasal subfields in the spatial frequency 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).

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

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.

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.

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.

Fig. 12

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

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.

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.

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.

Fig. 16

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

Fig. 17

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

### Tables (3)

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

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.

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.

### Equations (25)

$A ℓ ( κ , 0 ) = ℱ [ V ℓ ( ρ , 0 ) ] = 1 2 π ∬ − ∞ ∞ V ℓ ( ρ , 0 ) e − i κ ⋅ ρ d ρ .$
$V ℓ ( ρ , z ) = ℱ − 1 [ A ℓ ( κ , 0 ) e i k z ( κ ) z ]$
$k = k 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$
$κ 0 = ( k 0 x , k 0 y ) T = k ( s 0 x , s 0 y ) T$
$V ℓ ( ρ , 0 ) = V ℓ ′ ( ρ , 0 ) e i κ 0 , ℓ ⋅ ρ$
$A ℓ ( κ , 0 ) = A ℓ ′ ( κ − κ 0 , ℓ , 0 )$
$V ℓ ′ ( ρ , 0 ) = ℱ − 1 [ A ℓ ′ ( κ , 0 ) ] .$
$k z , ℓ ( κ ) = γ 0 , ℓ + γ 1 , ℓ ⋅ κ + γ ℓ ( κ − κ 0 , ℓ )$
$γ 0 , ℓ = k 0 z , ℓ + 1 k 0 z , ℓ ( k 0 x , ℓ 2 + k 0 y , ℓ 2 ) − k 0 x , ℓ 2 k 0 y , ℓ 2 k 0 z , ℓ 3 ,$
$γ 1 , ℓ = [ k 0 x , ℓ k 0 z , ℓ ( k 0 y , ℓ 2 k 0 z , ℓ 2 − 1 ) , k 0 y , ℓ k 0 z , ℓ ( k 0 x , ℓ 2 k 0 z , ℓ 2 − 1 ) ] T$
$k 0 z , ℓ = k 2 − k 0 x , ℓ 2 − k 0 y , ℓ 2 .$
$γ ℓ ( κ ′ ) = k z , ℓ ( κ ′ + κ 0 , ℓ ) − γ 0 , ℓ − γ 1 , ℓ ⋅ ( κ ′ + κ 0 , ℓ )$
$k z , ℓ ( κ ′ + κ 0 , ℓ ) = k 2 − ( k x ′ + k 0 x , ℓ ) 2 − ( k y ′ + k 0 y , ℓ ) 2 .$
$V ℓ ( ρ , z ) = V ℓ ′ ( ρ + z γ 1 , ℓ , z ) exp [ i ( γ 0 , ℓ z + γ 1 , ℓ ⋅ κ 0 , ℓ z + ρ ⋅ κ 0 , ℓ ) ]$
$V ℓ ′ ( ρ , z ) = ℱ − 1 [ A ℓ ′ ( κ , 0 ) e i γ ℓ ( κ ) z ]$
$d : = ∑ x , y | V ℓ , operator ( x , y , z ) − V ℓ , reference ( x , y , z ) | 2 ∑ x , y | V ℓ , reference ( x , y , z ) | 2$
$η : = Complexity of SPW algorithm Complexity of semi-analytical SPW algorithm$
$A ℓ ( κ , 0 ) = ∑ i , j = 1 N Φ i , j ( κ ) A ℓ ( κ , 0 )$
$Φ ( k x ) = 0.5 [ sin ( π a ( k x − a 2 ) ) + 1 ] for 0 < k x ≤ a , Φ ( k x ) = 1 for a < k x < b − a , Φ ( k x ) = 1 − 0.5 [ sin ( π a ( k x − 2 b − a 2 ) ) + 1 ] for b − a ≤ k x < b , and Φ ( k x ) = 0 anywhere else ,$
$Φ i , j ( κ ) = Φ ( k x − k i , x ) Φ ( k y − k j , y )$
$V ℓ ( ρ , 0 ) = V ^ ℓ ( ρ , 0 ) e i φ ℓ ( ρ , 0 )$
$φ ( ρ , 0 ) | ( ρ ≈ ρ 0 , ℓ ) = ε ℓ + κ 0 , ℓ ⋅ ρ + Δ ℓ ( ρ , 0 )$
$ε ℓ = φ ℓ ( ρ 0 , ℓ , 0 ) + ∂ 2 φ ℓ ∂ x ∂ y | ρ 0 , ℓ ⋅ ( x 0 , ℓ y 0 , ℓ ) − ∂ φ ℓ ∂ x | ρ 0 , ℓ ⋅ x 0 , ℓ − ∂ φ ℓ ∂ y | ρ 0 , ℓ ⋅ y 0 , ℓ$
$κ 0 , ℓ = ( k 0 x , ℓ , k 0 y , ℓ ) T = ( ∂ φ ℓ ∂ x | ρ 0 , ℓ − ∂ 2 φ ℓ ∂ x ∂ y | ρ 0 , ℓ ⋅ y 0 , ℓ , ∂ φ ℓ ∂ y | ρ 0 , ℓ − ∂ 2 φ ℓ ∂ x ∂ y | ρ 0 , ℓ ⋅ x 0 , ℓ ) T .$
$V ℓ ( ρ , 0 ) = V ^ ℓ ( ρ , 0 ) exp [ i ( ε ℓ + Δ ℓ ( ρ − ρ 0 , ℓ , 0 ) ] exp [ i ( ρ ⋅ κ 0 , ℓ ) ]$