Centroid propagation through optical systems with ABCD kernels and nonuniform or finite apertures

Jesús Lancis, Raúl Martínez-Cuenca, Jorge Ares, and Salvador Bará* Grup de Recerca d’Òptica de Castelló (GROC), Departament de Física, and Institut de Noves Tecnologies de la Imatge, Universitat Jaume I, 12071 Castelló, Spain Departamento de Física Aplicada, Universidad de Zaragoza, 50009 Zaragoza, Spain Área de Óptica, Departamento de Física Aplicada, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia, Spain *Corresponding author: salva.bara@usc.es

As an example of practical relevance, many wavefront sensing devices include a relay optical system in order to transport and possibly rescale onto the CCD camera the raw irradiance distributions present at the back focal plane of the microlens array lying at the core of the Hartmann-Shack sensors.In the design and implementation of these transport systems it becomes essential to ensure that the centroid position information is preserved along the light path.
To study the propagation of an irradiance centroid throughout any optical system, it is always possible to perform a direct calculation of the detailed structure of the diffracted wave field, and from it to compute first the irradiance distribution and subsequently the centroid.However, for a wide class of systems there is a more direct route: In this paper we show that if the propagation of a light field can be described by a diffraction integral with an ABCD kernel, then the rigorous propagation of its irradiance centroid will take place in the geometric regime and will be completely described by the corresponding ABCD ray-transfer matrix, exactly in the same way as if the centroid path were a conventional geometrical ray.For propa-gation through homogeneous or weakly inhomogeneous media, this result is a particular instance of the optical Ehrenfest's theorem [20][21][22][29][30][31].We also show, however, that potentially significant deviations from this geometrical propagation rule may arise in the presence of finite or nonuniform apertures truncating or otherwise modifying the input beam irradiance distribution.

CENTROID PROPAGATION IN ABCD SYSTEMS
The propagation of geometric rays between two planes through first-order optical systems can be described in terms of the well-known ABCD ray-transfer matrix formalism [32].A ray is characterized at each plane by its position (coordinate of the impact point) and momentum (essentially the angle that its direction of propagation makes with the optical axis).The ABCD matrix relates the position and momentum at the output plane with those at the input plane by means of a linear transformation.Since its determinant is equal to the ratio of the refractive indices at both ends of the system, this matrix possesses only three independent parameters.In particular, if these refractive indices are equal, the matrix turns out to be unitary (AD − BC ¼ 1).The ABCD matrix can be decomposed in different ways as the product of three simpler matrices, each one accounting for an elementary transformation: propagation through free space, magnification, and refraction by a lens [33].It has also been shown that the wave propagation between the input and output planes in an ABCD system can be accounted for by an integral transform with an ABCDdependent quadratic kernel, as a generalization of the usual Fresnel transform associated with diffraction in a homogeneous medium [32,33].There is a wide class of optical systems whose behavior can satisfactorily be modeled by an ABCD matrix within the paraxial domain.
Let us consider a generic ABCD optical system in which the diffracted field u o ðr o Þ at the ouput plane can be computed in terms of the input field uðrÞ by the superposition integral where the integration is extended to the whole input plane; d 2 r ¼ dxdy is the differential area element at that plane; and Kðr o ; rÞ is the quadratic diffraction kernel given by [33]: where λ is the wavelength and δðrÞ is a two-dimensional Dirac δ distribution.The first expression, for B ≠ 0, corresponds to a general system, whereas the second one, ðB ¼ 0Þ, applies to input and output planes optically conjugated; i.e., it corresponds to the imaging condition.
The transverse position x o and optical momentum p o of the irradiance centroid at the output plane are given by the expected value of the corresponding position (x o ≡ r o ) and momentum (p o ≡ ðikÞ −1 ∇ o ) operators as (see Appendix A) In order to deduce the expressions describing the transfer of the centroid from the input to the output plane in the general case ðB ≠ 0Þ, one can substitute Eqs. ( 1) and (2) into the third members of Eqs. ( 3) and (4), exchange the order of integration by performing first the integral over r o , take into account several basic properties of the Dirac δ distribution and its derivatives (see Appendix A), and, after a little algebra, one finally obtains: so that the output position and momentum of the irradiance centroid, given some initial x and p values, are just those that would correspond to a classical geometric ray with the same initial conditions.This simple and useful result provides a fast and accurate way to study the centroid propagation throughout any true ABCD optical system without the need to resort to the explicit calculation of the diffracted wave fields.Writing Eq. ( 5) as a vector-matrix relation where the elements of the vectors are vectors themselves is a convenient shorthand notation for two similar formulas, one applying to the Xand other to the Y -components of the positions and momentums.This equation is then to be interpreted with regard to its components.The same simplifying notation will be used henceforth.Equation ( 5) is also straightforwardly obtained for the imaging case by direct integration of Eq. ( 1) using the B ¼ 0 kernel, computing directly from it both the irradiance I o ðr o Þ and the optical path length W o ðr o Þ, and substituting them into the last members of Eqs. ( 3) and ( 4).See also Appendix A for details.

CENTROID PROPAGATION THROUGH THIN OPTICAL ELEMENTS WITH NONUNIFORM TRANSMITTANCES AND FINITE APERTURES
The presence of optical elements with nonuniform transmittances and/or finite apertures truncating or otherwise modifying the irradiance profile of the input beam in general will break this simple geometrical propagation rule.To see how they do modify the centroid propagation, let us consider an optically thin element characterized by the complex transmittance function tðrÞ ¼ jtðrÞj exp½ikW t ðrÞ such that the incident field passing through it is transformed as u o ðr o Þ ¼ tðrÞuðrÞ, with r ¼ r o .Consequently, the exit irradiance distribution becomes I o ðr o Þ ¼ TðrÞIðrÞ, where TðrÞ ¼ jtðrÞj 2 , and the local wavefront slopes are modified as Substituting these results into Eqs.( 3) and ( 4) and taking into account that now the radiant flux is in general not conserved, i.e., E o ¼ R TðrÞIðrÞd 2 r ≠ E, it is easy to obtain the following transformation formulas: The upper element of the last summand of Eq. ( 6) reveals the existence of an irradiance-and transmittance-dependent centroid shift, thus departing from the expected geometrical behavior x o ¼ x.Note that this shift does not depend on the phase of the incident beam nor on the phase added to it by the optical element.The integral expressions in the lower element of this last summand convey information about the bias in the output momentum: the first integral can be interpreted as an effective correction to the input momentum, due to the nonuniform weighting of the input irradiance by the factor TðrÞ, whereas the second one, in turn, accounts for the additional bending of the centroid ray path due to the refractive action ∇W t ðrÞ of the optical element, weighted by TðrÞIðrÞ.
In the case of a thin spherical lens of focal length f ðW t ðrÞ ¼ −r 2 =2f Þ, this last term corresponds to the bending that the lens would impart to a classical geometrical ray incident on it not at the input centroid location, x, but at the output shifted one, x o .These last statements can be made more clear by defining and substituting ∇W t ðrÞ ¼ −r=f into Eq.( 6) to obtain which supports the interpretation stated above.Equation ( 8) can alternatively be rewritten as the sum of the contributions of an ideal ABCD lens plus an additional output shift: where Δx o and Δp o , the output shifts at the exit plane of the lens, are explicitly given by After leaving the element, the shifted and deflected centroid will propagate again in free space according to the usual rule Given the linearity of Eq. ( 11), at any distance z from the element, the actual centroid will be offset from its classical geometric position and momentum (those corresponding to so that the bias in the momentum is constant, whereas the bias in the centroid position is a linear function of z. A particularly interesting case of the axial evolution of the centroid shift is that of a thin spherical lens illuminated by a quadratic wavefront.For such a lens Eq. ( 10) can be substituted into Eq.( 12).If, additionally, the incident wave has a quadratic phase such that W ðrÞ ¼ ðr − r 0 Þ 2 =2s, where s is the paraxial curvature radius of the wavefront at the input plane of the lens and r 0 is the transversal position vector of its curvature center, then by Eq. ( 7) we have Δp ¼ ð1=sÞΔx.Substituting this into Eq.( 10) and regrouping terms we finally get Although formally written Δp o ðzÞ, in this case the momentum bias does not depend on z, a not unexpected result, since after leaving the lens the centroid propagates in a homogeneous medium, and hence its propagation direction is preserved.The centroid shift Δx o ðzÞ, in turn, does depend linearly on z.Note that when the observation plane is optically conjugated to the plane containing the center of curvature of the incident wavefront, that is, when the object-image relationship ð1=z þ 1=s ¼ 1=f Þ holds, the centroid shift identically cancels out: no matter how big the centroid bias Δx may be at the output of the lens, the position of the actual and the ideal centroids will coincide in the image plane.On the other hand, Eq. ( 13) shows that both the centroid position and momentum shifts in the image space are directly proportional to Δx: if there is no position shift at the lens plane, the centroid of the quadratic beam will be transformed by the lens as a classical geometric ray.
A simple numerical example may be helpful to give some insight into the magnitude and potential significance of the Δx o and Δp o output shifts.Let us consider a thin unaberrated converging lens limited by a circular pupil at f =5 (focal length 10 units normalized to its aperture radius R) and a set of rotationally symmetric TEM 00 Gaussian wavefronts of variable waist size w 0 propagating in the direction of the optical axis (Z), whose centers can be located at any transversal position with respect to the lens center (see Fig. 1).Because of the rotational symmetry of the irradiance distributions of these Gaussian beams, the input irradiance centroid is always coincident with geometrical center of the beam.Besides, due to the rotational symmetry of their phases around the Z axis, the input momentum is in all cases identically zero.For the sake of simplicity, let us further restrict the analysis to the particular case in which the beam waist is precisely located onto the lens plane, so that the input wavefronts are flat at this particular axial location.Figures 2 and 3 show how the centroid position and momentum behave at the exit plane of the lens, in comparison with what would be expected from the classical geometric behavior of an unobstructed ABCD system.Since both the incoming beam and the lens are rotationally symmetric in amplitude and phase, these figures drawn as onedimensional (1D) plots are valid for any radial direction in a reference frame with its origin located at the lens center.All results have been obtained by direct numerical evaluation of the corresponding input (x; p) and output (x o , p o ) centroid parameters.
In Fig. 2 we plot the output versus input centroid positions for a set of five Gaussian beams with σ (half-width at the 1=e irradiance level) ranging from 0.2 to 1 pupil radius in equally spaced steps.The expected geometric behavior (x o ¼ x) is represented by the unlabeled diagonal straight line.Figure 3 shows the modulus of the output momentum p o as a function of the input centroid position for this same set of wavefronts, after being refracted by the f =5 converging lens.As in the figure above, the unlabeled straight line corresponds to the ideal ABCD behavior.All figures were drawn with the input centroid x ranging from 0 (i.e., located at the lens center) to 1 (i.e., located on the pupil rim).The magnitude of the potential biases is clearly apparent.With the beam and lens parameters used in this particular example and σ ¼ R, the output centroid at the exit plane of the lens can get shifted a distance greater than 0:6R toward the lens center, and its momentum (angle of propagation with the Z axis), can be more than 60 mrad smaller (in absolute value) than what would be expected had the lens aperture not truncated the input beam.Note that both the position and momentum biases are bigger for beams with higher values of σ, since the effects of truncation affect them in a more drastic way.Smaller σ beams, in turn, are less affected by truncation, especially when they are away from the pupil rim, and the finite-aperture thin lens behaves in these cases more closely to a true ABCD system for a wider range of input centroid positions.For the sake of simplicity, in this numerical example we have only analyzed the dependence of the output position and momentum shifts on the position of the input centroid.Similar calculations may be carried out to study the bias produced by variable input momentums.
Once the beam leaves the lens, its centroid-laterally shifted and angularly deflected with respect to the ideal one-propagates through free space according to the conventional evolution equation given by Eq. (11).As expected from Eq. ( 13), the actual centroid path crosses the ideal one at the image point of the center of curvature of the wavefront incident on the entrance plane of the lens.Since in this example the input wavefront is flat and propagating along Z, the   crossing occurs precisely at the back focal point of the lens (see Fig. 4).
This particular example helps to visualize in an intuitive way the main distortions that the presence of apertures and uneven irradiances introduces in the centroid propagation paths after refraction by a lens, quantitatively described by Eq. ( 8).On the one hand, the centroid at the exit plane of the lens gets laterally shifted due to the hard spatial limits imposed by the lens aperture and the corresponding change of the centroid integration area [resulting in the spatial shift Δx given in Eq. ( 7)].On the other hand, its momentum or direction of propagation is biased with respect to the one expected from geometrical optics, due to the combination of two factors: a bias independent of the lens power (−1=f ), arising due to the spatial limits imposed by the lens aperture on the local distribution of the input momentum [resulting in an effective correction, Δp, shown in Eq. ( 7)], and, additionally, a bias due to the refractive action of the lens, which bends the centroid path as it would bend a geometrical ray incident on it, not at the input position x, but at the shifted one, x þ Δx.The overall consequence is that the centroid propagation path in the image space will also be biased, excepting-if both the lens and the incoming wavefront are free from aberrations-at the image plane of the center of curvature of the incoming wave.

ADDITIONAL REMARKS AND CONCLUSIONS
The propagation of irradiance centroids through optical systems with general ABCD kernels can be satisfactorily described using the simple formalism of ray-transfer matrices, without resorting to the explicit computation of the detailed structure of the diffracted wave fields.Centroid paths in ABCD systems behave the same way geometrical rays would.This useful property drastically simplifies the first steps of the design of optical systems for centroid detection, steering and imaging.It also provides some rationale to support the description of wavefront slope measuring devices (such as Hartmann-Shack wavefront sensors and laser raytracing aberrometers) in terms of geometrical optics.This description, often met in the literature in terms of classical geometric rays, which is in itself only a rough approximation, can however be rigorously applied with the proviso that (a) the optical subsystems composing the wavefront sensor effectively behave as true ABCD systems, and (b) instead of referring to classical geometric rays, the description shall refer to the centroid propagation path, taking into account that its slope (optical momentum) is given by the irradiance-weighted spatial average of the local wavefront slopes.
The presence of optical elements with finite apertures truncating or otherwise modifying the irradiance distribution of the input beam gives rise to deviations from this simple geometrical picture.After crossing an element of uneven transmittance, the centroid position is generally shifted and its slope deflected with respect to the values expected by geometrical optics.This gives rise to a centroid path discontinuity at the element plane.Further propagation through free space translates into a centroid position bias linearly dependent on z.Interestingly enough, for incident beams with quadratic phases (spherical, Gaussian, etc.) refracted by thin unaberrated lenses, this bias cancels out at the image plane of the curvature center of the input wavefront.
These facts provide some interesting clues for the proper design of Hartmann-Shack wavefront sensors.These devices very often include a relay optical system to transport, rescale, and image onto the detector (usually a CCD camera) the irradiance spots formed at the back focal plane of the microlens array.In the case of no significant truncation of the focal beams by the finite aperture of the relay optical elements, no special caution is required to deal with the centroids of the beams leaving it.In the case of truncation, however, the possible existence of centroid shifts Δx o at the CCD plane has to be evaluated.In accordance with the results in Section 3, as far as the beams propagating through the relay keep a quadratic phase and the lenses composing it are aberration-free, there will be no centroid position bias at the CCD detector if the CCD plane is optically conjugated to the plane containing the center of curvature of the wavefronts incident on the relay.For a usual Hartmann-Shack setup, this means that truncation by the relay lenses will not be much of an issue as far as (a) the beamlets produced by each microlens of the array are sensibly spherical ones, (b) they all focus on a common plane, and (c) the CCD detector is optically conjugated to that plane.It is well known that this is not the only possible configuration for wavefront slope sensing: if there is no truncation, the wavefront slopes can be determined from the centroids measured at any plane located after the microlens array, knowing the precise distance to this plane.In case of potential truncation, however, imaging the focal spots onto the CCD is advisable in order to avoid the kind of bias described in this work.The potential existence of a momentum bias is also relevant for those slope detection schemes based on the measurement of the beam centroid at two axially displaced planes, in what is sometimes called a "differential wavefront sensor" configuration [34].In these cases, it is advisable to ensure in advance that all optical elements composing the relay system have enough aperture to allow the propagation of all foreseeable incident beams without being affected by significant truncation.It should be noted, though, that usual wavefront sensing configurations (i.e., those using a known reference wave to determine the initial position of the microlens foci and measure displacements relative to them) may be more resilient to the deleterious effects of finite apertures or nonuniform transmittances, due to the fact that some amount of partial cancellation of the biasing effects can reasonably be expected.
Whether or not a given optical system will behave as a true ABCD one or will give rise to beam truncation does not only depend on the kind and configuration of the optical elements composing it, but also on the wavefronts under study.Strictly speaking, there are no absolute ABCD systems: it is always possible to find a subset of beams whose centroids will be significantly affected by the finite apertures of their optical elements, departing from the expected geometrical propagation.However, if the beam interaction with the limiting apertures is sufficiently weak, a wide range of optical systems can be deemed to behave as true ABCD ones for many practical applications.

APPENDIX A
The demonstration of the last equality in Eq. ( 4 Equation ( 5) for the nonimaging case ðB ≠ 0Þ easily follows after some algebra by substituting Eqs. ( 1) and (2) into the first integrals of Eqs. ( 3) and ( 4), exchanging the order of integration of the resulting six-dimensional integral to perform in the first place the integration on r o , and taking into account that from the well-known equality for the 1D Dirac δ distribution, its nth order derivatives, δ ðnÞ ðxÞ, can be expressed as [22] δ ðnÞ ðxÞ ≡ d n δðxÞ dx n ¼ and hence the definite integrals of powers of α multiplied by linear exponentials with infinite integration limits appearing in the calculation of the centroid position x o and its optical momentum p o (with n ¼ 0; 1) can be formally evaluated as Note that there is a misprint affecting the c factor in the corresponding equation listed as (A6) in the appendix of Ref. [22].Once again, the fact that the resulting magnitudes are real-defined implies that any remaining imaginary term cancels out.
For the imaging condition ðB ¼ 0Þ, the derivation of Eq. ( 5) is simpler.After substituting the corresponding kernel into Eq.( 1) and integrating, we have so that the corresponding amplitudes and phases are related as and Eq. ( 5) straightforwardly results by direct substitution of Eq. (A7) into the last terms of Eqs. ( 3) and (4).

Fig. 1 .
Fig. 1.Output irradiance distribution for an input Gaussian beam of half-width σ ¼ R (measured at the 1=e irradiance level), centered at x ¼ 0, y ¼ R, and truncated by a lens with circular pupil of radius R.

Fig. 2 .
Fig. 2. Output versus input centroid position for Gaussian beams with the half-width σ values indicated in curve labels.The unlabeled diagonal straight line corresponds to the ideal ABCD behavior.All units normalized to the lens pupil radius, R.

Fig. 3 .
Fig. 3. Output momentum (mrad) versus input centroid position for the wavefronts in Fig. 2 after being refracted by an f =5 unaberrated converging lens.The unlabeled straight line corresponds to the ideal ABCD behavior.Length units normalized to the lens radius R.

Fig. 4 .
Fig. 4. Longitudinal section of the image space showing the centroid propagation paths after beam refraction by a f =5 converging lens assuming the incident wavefront has its input centroid located at the lens pupil rim.Horizontal line, optical axis; slanted dotted line, ideal ABCD behavior; solid lines, actual centroid paths.Labels correspond to the beam half-width values σ.All centroid paths intersect at the image of the center of curvature of the incoming wavefont (in this example the back focal point of the lens at z ¼ 10).All distances normalized to R.
n exp½−iαxdα ¼ 2πi −n δ ðnÞ ðxÞ: ðA3Þ Dirac δ derivatives operate under the integral symbol as[35] Z∞ −∞ f ðxÞδ ðnÞ ðxÞdx ¼ ð−1Þ n f ðnÞ ð0Þ;ðA4Þwhere f ðnÞ ðx 0 Þ ¼ d n f =dx n evaluated at x ¼ x 0 .By a change of variables we also getZ ∞ −∞f ðxÞδ ðnÞ ½cðx − x 0 Þdx ¼ ð−1Þ n jcj −1 f ðnÞ ðx 0 Þ: ) is immediate by substituting u o ðr o Þ ¼ a o ðr o Þ exp½ikW o ðr o Þ into the first integral and applying the gradient operator ∇ o ¼ ð∂=∂x o ; ∂=∂y o Þ.Take into account that p o ¼ hp o i is real by definition, so that the formal imaginary part ðikÞ −1 R a o ðr o Þ∇ o a o ðr o Þd 2 r o remaining at the final expression obtained after performing the calculations is identically zero.