General diffraction integral for converging freeform vectorial fields

Herein it is presented a new vector diffraction integral capable to compute the diffraction patterns for aspheric and freeform wavefronts, with different pupil shapes, and different polarization states. The aforementioned vector diffraction integral is based on the energy conservation theorem, the far field approximation, and the angular spectrum representation of optical fields. The integral is validated through illustrative examples whose results are as expected and coincides with the related works.


Introduction
RayleighSommerfeld and scalar Kirchhoff beam propagation models are convenient ways to describe the propagation of light in paraxial conditions.Nevertheless, these wavefront propagation methods cannot be used for high-numerical aperture optical systems.This is due to the scalar nature of these methods which do not account for the state of polarization of light.Therefore, the propagation of light phenomena should be extended to account for the vector character of the light.The first formalism that considers the vector character of the fields is Ignatowski's study of Maxwell's equations [1,2].Then, Kottler proposed a model considering the problem as an integral equation using Green's function techniques [3].Later, the method of Kottler was generalized by Stratton and Chu [4].Richards and Wolf obtained another diffraction integral to describe the focusing of light [5][6][7][8].The foundations of this integral are energy conservation and far-field approximation.It is important to remark that all these formalisms have been conceived for monochromatic light.
This paper introduces a new vector diffraction integral, not yet published according to the best authors' knowledge, capable of computing the diffraction patterns for aspheric, freeform wavefronts in different shapes of pupils and polarization states.This vector diffraction integral is based on energy conservation, the far-field approximation and the angular spectrum representation of optical fields.This manuscript is divided as follows: In section 2 the mathematical derivation of this new vector diffraction integral is presented.In section 3, examples of diffraction patterns computed with the aforementioned integral are shown and in section 4 are presented the conclusions.

Mathematical model
The model presented here is based on the angular spectrum representation of the optical fields.Angular spectrum representation mathematically describes optical fields on homogeneous media using concepts of the Fourier analysis.The idea is that optical fields are described as a superposition of plane waves and evanescent waves, which both are solutions to Maxwells equations.Those plane waves contribute to the field and their direction of propagation coincides with the direction of geometric optical rays.Acknowledging the aforementioned considerations, the electric field E(r, f, z), where (r, f, z) are cylindrical coordinates, at focus can be expressed as , , 0 exp cos sin d d 1 is the Fourier transform of the E(r, f, z) respect to k x and k y at z = 0. k x , k y and k z are the Cartesian transverse spatial frequencies.If we set the far field approximation and the energy conservation theorem [28] over equation (1) to eliminate the evanescent waves, equation (1) turns to, x y z x y 0 E 0 is a constant amplitude setting the far field approximation and the energy conservation theorem.k cos , sin , is the position vector of the focus.
x y x y x y 0 .q(k x , k y ) is the apodization factor, which is the factor obtained from the energy conservation theorem, l 0 (k x , k y ) is the amplitude distribution of the beam and ˆ( ) k k a , x y is a unit vector representing the polarization direction of the electric field near the focus.Equation (2) seems very abstract and cumbersome since the integration is over the frequency domain defined by k x , and k y .It will be more understandable if equation (2) will be in terms of the spatial domain since the wavefronts are represented in the spatial domain.In other words a wavefront is typical represented by a function w(x, y), that depends in two independent variables x and y instead of k x , and k y .The main goal of this manuscript is to express equation (2) in a more intuitive way in the spatial domain.In figure 1 we observe that equation (2) considers only plane waves, then [k x , k y , k z ] are given just by real values.Therefore, acknowledging the Malus-Dupin theorem [29], we can get the propagation vector of a wavefront w(x, y) as the normal vector towards w(x, y), where ∂ x w ≡ ∂w/∂x is the partial derivative of w respect to x and ∂ y w ≡ ∂w/∂y is the partial derivative of w respect to y. k x , k y and k z are the cosine directors in x, y and z directions, respectively, of the wave vector k.In order to express the integration in terms of x and y instead of k x , and k y .The Jacobian can be computed using equation (3), and changing the integrands [30], where ∂ xy w is the derivative of w respect to x and y.Let us focus on â.In fact â can adopt several forms depending on the polarization state of the wavefront; the only condition required is that â must be orthogonal to k.
According to Rodrigues' rotation formula, [31], all the possible vectors orthogonal to k can be represented as   Rodrigues' rotation formula rotates a vector given an axis and angle of rotation.The axis in this case is the propagation vector k and the vector rotated is the polarization state ˆ( ) a x y , , θ is the rotation angle along the plane perpendicular to k. Depending on the value of θ is the polarization state.For example, for a radially symmetric beam, if θ = 0, we have a radial polarization, or if θ = π/2 we get an azimuthal polarization, as shown in figure 2.More generally θ(x, y) can be a function that depends on the position and expresses the polarization aberration function.We recall that â in equation ( 5) is a unit vector.Notice that equations (3) and (5) are orthogonal.
Finally, if we replace equations (3), ( 4) and (5) in equation (2), we get the new vector diffraction integral proposed for this article, Equation (6) is the most important equation of the manuscript.It describes the diffraction pattern of a polarized freeform wavefront.The wavefront must be a continuous function w(x, y), for example, it can be expressed in terms of Zernike polynomials or any other real smooth continuous function.The polarization depends directly on the function θ(x, y).Without losing generalization, we added the mask P(x, y) in equation (6).P(x, y) can be the information on the mask's shape.The mask's shape can be elaborated as P(x, y).By default, equation (6) has a rectangular pupil shape since it is integrated along dx and dy.
In order to align an input normalized polarization vector [ ] e e , , 0 x y T with the polarization vector ˆ( ) x y a , , the polarization angle θ should be,  ( ) Replacing equation (7) in equation ( 6), ´- x y

Illustrative examples
In this section, we present several illustrative examples of diffraction patterns directly computed with equation (6).The caption of each figure indicates the input values of the example being considered.For all the examples q(x, y) = 1 and l 0 (x, y) = 1.
In order to validate equation (6), we show in figure 3 three examples of different radially symmetric wavefronts.In red and in dashed black is the normalized intensity of the diffraction pattern of spherical wavefronts according to [8] and equation (6), respectively.In orange is the normalized intensity of the diffraction pattern of an aspherical wavefront.Figures 3(a) and (b) are the profiles with respect to z and r, respectively.Notice that he red and dashed black lines overlap, which means that the prediction made by equation (6) consists of [8].The pattern denoted by the orange line is notably larger than the red and dashed black since it is not spherical.From the figure, it is clear that the red line and the dashed black line overlap.The reason for the overlapping is that with the given conditions of this particular example, equation (6) can be reduced to Richards and Wolf integral [8].All the wavefronts in this example are radially polarized.Other classical results can be also obtained with equation (6), for example, if the field is linearly polarized using equation (9).
In figure 4, we present the normalized intensity of the diffraction pattern of a freeform wavefront with a polarization given by ( ) q = + x y x y , 2 2 .The wavefront with its respective polarization states is presented in top.The diffraction pattern in the f f r r cos sin axis for the components along x, y, and z are presented in bottom left, bottom centre and bottom right, respectively.Due to the aberration suffered by the wavefront, the diffraction patterns along x and y differ as expected.Also at the centre of both diffraction patterns, the electric field is zero.Next is the diffraction pattern of a linearly polarized wavefront using equation (9).In other words , , 0 1, 0, 0 x y T T .The linearly polarized wavefront can be found in figure 5 (top).In figure 5 bottom left, bottom centre, and bottom right is the diffraction patterns in the image plane along x, y, and z, respectively.To the best of the authors' knowledge, this classic example, the linearly polarized wavefront, has not been presented with aberrations in the literature.If the aberrations are removed, the diffraction pattern converges to the reported results in the literature.In this example, the diffraction pattern is computed with equation (9) rather than equation ( 6), since it is straightforward to express a linearly polarized wavefront with equation (9).
Figure 6 it is presents an example inspired in [32].Figure 6 shows the diffraction pattern generated by an offaxis parabolic mirror with a rectangular aperture, thus the wavefront is spherical.For this example, we use radial polarization.The form of the pattern is not spherical since the aperture is rectangular, so the intensity profile respect to the z in the image plane.The same occurs with the intensity profiles with respect to the x and y.The rectangular shape of the intensity profiles resembles the results presented in [32].
Equation ( 6) is very robust since it can compute the diffraction pattern of freeform wavefronts with different polarization states and pupils.Thus it can be considered the general diffraction integral for freeform wavefronts, different polarization states, and arbitrary pupils.But it is important to see that equation ( 6) is valid if and only if the Jacobean of a w(x, y) is not zero.For example, a single plane waves its Jacobean is zero, thus equation (6) can not be used.This result is expected since a single plane wave does not form a focus.Equation (6) can be considered as a great bridge between the two paradigms, geometrical optics and diffraction theory since it is obtained from the angular spectrum representation of the optical fields and the vector k has the direction of propagation of geometrical rays.We observe that equation (6) guarantees that k and â are always orthogonal despite the possible values of θ(x, y).This aspect is crucial and usually ignored by other diffraction integrals [18,19,[33][34][35][36].The method presented here differs from the paradigm presented in [18,35,36].The method in [18,35,36] considers several approximations like small aberrations or small angles in order to simplify the mathematical expressions.[18,35,36] are limited to Seidel aberration under the aforementioned approximations.Our method works directly with an explicit wavefront function w(x, y), its normal vectors (propagation), and its tangential vectors (polarization).The derivation of the equations in [18,35,36] is cumbersome, and needs to use additional parameters like optical coordinates.These additional parameters are need it because it uses cylindrical and spherical coordinates.Our method is straightforward since it's only based on Cartesian coordinates.
Due to the reason of space, we do not include an example with a more complicated aperture shape.But from equation ( 9) is clear that the shape of the aperture can change if P(x, y) is defined to follow a particular aperture shape.Equation ( 9) is straightforward, it is an integral along the aperture shape denoted by P(x, y) and y min , x min , y max and x max .Also, it is important to remark that the elliptical or circular polarization states can be calculated with the sum of different diffraction integrals.In the same way that Richards and Wolff diffraction integral considers elliptical or circular polarization states.
The method presented here differs from the numerical approach presented in [37].The method in [37] is based on numerical approximations like the homeomorphic Fourier transform.The deduction of closed-form equations (6) and ( 9) is entirely analytical and it was done with no numerical approximations nor iterations.
Despite all the presented examples being quadratic wavefronts.There is no limitation for equations (6) and (9) of the shape of the wavefront, as long as it is a real smooth continuous function.But the more complicated the wavefront, the longer time is the computation.For the example presented in figure 5 the time needed for the computation using 36 hours in an MSI GF63 THIN 11UC laptop using Wolfram Mathematica 13.
Observe that we set q(x, y) = 1 and l 0 (x, y) = 1 in the examples because we are interested only behaviour of the integral rather than the illumination pattern or the apodization function.The illumination pattern can influence a lot the diffraction pattern.For example, the diffraction pattern of two equal wavefronts with different illumination patterns may change significantly.The same applies to the apodization function.In our last work, we presented an integral that computes the diffraction pattern of the freeform wavefront with rotationally polarized [38].The mentioned integral can be seen as a special case of equation (6) when θ = 0. Therefore all the properties of the mentioned integral are presented in equation (6) when θ = 0. Equation ( 6) is more robust than its predecessor since it not only considers radial polarization and azimuthal polarization but also consider any polarization state since in equation (6), θ may change point by point since it is a function, θ(x, y).In the same work, we mentioned that the extension of the method to other polarization states is straightforward which in fact is done with equations (6) and (9).The creativity and abstract reasoning to get equations (6) and (9) were not trivial, but now we can confirm that we have extended the method [38] to other states.
The sampling of w(x, y) is not related to numerical aperture optical systems since equations ( 6) and (9) directly work with the polarization state ˆ( ) x y a , and the propagation vector k.The numerical aperture can easy compute with the propagation vector k of the wavefront in the boundary of the aperture.In other words, the numerical aperture can be computed with the marginal ray.Equations (6) and (9) do not change their form where when is located the marginal ray in a given aperture.But it is important to remark the accuracy of equations ( 6) and ( 9) is proportional to the sampling of w(x, y).The higher the number of local maxima/minima in w(x, y), the higher the sampling of w(x, y).
Notice that equations (6) and (9) are the same integrals but expressed in different forms equation ( 6) is friendly to express polarization state depending on the angle θ(x, y).For example radial polarization or azimuthal polarization.Equation ( 9) is more accessible for the input polarization vector [ ] e e , , 0 x y T when the polarization is easier to desire in terms of the Cartesian components of the input polarization vector [ ] e e , , 0 x y T .Equations ( 6) and (9) can compute the diffraction pattern of commercial optical systems.The mask P(x, y) can fit the dimensions of a commercial optical system.w(x, y) and θ(x, y) can be obtained from commercial optical design software.This is the natural step to follow in this research and it is out of the scope of this article.

Conclusion
In this manuscript, it was introduced new diffraction integral that considers the shape of the wavefront w(x, y), the polarization state θ, and the form of the mask P(x, y).This integral is mathematically described by equations (6) or (9).Equations (6) and (9) have been tested for several scenarios and the results are as expected, and particularly coincide with [8].For example, the wavefront is spherical and radially polarized.Also when the wavefront is spherical and the field is linearly polarized.Equations (6) and (9) represent the same integral but the first equation is more friendly to be computed for the polarization that depends on the angle θ(x, y) like radially and azimuthally polarized wavefronts.Equation ( 9) is more friendly to be computed for input polarization vector [ ( ) ( ) ] e x y e x y , , , , 0 x y T .It is important to mention that both integrals are valid if and only if the wavefront w(x, y) is convergent and it is a real smooth continuous function.
For its robustness and versatility equation (6) can be considered as the general diffraction integral for freeform wavefronts, different polarization states, and arbitrary pupils.The next step to follow in this research is to compute the diffraction pattern of commercial optical systems with equation (6).

Figure 1 .
Figure 1.Diagram of the diffraction phenomena.The wavefront w(x, y) is placed in the coordinate system (x,y).The vector k is normal to w(x, y) and â is normal to k.The diffraction pattern is in the coordinate system ( ) f f r r cos , sin and z is the axis of propagation.

Equation ( 9 )
is similar to equation(6), but it is more accessible for the input polarization vector [ ] example the linear polarizations or polarization defined point by point, (6) is more practical for radial and azimuthal polarizations.

Figure 3 .
Figure 3. Normalized intensity of the diffraction produced by a spherical wavefront, in red and black is computed with equation (6) and[8], respectively.In orange is the Normalized intensity of the diffraction produced by w(x, y) = 0.35( − 1 + 2x 2 + 2y 2 ).Left is intensity profile with respect to the z axis.Right is the intensity profile with respect to the r axis.The mask function for all the examples in this figure is P(x, y) = {x 2 + y 2 < 1 ⟹ 1, otherwise ⟹ 0} and θ(x, y) = 0.

Figure 4 .
Figure 4. Normalized intensity of the diffraction produced by w(x, y) = cx 2 + (c/2)y 2 , with c = 0.05.Top is the shape of the wavefront.Bottom left, bottom centre, and bottom right are the intensity profile with respect to the x, y, and z axes, respectively in the axis of the image plane.The mask for all the examples in this figure is P(x, y) = {x 2 + y 2 < 8 ⟹ 1, otherwise ⟹ 0} and θ(x, y) = 0 and θ(x, y) = x 2 + y 2 .The domain of the wavefront is x ä {0, 5} and y ä {0, 5}.

Figure 5 .
Figure 5. Normalized intensity of the diffraction pattern produced by w(x, y) = cx 2 + (c/2)y 2 , with c = 0.05.Top is the shape of the wavefront.Bottom left, bottom centre, and bottom right are the intensity profile with respect to the x, y, and z axes, respectively in the axis of the image plane.The mask for all the examples in this figure is P(x, y) = {x 2 + y 2 < 8 ⟹ 1, otherwise ⟹ 0}, and the wavefront is linearly polarized.

Figure 6 .
Figure 6.Normalized intensity of the diffraction pattern produced by spherical wavefront radial polarization θ = 0. Top is the shape of the wavefront.Bottom left, bottom centre, and bottom right are the intensity profile with respect to the x, y, and z axes, respectively in the axis of the image plane.The mask for all the examples in this figure is P(x, y) = 1, and θ(x, y) = 0 and θ(x, y) = x 2 + y 2 .The domain of the wavefront is x ä {−0.7, 0.7} and y ä{−0.7,0.7}.