Analytic expressions and approximations for the on-axis , aberration-free Rayleigh and Debye integral in the case of focusing fields on a circular aperture

R. M. Aarts Philips Research Europe, HTC-36, NL 5656 AE Eindhoven, The Netherlands. Faculty of Electrical Engineering, Eindhoven University of Technology, NL-5600 MB Eindhoven, The Netherlands. J. J. M. Braat Optics Research Group, Faculty of Applied Sciences, Technical University Delft, Lorentzweg 1, NL-2628 CJ Delft, The Netherlands. P. Dirksen Philips Research Europe, HTC-04, NL-5656 AE Eindhoven, The Netherlands.


INTRODUCTION
The propagation of the optical or acoustical disturbance from the aperture or pupil towards the focal region admits various treatments. In this paper we consider the disturbance to be a scalar quantity, typically the pressure in the acoustical domain or the amplitude of the field in the optical domain. More refined propagation models including vector effects (propagation of density waves in solids, propagation of optical fields with a high degree of focusing) are outside the scope of this paper. Within the scalar approximation, we are left with two, or even three, possible representations, due to Rayleigh [1] and Kirchhoff [2]. In this paper we will use the Rayleigh integral that connects the pressure distribution on the aperture to the pressure distribution in the near and far field [3], the socalled Rayleigh-I integral [4]. We will also consider the special approximated diffraction integral that was proposed by Debye for focused fields [5]. For this integral to be a good approximation, the number of Fresnel zones in the aperture should be substantially larger than unity [6]- [8]. In the optical domain, where also the lateral dimension of the diffracting aperture is mostly many wave lengths large (typically in excess of 10 4 wave lengths), this condition is generally met. In that case, an approximated version of the Rayleigh integral leads to the Debye integral. This integral representation of the focused field was further developed for high-numerical-aperture optical diffraction problems including the vector character of the electromagnetic optical field [9]- [11] and a nonuniform amplitude mapping from the entrance to the exit pupil of a thick lens system. Possible aberration and transmission defects of the focusing wave field have been treated in [12]- [17]. It is important to note that the Debye integral also neglects the nearfield diffraction term that is present in the Rayleigh integral. In the optical domain, the distance from the aperture to the focal region is many wave lengths large and the near-field effects do not need to be considered. This condition even holds for microlenses [18]. In the acoustical domain, the near field can be important. To check the influence of the near-field term, we will also treat in this paper the so-called incomplete Rayleigh integral where the near-field term in the integrand has been omitted.
Regarding the boundary conditions at the edge of the aperture, throughout this paper we will adopt the so-called hard boundary conditions. These allow a discontinuity in the pressure distribution or in the optical field (Kirchhoff boundary condition). Although physically unrealistic, these boundary conditions are acceptable if the lateral size of the aperture is larger than, say, ten to twenty wave lengths. Although these boundary conditions are often applied for even smaller diameters, it should be born in mind that the real physical result might seriously deviate from the one in the Kirchhoff approx-imation. In general, one could say that in acoustic diffraction problems, the Kirchhoff conditions are not always applicable, but in ultrasound applications and certainly in the optical domain one is in a safe region for applying these hard Kirchhoff boundary conditions. An extensive literature is available on the subjects described above. The acoustic field of a focusing radiator was studied in [19] and a deviation of the location of highest acoustic power from the geometric focus was demonstrated, both theoretically and experimentally. Numerical calculations and measurements of the complex field in the focal region of a smallaperture microlens have been carried out by Farnell [20]. An analytic result for the axial field component behind an aperture illuminated by a converging spherical wave was given (without derivation) by Osterberg and Smith [21] in a seminal paper that also draws interesting physical conclusions from the nature of the analytic solution. The special effects arising at small aperture focusing, viz. the asymmetry of the axial intensity with respect to the geometrical focus, has been been treated in [22]- [26]. Special attention has also been paid to the transition from low to high aperture when the quadratic pathlength approximation starts to produce incorrect results. Extensions and improvements of the Fresnel appoximation have been proposed to adequately address the higher aperture case [27]- [29]. The asymmetry around focus, first remarked at low aperture, has also been analyzed at high aperture in the framework of high-resolution three-dimensional microscopy. An analysis of the focal shift as a function of Fresnel number and aperture is presented in references [30]- [32]. In reference [31], an analytic result for the on-axis field as predicted by the Kirchhoff diffraction integral is given, using the initial Osterberg and Smith result. The near field, of special interest in the acoustic domain, has been studied in [33]- [39]. From the numerical point of view, the computation of the near-field part of the diffraction integral is the most challenging because of its strongly oscillatory behaviour. The Debye integral, mainly used in optical diffraction problems, has been analyzed with respect to its domain of validity in [25] and [40].
In this paper, with respect to the existing literature, we present some new explicit analytic and semi-analytic results for the axial fields represented by the various diffraction integrals. These analytic results are then evaluated and compared with standard numerical calculation results. We start by presenting a proof of the Osterberg and Smith analytic result for the Rayleigh integral. Although the original diffraction problem presented by Osterberg and Smith applied to an integration over a sphere and ours to integration over a plane, the proof can also be directly applied to the integration over a spherical cap. Using the method we applied to obtain the Osterberg-Smith result, we address the incomplete Rayleigh integral and we develop an approximating analytic expression that will be compared with the result following from a numerical evaluation of the incomplete Rayleigh integral. We will especially consider the ability of the various integral expressions and analytic results, exact and approximated, to produce the correct value of the diffracted field relatively close to the diffracting aperture. In all cases, the analysis of the diffracted and focused field is limited to the axis perpendicular to the aperture and going through its center. We also suppose a uniformly focused wave field in the aperture and do not include amplitude or phase deviations from this ideal profile in the aperture. This special case can be treated analytically to a large extent, and thus provides a well-understood yardstick for the validity of the various approximations. From there one can extrapolate towards assessment of off-axis performance in the presence of aberrations.
The paper has been organized as follows. In Section 2, we present the general expressions for the complete and incomplete Rayleigh integral and for the Debye integral. In Section 3 we develop the special form of these integrals in the case of circular symmetry and on-axis field evaluation. In Section 4 we present the new derivation of the analytic expressions for the on-axis field after the diffracting aperture that are then used to obtain accurate and fast evaluations of the on-axis fields. In Section 5 we concentrate on the special cases that arise when calculating the field in the aperture itself (extreme 'near field') and close to the focal point. Finally, in Section 6, we present numerical examples illustrating the domain of applicability of the various integral representations.

COMPLETE AND INCOMPLETE RAYLEIGH INTEGRAL AND DEBYE INTEGRAL
We consider in this section fields on a circular aperture A of radius a and the associated integral expressions of Rayleigh and Debye for the (scalar) field due to focusing; we restrict attention to the field values on the axis through the center of the aperture and perpendicular to it. That is, we consider in the aperture plane z = z a field of the form where k is the wave number of the focusing wave, on the aperture to the on-axis field point F(x f , y f , z f ) with x f = y f = 0, and E 0 is an amplitude factor that vanishes outside A. The subscript f refers to 'focal'. The direction-dependent amplitude factor E 0 (x , y ) has the dimension of amplitude times meter and is an invariant quantity along a straight line joining Q(x , y , z ) in the field and the focal point F. We refer to Figure 1 for the configuration in the case of a circular aperture and for definitions and notations.
It follows from the Huygens-Fresnel principle that the field E f (x, y, z; z ) in the plane with axial coordinate z is basically given by each of two integrals of the Rayleigh type. In this paper, we choose the Rayleigh integral (usually called the Rayleigh-I integral) that uses the field in the aperture in the integrand and not the z-derivative of the field. This allows us to have a direct relationship between comparable quantities in the aperture domain and in the near and far field. We thus  write 1/2 the distance between the field point P(x, y, z) in the observation plane and the point Q(x , y , z ) on the aperture A. Eq. (2) takes a simpler form in the case that P is on the axis (x = y = 0) and the amplitude factor E 0 in Eq. (1) for E A equals unity throughout the aperture, pertaining to a spherical wave with uniform amplitude. We shall turn to this case later. Carrying out the Now it is customary in optics, where in many cases kR QP >> 1, to ignore the −1 in the term (ikR QP − 1) occurring in the integrand in Eq. (3) so that one arrives at, what is sometimes called, the incomplete Rayleigh integral It is interesting to see how close the field point P should get to the aperture and how small k should get so that the incomplete expression (4) ceases to be an accurate approximation to the complete expression (2) or (3). We shall consider this question for the case that P is on the axis and E 0 is uniform on A.
A second approximation, essentially due to Debye, for the case of a focused field E A as in Eq. (1), is obtained by adopting a spectral approach. The angular spectrum of the field in the aperture is given bỹ with nonnegative square roots at the right-hand side of Eq. (7) in either case. In the case of the purely imaginary value for k z , this choice assures an exponentially decreasing field in the propagation direction. The integral in Eq. (6) extends over all real values of (k x , k y ). In order to make the integration range into a bounded set,Ẽ(z ; k x , k y ) is approximated bỹ inside Ω andẼ(z ; k x , k y ) = 0 outside Ω. Here E A has been assumed as in Eq.
and Ω denotes the solid angle that the aperture A subtends at F. The quantity E 0 , with dimension of field strength times meter, is commonly denoted as the 'ray strength' because it is an invariant quantity along a fixed propagation direction in the focusing field distribution. E 0 is obtained by multiplying the field strength in a point on a spherical surface, having its center at the focus, with the radius of curvature of that surface. The value of E 0 in the specific case of a telecentric system with the exit pupil at infinity or of a telescopic system with both pupils at infinity is treated in [17].
Eq. (8) results upon asymptotic expansion ofẼ in Eq. (5) for the case that E 0 is uniform on the aperture by using the stationary phase method [8] in which only the contribution of the dominant stationary point is retained and the rim contribution to the integral is neglected.

COMPLETE AND INCOMPLETE RAYLEIGH INTEGRAL AND DEBYE INTEGRAL FOR THE ON-AXIS, ABERRATION-FREE CASE
In this section we develop the special form of the three integrals to be considered for the case of circular symmetry. The circular symmetry is present due to the fact that we choose our observation point on the axis of symmetry of the aperture; moreover, we limit ourselves to a perfect focusing field without any amplitude variation or phase aberration. The constant amplitude approximation is generally not allowed in highnumerical-aperture large-field imaging with optical lenses. An amplitude roll-off proportional to k 1/2 z is found in this case [10], but for the relatively low values of the aperture s 0 that we consider in this paper, the amplitude can be assumed to be uniform. We then take a unit radius of curvature for the spherical surface on which E 0 is measured and assume that the ray strength E 0 (x , y ) has unit value for (x , y ) ∈ A and is zero for (x , y ) outside A. We obtain for the complete Rayleigh integral in Eq. (2) the expression (10) The on-axis value of E f results upon taking x = y = 0. Then the integrand becomes radially symmetric and it follows that where (12) again, see Figure 1. Next, we introduce normalized variables ρ, S and T according to with T being the normalized distance from the aperture plane to the focal point F of the converging beam. It then follows that Note that the quantity on the right-hand side of Eq. (15) between [ ] depends on z through S = (z − z )/a, see Eq. (13).
In a similar fashion, there follows for the incomplete Rayleigh integral in Eq. (4) the on-axis, aberration-free expression We finally consider the Debye integral in Eq. (9) on axis for the aberration-free case, with with Ω, for the circularly symmetric case, given by a conical solid angle that is delimited by the circular aperture A. Thus, the integration range is here With the upper option in the definition of k z in Eq. (7), it then follows that Using polar coordinates (k x , k y ) = kκ(cos θ, sin θ), where 0 ≤ κ ≤ 1, 0 ≤ θ ≤ 2π, we then find where we have substituted τ = (1 − κ 2 ) 1/2 ∈ 1 − s 2 0 , 1 . Finally, in terms of the variables introduced in Eq. (13), we have and there results

ANALYTIC AND SEMI-ANALYTIC EXPRESSIONS FOR THE COMPLETE AND INCOMPLETE RAYLEIGH INTEGRAL AND DEBYE INTEGRAL ON AXIS
In this section we show analytic or semi-analytic expressions for the three diffraction integrals that were derived in the previous section. The analytic expressions allow a fast and accurate calculation of the near and far field everywhere behind the diffracting aperture on axis. Without an explicit derivation, the analytic result for the complete Rayleigh integral was already given in [21]. We first present such a derivation for the complete Rayleigh-I integral and then use the same procedure to produce analytic or semi-analytic results for the other integrals of interest.

Derivation of the analytic expression for the complete Rayleigh integral
We have To prove this result, we start from Eq. (15) and substitute in the integral Noting that and that so that ρ dρ we get Now differentiate Eq. (29) with respect to z, with y(0) and y(1) depending on z through Eq. (26) and S = (z − z )/a, to get This then gives Eq. (23) as required.

Integral expression for the incomplete Rayleigh integral
In this subsection we derive an adapted and simplified integral expression for the incomplete Rayleigh integral E f ,inc . In contrast with the complete Rayleigh integral E f , it was not possible to develop an exact analytic expression for E f ,inc . The analysis could not be continued beyond an expression in terms of sine-and cosine-integrals. With a further approximation to the value of the integrand in E f ,inc , an analytic expression becomes feasible and, further on in Section 6, we will study the influence of this approximation of E f ,inc on the field obtained close by and further away from the diffracting aperture.
We shall now show that To prove this result, we start from Eq. (16) and use the substitution (24). It follows now that with ρ(y) the inverse function of y(ρ) in Eq. (24). It follows from Eq. (24) that i.e., that 2y S 2 + ρ 2 (y) = y 2 − T 2 + S 2 . This completes the proof of Eq. (31).

The integral expression (31) can be brought into a form involving the sine and cosine integrals
Si(z) = z 0 sin t t dt, with γ = 0.5772 . . . (Euler's constant), see Section 5.2 of [41]. Indeed, letting and writing one readily finds However, this form is not very convenient when one is interested in simple and effective approximations of E f ,inc . Furthermore, the integration limits ka(V ± W), ka(U ± W) can become very large in modulus and are complex in the case that S > T.
The following observation can be made. The integral expression in Eq. (31) for E f ,inc is such that in many cases the function (y 2 − T 2 + S 2 ) −1 can be considered as being nearly constant. Letting U and V as in Eq. (35), we have displayed in Table 1 the ratio of the extreme values of (y 2 − T 2 + S 2 ) −1 on the integration range as a function of S ∈ [0, 5T], where T is fixed at the value 8. Note that the right-hand side of Eq. (38) equals which is a ratio of the derivative and a differential quotient of the very smoothly behaved function f . The table shows that somewhat beyond the value S = 1 the function (y 2 − T 2 + S 2 ) −1 can be regarded as almost constant. Hence, accurate approximations to E f ,inc can be found in the form whereŷ is a number chosen between the integration limits S − T and √ S 2 + 1 − √ T 2 + 1.

Analytic expression for the Debye integral
We have immediately from Eq. (20) that Recall that, also see Figure 1, Thus there results Note that the choiceŷ = S − T in Eq. (40) gives for the right-hand side and this comes quite close to the righthand side of Eq. (44) since for some X between S and T, see Eqs. (38) and (39).
We now have at our disposal exact analytic expressions for the complete Rayleigh integral E f , see Eq. (23), and for the Debye integral E D , see Eq. (44). Note that the expression for the Debye integral is even with respect to the quantity S − T which leads to (conjugate) symmetry with respect to the geometrical focus at the position S = T. For the incomplete Rayleigh integral, we either use the basic expression in Eq. (16) or we use the semi-analytic expression for E f ,inc in Eq. (31), that can be further evaluated by numerical integration. These expressions will be further used in Section 6 of this paper when we discuss various settings that are characteristic for acoustic or optical focused wave fields.

Value in focus ( z = z f )
We have at z = z f , so that S = T,

Value and behaviour at the aperture ( z = z )
We have at z = z , so that S = 0,
The analysis of E f ,inc at z = z is more awkward. For instance, the limiting value and behaviour of the right-hand side of Eq. (40) as z → z (S → 0) depends on the choice ofŷ. The choiceŷ = S − T and S → 0 yields the non-zero limit On the other hand, when we chooseŷ = , the right-hand side of Eq. (40) tends to 0 as S → 0. More precisely, we have the limiting behaviour as S → 0, and the factor in front of S does depend on A.
The behaviour of the exact expression for E f ,inc as z → z differs from the approximations just given. This is so since E f and E f ,inc differ by a term that becomes only non-negligible when kaS is of the order of unity or smaller while the approximations (51) and (52) are valid in a much larger range. To better understand this, we note that from Eq. (16) one has while from Eq. (15) (with S = (z − z )/a) where, for ease of notation, we have omitted the coordinate dependence of the expressions for E f and E f ,inc . Thus, the difference between the two is the term on the second line of Eq. (54). The integrals in Eq. (54) are of a very similar nature, but the factors by which they get multiplied are vastly different in order of magnitude. Therefore, the second term in Eq. (54) cannot be neglected anymore only when the additional factor 1/(S 2 + ρ 2 ) 1/2 in the second integral becomes of the order ka on a substantial part of the integration range, i.e., roughly when kaS ≤ 1.
Now E f tends to a finite limit (1/aT) exp{−ikaT} = 0 as S → 0 while E f ,inc → 0 as S → 0. More precisely, we have as an approximation from Eq. (37) when 0 < kaS ≤ 1 in which C is a constant of order unity that we shall ignore below. The quantity k ln(kaS)/T by which S is being multiplied, becomes very large compared to the limit value 1/aT of |E f |. Therefore, we see an extremely steep decay to zero of |E f ,inc | when kaS drops below 1.

NUMERICAL EXAMPLES AND DISCUSSION OF RESULTS
In this section we discuss some numerical results that have been obtained with the aid of our analytic expressions to illustrate the accuracy of the analytic results and the range of applicability of the approximations involved. We have the following expressions • Complete Rayleigh integral E f , with its on-axis analytic solution according to Eq. (23), • Incomplete Rayleigh integral E f ,inc with the near-field term omitted, numerical integration using Eq. (16); for this we use a low-order method using an adaptive recursive Simpson rule ('quad' implemented in MATLAB), • An approximated analytic expression of the incomplete Rayleigh integral E f ,inc according to Eq. (40) with the liberty to substitute an arbitrary value forŷ between the integration limits S − T and • An analytic result Eq. (44) for the Debye integral E D , an approximation of the Rayleigh integral for large values of the wave number k.
The geometry of the diffraction problem has been chosen such that, with an aperture diameter 2a = 37 mm and a radius of curvature of the focused wave front R = 71 mm, the numerical aperture of the focusing field equals 0.2606. The wave number of the radiation field has been chosen k = 12736 m −1 , corresponding to a wave length λ of approximately 0.5 mm. These values with typically 100 wave lengths fitting into the aperture are common for acoustic diffraction problems and they also occur in the optical domain when dealing with micro-lenses.
In Figure 2 we have plotted the modulus of the axial field pertaining to the complete Rayleigh integral (blue solid curve), using the analytic expression of Eq. (23). The incomplete Rayleigh integral E f ,inc (dotted green curve), approximated according to Eq. (40) with the choiceŷ = √ S 2 + 1 − √ T 2 + 1, closely follows the exact result in the focal region (the normalized distance T towards the geometrical focus corresponds to a value S = 3.70527). Close to the aperture, for S < 1, substantial deviations are observed. These are mainly due to the approximation that was used in deriving Eq. (40). The third curve, the Debye integral E D (dashed red), shows significant deviations from the exact result E f . The approximation involved in deriving the Debye integral induces symmetry with respect to the geometrical focus. It can be seen in the figure, that for the relatively small k-value we used, the Debye approximation is inadequate to reproduce the correct positions of the axial field zeros, the shift of the maximum field value with respect to the geometrical focus and the inherent asymmetry of the diffracted field with respect to the optimum focus, see also [17], [25], [31] for an analysis of this phenomenon and further examples.
In Figure 3 we highlight the difference between the complete and the incomplete Rayleigh integrals, E f and E f ,inc , respectively. In [39], the difference between E f and E f ,inc has also been studied, in this case for plane wave illumination. A numerical evaluation of the incomplete Rayleigh integral was compared with the analytic solution for the complete integral. In our case, for converging wave illumination, E f (solid blue curve) has been obtained using the analytic result of Eq. (23), E f ,inc (dotted magenta curve) has been calculated by numerical integration of the expression in Eq. (16).
It is seen from the figure that the differences between the complete and incomplete integral are negligible as soon as the observation point has moved out of the near-field region. Very close to the aperture, typically for |z − z | ≤ λ, we observe an important difference. The incomplete integral approaches zero field value while the complete integral correctly reproduces the field at the aperture itself. The steep drop to zero of E f ,inc agrees with the observation made at the end of Subsection 5.2. Note that this behaviour for the converging wave illumination is also visible in Figure 9(a) of [39] for plane wave illumination of a circular aperture. The insert in Figure 3 illustrates this phenomenon that in our example occurs in the axial range S ≤ 0.03.
In Figure 4 we compare the numerically obtained value of E f ,inc and the analytic expression of Eq. (40) that was obtained by introducing an approximation in the integrand of the ex- pression for E f ,inc . A free parameter in obtaining the analytic solution of Eq. (40) isŷ, that should be chosen between the integration limits S − T and √ S 2 + 1 − √ T 2 + 1. From the figure we deduce that the approximated results are reliable for values of S larger than typically 1.5. For smaller values of S, the choice ofŷ strongly influences the value of the diffracted field and strong oscillations remain, either around the correct average value (ŷ = S − T) or with an off-set towards zero (see the curves withŷ = from Eqs. (51) and (52) given the choice of the parameterŷ. The graphs in Figure 4 with the various parameter settings for y confirm these results for small S-values but also show the substantial difference between the approximated integral and the numerically obtained result for E f ,inc . We conclude that the approximated analytic result of Eq. (40) should only be used relatively far from the aperture.
Finally, in Figure 5, we keep the aperture and focus geometry fixed but vary the wave number of the radiation and inspect the behaviour of the complete Rayleigh integral E f (solid blue curve, Eq. (23)), and the Debye integral (dashed red curve, Eq. (44)). In the upper left graph, the number of wave lengths that fits across the aperture is less than 6, in the lower right graph, this value amounts to 588000. In the first case, we have a situation that is frequently encountered in acoustical problems while the very high k-values are typical for optical diffraction problems. For instance, the example with k = 10 5 might correspond to an optical microlens with a diameter of 0.3 mm, used in the visible domain of the spectrum; the last example with k = 10 8 mm −1 may correspond to a telescope objective with a typical diameter of 30 cm. Starting with the low k-value of 10 3 , we immediately see that the (symmetric) Debye integral is a very poor approximation to the field on axis. With the numerical aperture value of 0.2606 for the focusing beam, the agreement is still relatively good, at lower aperture, the divergence between the Rayleigh result and the Debye approximation becomes very important. Earlier studies of the diffraction of a focused beam at low Fresnel number (in our case ≈ 10) using numerical evaluation of the diffraction integral already revealed the asymmetry of the intensity distribution with respect to the geometrical focal point and the shift of the intensity maximum from the focus towards the aperture, see [7], [24], Figure 1 of [25], [30], [31]. At larger k-values, the Debye result gradually approaches the exact Rayleigh result and, in the middle right graph, the difference has become imperceptible. This is due to the fact that we are now many wave lengths away from the diffracting aperture and, in our example, the numerical aperture is also sufficiently large with many Fresnel zones fitting in the aperture. We remark that with the very small wave length in the lower right graph (λ = 63 nm), the central maximum of the focal field is confined to an axial range of typically 3.5 µm. In many applications, the exact location of this maximum should be known and controlled to within 1% of this range. Well-known examples of such a precise focal setting are found in highresolution imaging for optical microlithography and inspection microscopy. For that reason, we have checked the possible divergence between the Rayleigh and the Debye integral up to this very high k-number of 10 8 m −1 . Figure 5 convincingly illustrates why the Debye integral can be used without any problem when imaging with classical optical systems. The lateral extent of the lens aperture is always much larger than the wave length of the incident radiation and the aperture size and the numerical aperture should be sufficiently large to achieve the required image size and image resolution. Together, these conditions lead to a very large number of Fresnel zones in the aperture and allow a safe application of the Debye integral.