High numerical aperture vectorial imaging in coherent optical microscopes

Imaging systems are typically partitioned into three compo nents: focusing of incident light, scattering of incident l ight by an object and imaging of scattered light. We present a model of high Numeri cal Aperture (NA) imaging systems which di ffers from prior models as it treats each of the three components of the imaging system rigorously. It is well known that when high NA lenses are used the imaging system must be tr eated with vectorial analysis. This in turn requires that the scat tering of light by the object be calculated rigorously according to Maxwell ’s equations. Maxwell’s equations are solvable analytically for only a sm all class of scattering objects necessitating the use of rigorous numer ical methods for the general case. Finally, rigorous vectorial di ffraction theory and focusing theory are combined to calculate the image of the scattered l ight. We demonstrate the usefulness of the model through examples. © 2008 Optical Society of America OCIS codes:(180.0180) Microscopy; (050.1960) Di ffraction theory; (260.2110) Electromagnetic theory; (10.1758) Computational imaging References and links 1. B. Richards and E. Wolf, “Electromagnetic di ffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. R. Soc. London, Ser. A 253, 358–379 (1959). 2. J. Goodman, Introduction to Fourier Optics(McGraw-Hill, New York, 1988). 3. V. S. Ignatowsky, “Di ffraction by a lens of arbitrary aperture,” Trans. Opt. Inst. P etr.1, 1–36 (1919). 4. R. Luneburg, Mathematical Theory of Optics (University of California Press, Berkeley and Los Angeles, 1966). 5. C. J. R. Sheppard and T. Wilson, “The image of a single point i n microscopes of large numerical aperture,” Proc. R. Soc. London, Ser. A379, 145–58 (1982). 6. P. T̈orök, P. Higdon, R. Jǔ skaitis, and T. Wilson, “Optimising the image contrast of conv entional and confocal optical microscopes imaging finite sized spherical gold scatt erers,” Opt. Commun. 155, 335–341 (1998). 7. A.S. van de Nes, P. T̈ orök: “Rigorous analysis of spheres in Gauss-Laguerre beams,” Opt. Express15, 13360– 13374 (2007) 8. J. Judkins and R. Ziolkowski, “Finite-di fference time-domain modeling of nonperfectly conducting metall ic thinfilm gratings,” J. Opt. Soc. Am. A. 12(9), 1974–1983 (1995). 9. J. Liu, B. Xu, and T. Chong, “Three-Dimensional Finite-Di fference Time-Domain Analysis of Optical Disk Storage System,” Jpn. J. App. Phys. 39, 687–692 (2000). 10. L. Liu, Z. Shi, and S. He, “Analysis of the polarization-d ependent di ffraction from a metallic grating by use of a three-dimensional combined vectorial method,” J. Opt. Soc. Am. 21, 1545–1552 (2004). 11. J. Stamnes, Waves in Focal Regions (Adam Hilger, Bristol and Boston, 1986). 12. P. T̈orök, “An imaging theory for advanced, high numerical aperture o ptical microscopes,” (2004). DSc thesis. 13. P. T̈orök, P. Varga, Z. Laczik, and G. R. Booker, “Electromagnetic di ffraction of light focused through a planar interface between materials of mismatched refractive indices : an integral representation,” J. Opt. Soc. Am. A 12, 325–332 (1995). #88564 $15.00 USD Received 15 Oct 2007; revised 30 Dec 2007; accepted 4 Jan 2008; published 7 Jan 2008 (C) 2008 OSA 21 January 2008 / Vol. 16, No. 2 / OPTICS EXPRESS 507 14. P. T̈orök and P. Varga, “Electromagnetic di ffraction of light focused through a stratified medium,” Appl. Op t. 36(11), 2305–2312 (1997). 15. J. Jin,The finite element method of electromagnetics (Wiley Interscience, 2002). 16. H. Pocklington, “Electrical oscillations in wires,” Pr oc. Cam. Phil. Soc. 9, 324–332 (1897). 17. S. Yee, “Numerical solution of initial boundary value pro blems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). 18. O. Martin, A. Dereux, and C. Girard, “Iterative scheme for c mputing exactly the total field propagating in dielectric structures of arbitrary shape,” J. Opt. Soc. Am. A 11, 1073–1080 (1994). 19. R. Luebbers, F. Hunsberger, K. Kunz, R. Standler, and M. S chneider, “A frequency-dependent finite-di fference time-domain formulation for dispersive materials,” IEEE Trans . Electromag. Compat. 32, 222–227 (1990). 20. P. Johnson and R. Christy, “Optical constants for noble me tals,” Phys. Rev. Lett. 6, 4370–4379 (1972). 21. K. Yee, “Numerical solution of inital boundary value prob lems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). 22. A. Taflove and S. Hagness, Computational electrodynamics, second edition (Artech House, 2000). 23. A. Bayliss and E. Turkel, “Radiation boundary condition s for wave-like equations,” Commun. Pure App. Math. 23, 707–725 (1980). 24. B. Engquist, “Absorbing boundary conditions for the nume rical simulation of waves,” Math. Comput. 31, 629– 651 (1977). 25. G. Mur, “Absorbing boundary conditions for finite-di fference approximation of the time-domain electromagneticfield equations,” IEEE Trans. Electromag. Compat. 23, 377–382 (1981). 26. J.-P. Berenger, “A perfectly matched layer for the absorp ti n of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). 27. A. Poggio and E. Miller, Computer techniques for electromagnetics , hap. Integral equation solutions of threedimensional scattering problems, pp. 159–264 (Pergamon Press , 1973). 28. P. T̈orök and C. Sheppard, High numerical aperture focusing and imaging (Adam Hilger). 29. P. T̈orök, P. R. T. Munro, and E. Kriezis, “A rigorous near to farfield transformation for vectorial di ffraction calculations and its numerical implementation,” J. Opt. Soc. A m 23, 713–722 (2006). 30. G. Arfken,Mathematical Methods for Physicists , 3rd ed. (Academic Press, Boston, 1985). 31. P. T̈orök, “Propagation of electromagnetic dipole waves through di electric interfaces,” Opt. Lett. 25, 1463–1465 (2000). 32. B. Karczewski and E. Wolf, “Comparison of three theories o f electromagnetic di ffraction at an aperture Part I: coherence matrices, Part II The far field,” J. Opt. Soc. Am. 56, 1207–19 (1966). 33. P. R. T. Munro and P. T̈ orök, “Calculation of the image of an arbitrary vectorial elect romagnetic field,” Opt. Express15, 9293–9307 (2007). 34. L. G. Schulz and F. R. Tangherlini, “Optical constants of silver, gold, copper, and aluminum. II. The index of refraction n.” J. Opt. Soc. Am. 44, 362–368 (1954). 35. P. T̈orök, P. Higdon, and T. Wilson, “Theory for confocal and conven tio al microscopes imaging small dielectric scatterers,” J. Mod. Opt. 45, 1681–1698 (1998). 36. P. R. T. Munro and P. T̈ orök, “Vectorial, high-numerical-aperture study of phase-co ntrast microscopes,” J. Opt. Soc. Am. A.21, 1714–1723 (2004). 37. P. R. T. Munro and P. T̈ orök, “Vectorial, high numerical aperture study of Nomarski’s d ifferential interference contrast microscope,” Opt. Express 13, 6833–6847 (2005). 38. M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave a nalysis of planar-grating di ffraction,” J. Opt. Soc. Am. 71, 811–818 (1981). 39. E. Moreno, D. Erni, C.Hafner and R. Vahldieck, “Multiple multipole method with automatic multipole setting applied to the simulation of surface plasmons in metallic nanos tructures,” J. Opt. Soc. Am. 19, 101–111 (2002). 40. M. Besbes, J. P. Hugonin, P. Lalanne, S. van Haver, O. T. A. Janssen, A. M. Nugrowati, M. Xu, S. F. Pereira, H. P. Urbach, A. S. van de Nes, P. Bienstman, G. Granet, A. Morea u, S. Helfert, M. Sukharev, T. Seideman, F. Baida, B. Guizal and D. van Labeke, “Numerical analysis of a slit-groove di ffraction problem,” J. Eur. Opt. Soc.2 07022 (2007).


Introduction
High NA optical microscopes are widely used both in scientific research and industrial applications, for example, in the fields of medicine, biology, chemistry and materials science. Since Abbe established the rigorous basis for the analysis of microscopes, low NA systems have been successfully modelled.
For Abbe's theory the first assumption is that scattering and diffraction can be modelled by the scalar approximation. This sets a limit to the NA of the objective lens at ≈ 0.5 [1]. The second assumption is that Fourier optics can be used to model the propagation of light through optical system which in turn invokes the Fresnel approximation and a thin lens model [2]. Much work has been done in this regime, however an exhaustive survey of the literature on the modelling of low NA imaging systems is beyond the scope of this article.
This article is concerned with rigorous vectorial imaging theory. Vectorial theory for calculating fields in the focal region is now well established [3,4,1]. Approximate imaging models have been proposed for restricted scattering objects such as Rayleigh spheres [5] as the scattered field is conveniently calculated analytically. Finite sized spheres have also been modelled using the method of the angular spectrum of plane waves [6] and expansion of the focused beam in terms of Mie modes [7]. Note that the latter paper contains a number of references on other was of describing a tightly focused Gaussian beam scattered by spherical particles. A rigorous vectorial imaging model for general objects has not yet been proposed.
We thus propose a general imaging model of the microscope, composed of four components as depicted in calculation of incident light as described in Section (2). The second component, (blue) is the calculation of the light scattered by the object as discussed in Section (3). The third component (white) is the Stratton-Chu integral, discussed in Section (4), used to resample the field in preparation for the final component of the model. The final element of the model (green) calculates the image of an arbitrary electromagnetic field as discussed in Section (5). Models similar to ours have previously been proposed however each lacks a degree of rigor in one or more of its components. Judkins and Ziolkowski [8] introduced a two-dimensional model which employed a rigorous numerical method, the Finite-Difference Time-Domain (FDTD) method, to model a Gaussian beam scattered by conducting thin film gratings. The scattered field was analysed in the far-field zone after application of a near-to far-field transform based upon the Fraunhoffer approximation. Lui et al. [9] extended the previous work to 3D with the incident beam calculated using an angular spectrum approach. The scattered field was analysed in the far-field zone by applying the Fraunhoffer diffraction integral. Lui et al. [10] made the most recent contribution, calculating the incident light and near-to far-field transform using Fresnel diffraction. The model that we propose here employs a more rigorous method for calculating the incident field and also for calculating the near-to far-field transform. We also use a technique not seen in previous models which allows the scattered field to be calculated at the detector itself, not just in the far-field zone.

Focused illumination
The Debye-Wolf integral is a rigorous solution to a specific electromagnetic boundary value problem which is particularly well suited to the analysis of optical focusing systems. It has been derived in different forms [3,4,1,11]. Although these derivations differ, their results are consistent. An important feature of the Debye-Wolf integral is that it requires the pupil plane to be in infinity which is the case for the vast majority of high NA objective lenses as these are designed to be telecentric from the image side.
Consider such a convergent spherical wave with its centre of convergence at the focus of a high NA lens. Assume the wavefront subtends a solid angle Ω, limited by the aperture of the lens positioned in the first principal focus. Suppose also that the coordinate system in which the focused field will be calculated has its origin at the second principal focus. If the electric and magnetic fields on a particular wavefront are given by E a and H a respectively, the field at a point P described by position vector r p = (x p , y p , z p ) is given by [12]: where s = (s x , s y , s z ) is a unit vector normal to the spherical wavefront and k = 2π/λ is the wavenumber of the incident illumination and f is the focal length of the lens. It can be shown [12] that in order for this to be valid one requires z p ≪ f and that the Fresnel number of the lens must be infinite. An extension of these equations is required on order to model focusing into a stratified medium. This problem has been solved by Török et al. [13] and it is this solution which is employed here. The unit vector of a typical ray propagating through the stratified medium in the Nth medium is denoted s N = (s N x , s Ny , s Nz ) = (cos φ sin θ N , sin φ sin θ N , cos θ N ). The field is to be calculated at position P with position vector r p = (x p , y p , z p ) = (ρ cos ϕ, ρ sin ϕ, z p ). The media are labelled sequentially from left to right. The wave numbers in each medium are labelled similarly as k 1 , k 2 , k 3 , . . . , k N . The interfaces between each medium are positioned at The field in the Nth medium is then given by [14]: (2) where κ = k 1 sin θ 1 cos(φ − φ p ), Ψ i = h N−1 n N s Nz − n 1 h 1 s 1z , Ω 1 is the solid angle of convergence of rays focused by the lens as measured in medium 1 and the substituion ds 1x ds 1y /s 1z = sin θ 1 dθ 1 dφ has been made. In the above formulae E and H are the electric and magnetic strength vectors associated with the geometric optics solution of Maxwell's equations [4,12].

Finding the geometrical optics electric field vector for a focusing lens
The geometrical optics electric field vector is found by vectorial ray tracing. This is most elegantly done using the generalised Jones matrix formalism introduced by Török et al. [6,13] The formalism considers the electric field vector associated with a ray traversing an optical system. Note that what follows is, unless otherwise stated, valid for the calculation of the magnetic field vector as well. For generality a Babinet-Soleil compensator is employed in the first principle focal plane of the objective. If collimated light, linearly polarised in the x-direction is incident normally upon the Babinet-Soleil compensator, the geometrical optics approximation to the field on the converging spherical wavefront is given by: where: where are the transmission coefficients of the stratified medium are defined in [14]. Also, A ± = cos(δ/2) ± i cos(2φ BS ) sin(δ/2), B = sin(δ/2) sin(2φ BS ), φ BS denotes the orientation and δ is the retardation of the Babinet-Soleil compensator. ǫ 0 is the permittivity of free space and µ 0 is the permeability of free space. T is the transpose operator.
Note that in Eqs. (3) and (4), E 0 and H 0 are found for the propagation vector s(φ + π, θ) instead of s(φ, θ) as may be envisaged. This is because when the initial propagation direction (0, 0, 1), of a collimated ray, is transformed in the same way that the geometrical field approximation is, s is obtained as: This is a byproduct of attempting to relate the coordinate systems used to describe coordinates within the first principal focal plane and the directions of rays following the lens, which originated from those positions. Thus, if the polar coordinate φ is used to describe the coordinates within the first principal focal plane, then the propagation vector after refraction by the lens is described by polar angle φ + π as shown in Fig. (2). This has some important implications when the Debye-Wolf integral is finally evaluated. Expanding Eqs. (3) and (4) yields: where we have substituted C ± = T (N−1) s cos θ N . Equations (8) and (9) may be substituted into Eq. (2) respectively to obtain: where (11) where Ψ d = k N z p cos θ N and α 1 is the semi-convergence angle of the objective lens in medium 1.

Interaction of light with specimen
Having calculated the form of the incident light in Section (2), it now remains to calculate how it is scattered by a specimen. There are very few scattering objects for which the scattered field may be calculated analytically. Thus, in a majority of cases a rigorous numerical method must be employed.
The It is beyond the scope of this paper to discuss the strengths and weaknesses of each method however we note that an advantage of our imaging model is that any rigorous method may be used to calculate the light scattered by the specimen.
We employed the FDTD method for several reasons. Firstly, it is easy to implement when compared to the other methods. The simplicity stems from the very logical mapping of the algorithm into a computer program. Furthermore, the algorithm itself is rather simple. Secondly, because the method is very memory efficient it is capable of modelling larger aperiodic scatterers than the MOM and Green's tensor method and at least as big as that possible using the FEM. The FDTD method however does not suffer from the matrix inversion problem which the FEM sometimes has a problem with. The FDTD method is in fact very stable. Thirdly, the FDTD method calculates both the electric and magnetic fields simultaneously with equivalent accuracy. Finally, the FDTD method is a very mature method. It has been used in numerous applications and so a great deal of knowledge has been acquired making the FDTD method a very powerful tool for performing electromagnetic scattering calculations.
The FDTD method does of course have some non-ideal attributes. To start with, it models dispersive materials poorly. This is because it is necessary to take a frequency domain description of a material and convert it to a time domain model. Numerous methods have been used (see Luebbers et al. [19] as an example) successfully however it must still be considered a weak point of the model. We implemented the Drude model using the Auxilliary Differential Equation (ADE) method. Note that the Drude model is incorrect for many metals at frequencies beyond the infrared due to interband absorption [20]. As a result, permitivity is correctly modelled only at the centre wavelength of the simulation.
The regular orthogonal grid employed by the FDTD method is not suited to modelling complex objects. In general, the resolution of object features is limited to the grid cell size. Also, curved and sloping surfaces must be modelled by a stair case approximation. Only the FEM and MOM are able to avoid stair casing as they may employ an irregular mesh. Despite these disadvantages, the FDTD method is probably the most practical and broadly useful method for our purposes. Note that any numerical method could be employed for this component of the imaging model as each component of the imaging model is decoupled from the others. This is a particular strength of the model.
It is beyond the scope of this paper to give a detailed explanation of the FDTD method however a brief introduction is given here. Maxwell's equations result in a set of coupled partial differential equations, one example from this set, for a source free region, is: where σ is the conductivity and ǫ is the permittivity of the region. Most numerical methods require discretisation of field values throughout space. The FDTD method employs a discretisation reported by Yee [21]. The field quantities on the Yee cell are described by an indexing system of the form (i, j, k) which corresponds to a position (i∆x, j∆y, k∆z) where ∆x, ∆y and ∆z are the physical dimensions of the Yee cell.
The field values must also be discretised in time and they are calculated at intervals of a specifically chosen time step ∆t. The electric and magnetic field values are however known half a time step apart. This allows an indexing scheme for time such that a time of index n refers to real time n∆t. Using this system, Yee showed how each partial differential equation of the form of Eq. 12 can be approximated to second order accuracy by a difference equation of the form [22]: and similarly for other components of E and H. Note that α i, j+1/2,k+1/2 and β i, j+1/2,k+1/2 are functions of ∆t and material properties at location (i, j + 1/2, k + 1/2) and that the superscripts n + 1/2, n and n − 1/2 are time indices. The set of difference equations allows an incident field to be introduced to the computational grid and the fields leap frogged in time. The incident field is introduced as a pulse with a Gaussian profile to limit its spectral width. The scattered field at the center wavelength of the pulse may then be found from the time domain data through an application of a discrete Fourier transform. The FDTD algorithm must execute for sufficient iterations so that the scattered field decays to a negligible amplitude. The incident field is introduced at a plane above the sample. This is an approximation however care is taken to ensure that the plane is sufficiently wide and the beam introduced adequately close to its waist such that the incident field is introduced accurately into the computational grid. This is usually described as the "total-field/scattered-field" technique [22]. It would be more accurate to use the "pure scattered field formulation" [22] where the incident field is calculated analytically, everywhere within the computational grid, at each time step. This is however computationally burdensome and so we opted to implement the 'total-field/scatteredfield" technique.
The The primary limitation of the FDTD method is the type of objects which can be accurately modeled. For example, the Yee cell size must be no larger than the smallest feature to be modeled. This can lead to enormous memory requirements when modeling large objects with fine details. Finally we note that in the case where an object is embedded within a stratified medium, which is not included in the FDTD simulation, we assume that it is sufficiently far from an interface such that fields reflected back from the interface are negligible. This is a limitation of the FDTD method not the general imaging model.

Resampling of scattered light
Most rigorous numerical methods for calculating electromagnetic scattering calculate the scattered field on a dense grid. The FDTD method, for example, employs a grid spacing of no larger than λ/20. This means that the scattered field is calculated only in close proximity to the scattering object. This is because too much computer memory is required to include large quantities of homogeneous space within the computational domain. Our objective is, however, to calculate the image of the scattered light and so it is necessary to resample the numerical data rigorously onto a less dense grid with larger lateral dimensions.
Fortunately, this may be done by employing the modified Stratton-Chu integral theorem [27, 12, 28, 29]. The modified Stratton-Chu integral theorem may be stated as follows. Consider a closed surface S 0 containing all sources and sinks of radiation, the exterior of which is composed entirely of homogeneous space. The field anywhere outside the S 0 may be found according to: wherem is an outward surface normal, G is the free-space Green's function, G = exp(ikr)/r and r = r s − r p with r s = (x, y, z) the co-ordinates of an infinitesimal surface element dS on S 0 and r p = x p , y p , z p being the (fixed) observation point outside S 0 . It is reasonable to question how the modified Stratton-Chu integral theorem may be applied when the scattering object is contained within a stratified medium. We are able to do this since we limit ourselves to situations where the scattering object is sufficiently far away from the stratified medium such that light reflected back from the interfaces to the scattering object may be neglected. The propagation of scattered light through a stratified medium is however calculated rigorously in Section (5).
Since the near-field data is known numerically, the Stratton-Chu integral must be evaluated numerically. The surface of integration and associated complex amplitudes are defined using a mesh of triangles. Such a mesh is represented by a set of vertices V = r s,i = r 1 Each element of each triplet in F is an index into V, E and H. In this way, each triangle is constructed from three vertices and the field at each vertex is also known. The orientation of the surface is stored according to the order in which the facet indices are stored. Such a representation minimises computer storage and provides an efficient way to traverse the surface. This representation can be used to represent any polyhedral surface and so is very general.
Integration is performed over each facet and the results summed to give the final result. Gaussian quadrature is commonly used for integration over a triangle. In general, a high order Gaussian quadrature would be employed to improve integration accuracy. However, since the field is known only at the triangle vertices, only first order Gaussian quadrature may be employed. The first order scheme provides an exact result if the integral kernel varies no worse than a polynomial of first order [30]. This is why it is important to use a fine mesh to represent the closed surface of integration.
By employing first order Gaussian quadrature integration, the Stratton-Chu integrals may be evaluated according to: where U is the field of interest, N f acets is the number of facets, I is the kernel of the integral being evaluated, r p = (x p , y p , z p ) is the observation point,m i is the surface normal of facet i, r s,v j i is the jth vertex of facet i, E v j i is the complex electric field at the jth vertex of facet i, H v j i is the complex magnetic field at the jth vertex of facet i, and ∆ i is the area of facet i.

Calculating the image of an arbitrary field
In order to calculate the image of scattered light it is necessary to calculate the geometrical optics approximation to the scattered field on the Gaussian reference sphere. In the absence of a stratified medium this could be done using the modified Stratton-Chu integral theorem, as explained in Section (4), to calculate the scattered field on a section of a sphere with large radius, centered on the focus, subtending the same solid angle as the Gaussian reference sphere of the objective. Having found the geometrical optics approximation of the field on the Gaussian reference sphere, the Debye-Wolf equation could then be applied to find the field at the detector.
The calculation is not as straight forward in the presence of a stratified medium. This is because the modified Stratton-Chu integral theorem cannot be used to propagate the scattered field beyond the stratified medium. Instead, we use an approach grounded in the theories used for calculating the image of a harmonically oscillation electric dipole in a stratified medium [31] and the m-theory [32,33]. It is assumed that the scattered field is known on a plane, S 0 , normal to the optical axis and in the vicinity of the focus of the objective. It is also assumed that this plane is sufficiently large in lateral extent such that the majority of the power radiated by the scatterer into the half space containing S 0 propagates through it. Consider initially the scenario when the stratified medium is not present and S 0 represents an infinite plane separating two half spaces, one containing all sources and sinks of radiation and the other composed of homogeneous space. Then, it can be shown [32,33] that the field at any point P in the far field, situated in the half space which extends to z = ∞, is given by: where Q = (x Q , y Q , z Q ) is a point on S 0 , n is a surface normal directed towards the half space containing P and r is a vector pointing from Q to P. This reveals that each point Q can be considered the source of a spherical wave with polarisation given by r × (n × E(Q)) and this is exactly the type of wave required for use by the Debye-Wolf integral. By applying the theory for imaging dipole waves in stratified media[31] to Eq. (16) it is possible to calculate the image of an arbitrary field known within a stratified medium. Consider first the image of an on-axis equivalent magnetic dipole of moment −k × p embedded in a stratified medium in reflection mode as shown in Fig. (3) a). Then, assuming that the equivalent magnetic dipole has axial position z d p , the field at position r d = (r dt cos ϕ d , r dt sin ϕ, z d ) is given by: where: where we have introduced κ = n d z d cos θ d − n N z d p cos θ N , α d is the semi-convergence angle of the detector lens, θ 1 , θ 2 , . . ., θ N and θ d are shown in Fig. (3) a). The transmission coefficients T (N−1) s,p are as defined in Section (2) except that they are defined for rays propagating from material N into material 1. All other notation is the same as in Section (2). The image of an on-axis equivalent magnetic dipole in transmission, as depicted in Fig. (3) b) is given by the expressions in Eq. (17) however the I eq,s functions are defined as: I eq,s It is then simple to calculate the image of an array of equivalent magnetic dipoles. Suppose that the field due to an on-axis equivalent magnetic dipole of moment n × p in transmission or reflection is given by E oa,t/r (n × p, r d ), then the field at the detector due to an arbitrary field E i defined on a plane S 0 , in the vicinity of the first principal focus of the objective is given by: where q is a vector pointing from the intersection of S 0 and the z-axis to the point Q, the location of a particular equivalent magnetic dipole. In practice the plane S 0 is chosen to be large enough so that "most" of the power scattered by the scattering object propagates through it unperturbed. A degree of approximation is inevitable here. A further approximation is applied to the integral in Eq. (20) as it is evaluated as a sum over discrete equivalent magnetic dipoles. This is necessary for practical implementation.

Evaluation and analysis
We consider the imaging properties of two closely spaced gold spheres in reflection in order to demonstrate the usefulness of our general imaging model. We consider, in particular, how closely two such spheres may be positioned and still be resolved. There are numerous criteria used to assess resolution however we have opted for an adaptation to that of Rayleigh. This adaptation requires the definition of saddle-to-peak ratio. This is defined as the value of the normalised intensity minimum (dip) between the two adjacent intensity maxima due to the two particles. Two scatterers are then said to be just resolved when the saddle-to-peak ratio is approximately 0.735, which matches the saddle-to-peak ratio of the classical definition of Rayleigh's two point resolution criterion. Whereas it is straightforward to calculate the image of a single point scatterer [6], the image of two subwavelength spheres is considerably more complex as multiple scattering must be considered which, in general, necessitates the use of numerical methods. It is in this type of problem that our general imaging model is most useful. We note also that the imaging model is particularly useful for simulating other imaging modes, particularly those which are polarisation sensitive. As an example we refer to a previous publication which shows how the imaging model would be modified to model Nomarski's differential interference microscope [37].
We consider gold spheres of radius 60nm embedded in a medium with refractive index 1.52. A point light source of wavelength 632.8nm is employed at which wavelength gold has a refractive index of 0.2 + 3.32 i [34]. Experimental results [20,8] were used to obtain the Drude model parameters as ω p = 1.114×10 16 sec −1 , ν c = 2.8496×10 14 sec −1 and ǫ ∞ = 2.89. A Yee cell width of λ/40 was employed where λ is the wavelength of light in the embedding medium and the FDTD dimension was 300 Yee cells wide in each lateral dimension and 25 Yee cells deep with a PML 10 cells think around each edge. The Gaussian modulating pulse of the incident field had duration leading to a wavelength width of 60nm. A time step of δt = 1.9042 × 10 −17 sec, or, 0.95 of the Courant stability limit was employed.
All FDTD simulations were performed on our 36-node Beowulf cluster, each node having 1Gb of memory and a 2 GHz processor. The FDTD algorithm was implemented in serial however the cluster permitted fast evaluation of many scattering configurations. The Stratton-Chu integral was implemented in parallel. This results in a significant speedup [29] in calculation of the field on the equivalent magnetic dipole plane.
In order to establish the validity of the model we calculated the image of a single sphere. The image of a small scatterer may be calculated analytically [6] thus allowing the accuracy of our imaging model to be evaluated. Figure 4 shows a line scan of a single sphere along the x-axis  for the case where a focused x-polarised beam is used to illuminate the sphere. The cases of confocal scanning and conventional scanning are plotted. It is evident from the plots that our general imaging model is consistent with the analytic result for the image of a Rayleigh scatterer. The small deviation shown in the image is likely to be caused by inaccuracies stemming from the FDTD method such as the stair-casing approximation used to build the model of a sphere.
Next we consider the image of two spheres under a wide field microscope. The spheres are illuminated by a linearly polarised coherent, monochromatic plane wave due for example to Köhler illumination with sufficiently small field and source apertures, and that the detection optics are the same as in the scanning microscope. Note that in order to model Köhler illumination from an extended coherent source, the source would have to be decomposed into a number of point sources and a separate images, which are later summed, calculated for each such point source.
Assume also that the spheres are positioned on the x-axis. Figure 5 a) shows the two-point resolution calculated using the general imaging model and also using a model where the scatterers are treated as dipoles which do not mutually interact. It is interesting to note that the classical scalar approximation to the two-point resolution limit of .82λ/NA takes the value 0.42 µm. It is evident from the plot that the much simpler dipole model provides a reasonable approximation in the wide-field case. The reason for this will be explained later. By calculating the scattered light for x-and y-polarised focused beams it is possible to calculate the scattered light for linearly polarised incident light of any orientation. The two-point resolution for a particular sphere separation is found by interpolating the saddle-to-peak ratio from the ratios calculated directly using the general imaging model.
Next we consider the two-point resolution under a scanning microscope. Figure  the two-point resolution of gold spheres for incident focused linearly polarised light of differing angle of polarisation. The general imaging model is compared to a simple dipole model where the dipoles are assumed not to interact. Each of these cases is calculated for both confocal and conventional detection modes. Consider first confocal detection. The plots show that there is a correspondence between the dipole and general imaging models however the dipole model differs significantly from the general model. Both two-point resolutions follow the trend of increasing as the incident light tends towards being y-polarised. In the case of conventional detection, both models predict an improvement in the two-point resolution as the incident light becomes y-polarised. The dipole model however predicts a far more significant change than the full imaging model.  shows that there appears to be a particle separation at which all angles of incident polarisation and y-polarised incident light. It is evident from these diagrams that when the scatterers are scanned along the x-axis (indicated by a blue line in the diagrams) they encounter only a y-field component in the case of y-polarised incident light. In the case of x-polarised incident light they encounter both x-and z-field components. It is possible to show mathematically using Eqs. (21) that this is the cause of the variation in two-point resolution with angle of linearly polarised light. We will however only mention that it is the presence of the p z term in Eqs. (21) coupled with the variation of longitudinal field component with angle of linearly polarised light which leads to the variation in two-point resolution.
The difference between the dipole and general models can be explained by considering the interation between the two spheres. Consideration of the field close to the spheres shows that the scattered field differs significantly between the cases of x-and y-polarised incident illumination. Figure 8 shows the first frame of an animation of the scattered electric field intensity in the vicinity of the faces of both spheres. The angle of incident linearly polarised, focused light is varied in the animation and the angle is indicated by a white arrow in the top left hand corner of each frame. The animation shows that for y-polarised incident light, interaction between the spheres leads to a strong electric field between the two spheres.
Detection has hitherto been restricted to ideal confocal (infinitely small pinhole) and ideal conventional (infinitely large pinhole) detection. Previous publications have studied the effect of having a finite sized pinhole in a variety of imaging modes [35, 36, 37]. Figure 9 shows the variation of the two-point resolution with detector radius for three different angle of linearly polarised incident light. As can be seen, the two-point resolution worsens rapidly up until a detector radius of approximately 50µm. Also included in this plot is the full width at half maximum (FWHM) of the microscope point spread function along the x scan direction for xpolarised incident light. This has been calculated using the theory of Török et al. [36] and is included for the sake of comparison. The plot shows that measurement of two-point resolution predicts a limit to confocal behavior which is consistent with that predicted by the point spread function. The model could also be used to simulate a polarised light microscope. One could, for example, position an analyser in the detection path and detect the cross-and co-polarised scattered light. This measurement is however of little interest in this example as the peak ratio of cross-to co-polarised integrated intensity has a peak of approximately 5% for the range of sphere separations and detector radii considered.
The imaging model is also capable of modelling wide field microscopes. In this case the incident illumination is no longer tightly focused. As an example we have calculated the wide field image of a simple pattern etched into a slab of gold as shown in Fig. (10) inset. The pattern was etched into a square slab of width 50λ and each letter had a depth of λ/4. The slab itself had a depth of 2λ. Note that a wavelength of λ=405nm was used along with a 100×, .85 NA objective. The letters were sized such that the central part of the "I" had a width of λ/2. Any wave satisfying Maxwell's equations may be used as the incident illumination, however, we used a plane wave for simplicity. The imaging model is largely unchanged when calculating a wide field image. The only difference is that the plane of equivalent dipoles must be physically larger than in the scanning case due to the larger field of view. Figure 10 shows the first frame of an animation showing the image of the "ICL" structure as the angle of linearly polarised light is varied from horizontal to vertical. This animation shows the importantance of the the vectorial nature of scattering by sub-wavelength structures. Waveguiding appears to occur in the vertical segments of each of the letters as the incident polarisation angle approaches horizontal. The smaller features such as the horizontal part of the "L" are barely visible due to their substantially sub-wavelength width. Polarisation sensitive microscopes are often used to image sub-wavelength structures by looking at the scattered light with polarisation crossed with that of the incident light. Figure  11 a) shows the first frame of an animation which illustrates the evolution of such an image as the angle of the mutually orthogonal incident polarisation and analyser is rotated by π/2. focus of the objective. This is easily modelled by offsetting the axial position of the equivalent dipoles. Such an example demonstrates the ease with which realistic imaging model parameters are incorporated into simulations.

Conclusions
We have shown how realistic imaging systems can be modelled with high accurately. The model permits a wide variety of imaging conditions to be used, including for example, aberrations calculated by lens design software or measured experimentally, or vectorial beams to be used as incident illumination. Arbitrary scattering objects, embedded within a stratified medium, may be modelled using rigorous numerical methods. We have demonstrated the model's practicality through examples.
An additional strength of this method is the way in which the three parts of the imaging system are decoupled, allowing the substitution of alternative techniques to treat one component without affecting the others. For example, the incident illumination considered in this paper is restricted, for convenience, to that which may be obtained using a Babinet-Soleil compensator, or the rigorous numerical block to be substituted by analogous techniques (Finite Element Method, Fourier Modal Method, Volume Integral Method, etc). The method described in this paper can readily be extended to include pulsed laser illumination.