Open Access
24 September 2018 Tool for simulating the focusing of arbitrary vector beams in free-space and stratified media
Peter R. T. Munro
Author Affiliations +
Abstract
Vectorial models of focused beams are important to a variety of fields including microscopy, lithography, optical physics, and biomedical imaging. This has led to many models being developed, which calculate how beams of various profiles are focused both in free space and in the presence of stratified media. The majority of existing models begin with a vectorial diffraction formula, often referred to as the Debye–Wolf integral, which must be evaluated partially analytically and partially numerically. The complexity of both the analytic and numerical evaluations increases significantly when exotic beams are modeled, or, a stratified medium is located in the focal region. However, modern-day computing resources permit this integral to be evaluated entirely numerically for most applications. This allows for the development of a vectorial model of focusing in which the focusing itself, interaction with a stratified medium, and incident beam specification are independent, allowing for a model of unprecedented flexibility. We outline the theory upon which this model is developed and show examples of how the model can be used in applications including optical coherence tomography, high numerical aperture microscopy, and the properties of cylindrical vector beams. We have made the computer code freely available.

1.

Introduction

The calculation of electromagnetic fields in focal regions has been the subject of an immense amount of work over, at minimum, the last century. A formalism allowing the solution of this problem has been developed by Ignatowsky,1 Luneburg,2 and Richards and Wolf.3 This formalism has been applied to a variety of focusing problems, from studying the structure of tightly focused fields obtained by focusing linearly polarized light3 to more complex applications such as focusing Gauss–Laguerre4 and azimuthally and radially polarized5 beams, for example. The formalism has also been applied to focusing through stratified media (see, for example68). This is an important extension to focusing theory due to the large number of applications where this is relevant, such as image formation in high numerical aperture confocal microscopes,9 optical coherence tomography (OCT),10 and photoacoustic imaging,11 to name just a few.

Many contributions have been made to the problem of focusing vector beams through stratified media (see, for example7,8,1215). A barrier to researchers in the field putting these solutions to use is their implementation in computer software. Although the integral equations are more detailed than complex, their implementation may appear daunting to researchers from outside the field. Furthermore, even slight modifications to the beam type require a good understanding of the underlying formalism to correctly modify the equations governing the form of the focused field. One of the reasons for this is that the computational complexity of calculating the form of focused fields has hitherto dictated that as much as possible of the computation be performed analytically. The result of this is that it can be difficult for nonspecialists to modify the theory for their application and the final equations can lose connection with the starting point of the focusing formalism, which can be understood conceptually.

Calculating the field at each point in the focal region requires independent evaluation of a two-dimensional (2-D) diffraction integral. Early calculations of fields focused through a stratified medium employed an Intel Pentium processor operating at a frequency of 90 MHz.12 Numerical solution of the integrals presented in that work took advantage of the rotational symmetry of the optical system. In particular, integration about the axis of revolution of the optical system (i.e., the optical axis) was performed analytically, leaving only one-dimensional (1-D) integrals to be numerically evaluated. Aside from the elegance of the analytic solution, this was also necessary to keep the computational cost of evaluating the focused fields feasible. The downside of using this semianalytic approach is that the process of making the required coordinate system transformations and manipulations results in a series of equations that are not readily physically interpreted and are not easily modifiable to treat different beam types or optical systems.

The results contained within this paper were calculated on an Intel Xeon processor with 12 cores operating at a frequency of 2.3 GHz, thus reflecting a significant increase in processing power over the past two decades. This increase in computational power makes it feasible to solve the problem of tight focusing in the presence of a stratified medium largely numerically. This carries several advantages. First, it enables a more direct mapping between the equations governing the solution and the numerical implementation, making the solution easier to extend to alternative beam types and optical systems. Second, the focusing problem can be considered as a vectorial angular spectrum of plane waves, thus allowing the focusing and stratified media problems to be naturally decoupled. This provides for a very convenient way of calculating the electromagnetic field throughout a stratified medium, not just the transmitted field as has previously been the case.6 The aim of this work is to provide a formalism, along with freely available computer software,16 which is able to be used, with little prior knowledge, to calculate a range of focusing problems, possibly involving a stratified medium. Furthermore, the formalism is structured in a manner that prioritizes conceptual understanding of the code rather than computational efficiency. This allows the formalism to be extended by users with an understanding of the underlying theory and not necessarily a detailed knowledge of scientific computing, which is why we have presented this work in the form of a tutorial. Furthermore, it is our intention that this tutorial, along with the associated computer code, will enable researchers from a broad range of fields to easily solve problems related to high numerical aperture focusing.

In the rest of this paper, we introduce the mathematical basis of our method for calculating the focused field in the vicinity of a stratified medium that is adapted from that previously proposed by van de Nes et al.7 We begin with the Debye–Wolf integral before considering the matrix formalism for treating the propagation of plane waves within a stratified medium. We then show how these two formalisms can be combined before considering the numerical implementation of the resulting formalism. We illustrate the use of the model through three examples: image formation in OCT, the focusing of cylindrical vector beams through a dielectric interface, and a simulation of aberrations typically experienced in confocal microscopy. We note that the stratified media focusing code, along with code to generate the results presented in this paper, is freely available for download.16

2.

Mathematical Foundation

Throughout this paper, we represent the electric field vector at position r, as a time-harmonic quantity, E(r), such that the time-dependent field is given as

Eq. (1)

E˜(r,t)=R{E(r)exp(iωt)},
where ω is the angular frequency of time-harmonic oscillation and R refers to the real part. We shall also use the convention whereby a three-dimensional (3-D) vector, for example, a=(a1,a2,a3) can be written as a=(a,a3), where a=(a1,a2). This allows, for example, for the transverse and axial components in Cartesian space to be conveniently expressed independently. We shall denote unit vectors with a hat, i.e., a^, and where unit vectors are plotted (such as in Fig. 1), it is understood that no attention is paid to scale.

Fig. 1

Diagrams illustrating the coordinate systems used in the evaluation of the final form of the Debye–Wolf integral [Eq. (21)], upon which this paper is based. In particular: (a) illustrates the coordinate system used within the lens aperture and the focal region, (b) is a cross-sectional diagram corresponding to the shaded meridional plane in (a), (c) shows the coordinate system within the lens aperture, and (d) is the coordinate system in the focal region.

JBO_23_9_090801_f001.png

2.1.

Debye–Wolf Integral

The field at position (r,z), in the vicinity of the focus of a lens of arbitrary numerical aperture (NA) and focal length f, which satisfies the sine condition, can be found by evaluating the Debye–Wolf integral13 as

Eq. (2)

E(r,z)=in0k0f2π|s|<NA/n0ϵ(s)exp[in0k0(s·r+1|s|2z)]d2s1|s|2.
The axial position [i.e., z in Eq. (2) and Fig. 1] of the point of observation is measured from the location of the nominal focus. The lateral position of the point of observation [i.e., r in Eq. (2) and Fig. 1] is measured from the optical axis. The point of observation is assumed to be located within a homogeneous region of refractive index n0. As shown in Fig. 1, upon multiplication by f, the 2-D vector s uniquely defines the position within the aperture from where a typical ray, assumed to propagate parallel to the optical axis, originates before being refracted by the lens. After refraction, the ray propagates parallel to (s,1|s|2). The field on a Gaussian reference sphere located in the exit pupil of the configuration shown in Fig. 1 is denoted as ϵ(s). The Gaussian reference sphere is shown in Fig. 1 and we note explicitly that ϵ(s) corresponds to the field at the location on the reference sphere intersected by a line parallel to (s,1|s|2) and intersecting the nominal focus. It is important to note that as the aperture is located in the back focal plane of the lens, the exit pupil appears at infinity as viewed from the space containing the focal region. It can be shown17 that this means that ϵ(s) can, in fact, be found using geometrical optics.

The first task in evaluating Eq. (2) is thus to find, using geometrical optics, ϵ(s) for an arbitrary field in the back focal plane, which we denote by ϵfp(fs), where the superscript ϵfp stands for the focal plane. It is assumed that ϵfp(fs) has no axial component. We consider a ray propagating in the sample space, parallel to unit vector (s,1|s|2). We define θ as the angle that this ray makes with the optical axis, as shown in Fig. 1(a). This ray propagates in a particular meridional plane that makes an angle ϕ to the x axis in the sample space as shown in Fig. 1(c). Suppose ϵfp(fs)=ϵxfp(fs)i^+ϵyfp(fs)j^, where i^ and j^ are the unit vectors in the x- and y-directions, respectively. The generalized Jones matrix formalism18 provides a systematic way to find ϵ(s) from ϵfp(fs). We start by defining a coordinate system transformation matrix R, which transforms from a global Cartesian coordinate system to one where the x and y axes have been rotated by an angle ϕ around the z axis. We also define a lens rotation matrix L that rotates a vector by an angle of θ about the y axis of the coordinate system. These matrix transformations are defined as18

Eq. (3)

R=[cosϕsinϕ0sinϕcosϕ0001]L=[cosθ0sinθ010sinθ0cosθ].
It is then straightforward to show that

Eq. (4)

R1(ϕ+π)L(θ)R(ϕ+π)[001]=[cosϕsinθsinϕsinθcosθ]=[s1|s|2],
and furthermore that18

Eq. (5)

ϵ(s)=cosθR1(ϕ+π)L(θ)R(ϕ+π)[ϵxfp(fs)ϵyfp(fs)0]=12cosθ[1+cosθ+cos2ϕ(cosθ1)sin2ϕ(cosθ1)sin2ϕ(cosθ1)1+cosθcos2ϕ(cosθ1)2cosϕsinθ2sinϕsinθ][ϵxfp(fs)ϵyfp(fs)].
We emphasize that when we write an equation involving one or more of ϕ, θ, and s, it is implied that the values θ and ϕ are determined by s such that Eq. (4) is satisfied.

2.2.

Stratified Medium Matrix Formalism for Plane Waves

As will be explained in further detail in this section, the Debye–Wolf integral expresses the electromagnetic focusing problem as a superposition of plane waves. In this section, we thus consider the interaction of plane waves with stratified media.

2.2.1.

Preliminary considerations

We follow a stratified medium matrix formalism as proposed by Azzam and Bashara.19 We consider a stratified medium as shown in Fig. 2 with interfaces located at z=h1,h2,,hN, therefore possessing up to N+1 different media having refractive indices n0,n1,,nN, the first and last layers being infinite half-spaces. We allow each refractive index potentially to be complex valued as n=n+in and assume that both n and n are positive. We assume that the first medium is lossless and so n0 is purely real valued. We decompose the electric field vector of waves propagating in the positive z-directions onto transverse electric (TE) and transverse magnetic (TM) components with unit vectors a^TE+ and a^TM+, respectively. We follow a similar convention for waves propagating in the negative z-direction, except we interchange the superscript + with .

Fig. 2

The assumed structure of the stratified medium and associated notation as used in the evaluation of the final form of the Debye–Wolf integral [Eq. (21)] and its constituent terms.

JBO_23_9_090801_f002.png

2.2.2.

Incident, reflected, and transmitted fields

We describe a plane wave, incident from the objective lens, propagating in the first medium as

Eq. (6)

E+(s;r,z)=ϵ(s)exp(in0k0s·r)1|s|2exp(in0k0z),

Eq. (7)

=[aTE+(z)a^TE++aTM+(z)a^TM+]exp(in0k0s·r),
where we decompose ϵ(s)exp(in0k0z)/1|s|2 into TE and TM components denoted by aTE+(z)a^TE+ and aTM+(z)a^TM+, respectively. Following the convention of Azzam and Bashara,19 we define these vectors as

Eq. (8)

a^TE+=k^×(s,1|s|2)a^TM+=a^TE+×(s,1|s|2)a^TE=a^TE+a^TM=a^TE×(s,1|s|2),
where k^ is the unit vector parallel with the optical axis. Without loss of generality, we henceforth assume that each plane wave component is either TE or TM polarized, thus allowing us to drop the TE and TM subscripts from Eq. (7). Consider such a plane wave incident upon the stratified medium shown in Fig. 2. The wave incident upon the first interface, at the first interface, takes the form E+(s;r,h1)=a+(h1)a^+exp(in0k0s·r), where we have added the superscript to h1 to indicate that this is the limiting value of the field as z approaches h1 from the negative z-direction (see Fig. 2). Note that the z-dependence of E+(s;r,h1) is expressed entirely by the scalar a+(h1). Following the matrix formalism presented by Azzam and Bashara,19 we start by defining E(s;r,h1), E+(s;r,h1+), and E(s;r,h1+) as shown in Fig. 2, meaning that we have defined plane waves propagating in the positive and negative z-directions, evaluated on both sides of the interface at z=h1. We can then define a matrix relationship between a+(h1), a(h1), a+(h1+), and a(h1+) as19

Eq. (9)

[a+(h1)a(h1)]=F01[a+(h1+)a(h1+)]=1t01[1r01r011][a+(h1+)a(h1+)],
where F01 is implicitly defined and t01 and r01 are the Fresnel transmission and reflection coefficients, respectively, for the interface at z=h1. Next, we consider relating the fields within the second medium using a similar matrix formalism. To do this, consider plane wave components at the beginning and end of the second medium, i.e., E+(s;r,h1+), E(s;r,h1+), E+(s;r,h2), and E(s;r,h2). We can then define a relationship between their plane wave amplitudes, a+(h1+), a(h1+), a+(h2), and a(h2) as19

Eq. (10)

[a+(h1+)a(h1+)]=P1[a+(h2)a(h2)]=[exp(iβ1)00exp(iβ1)][a+(h2)a(h2)],
where β1=k0(h2h1)n11|(n0/n1)s|2 and P1 is implicitly defined. It is then possible to construct a matrix that relates the fields at z=h1 and z=hN as

Eq. (11)

[a+(h1)a(h1)]=S(h1;hN+)[a+(hN+)a(hN+)],
where S(h1;hN+) is obtained by successive application of Eqs. (9) and (10) as

Eq. (12)

S(h1;hN+)=[s11(h1;hN+)s12(h1;hN+)s21(h1;hN+)s22(h1;hN+)]=F01P1F12P(N2)F(N2)(N1)P(N1)F(N1)N.
We will assume in this work that a(hN+)=0, meaning that all light reaches the stratified medium after transmission through the objective lens, thus giving

Eq. (13)

a+(hN+)=1s11(h1;hN+)a+(h1)a(h1)=s21(h1;hN+)s11(h1;hN+)a+(h1).
The only thing remaining is to find the vectors a^0 and a^N+, the polarization vectors in media 0 and N, respectively. This is trivial in the TE case as they remain fixed in all media. In the TM case, the transverse component of the wave propagation vector in each medium (i.e., s,0 and s,N) is first calculated using the Snell’s law (see Sec. 2.2.5). The polarization vectors may then be found according to Eq. (8). The reflected field at the first interface and transmitted field at the last interface can thus be found as a(h1)a^0exp(in0k0s·r) and a+(hN+)a^N+exp(in0k0s·r), respectively. Note that the dependence upon z is entirely embedded within the matrix formalism. In general, this step must be performed for the TE and TM cases independently and the results superimposed to obtain the total field. We have thus shown how to calculate the reflected and transmitted fields at the first and last interfaces, respectively, for an arbitrary stratified medium and an incident plane wave.

2.2.3.

Internal fields at internal interfaces

Suppose it is desired to know the field at an interface in the interior of the stratified medium, say, for example, the field at hq, for some positive integer q. As in Sec. 2.2.2, we can construct matrix relationships relating the field components at hq to the incident, reflected, and transmitted fields as

Eq. (14)

[a+(hq)a(hq)]=S(hq;hN+)[a+(hN+)a(hN+)],

Eq. (15)

[a+(h1)a(h1)]=S(h1;hq)[a+(hq)a(hq)],
where

Eq. (16)

S(hq;hN+)=F(q1)qPqFq(q+1)PN2F(N2)(N1)PN1F(N1)N
and

Eq. (17)

S(h1;hq)=F01P1F12P(q2)F(q2)(q1)P(q1).
This means that Eq. (14) can be used in conjunction with Eq. (13) to yield

Eq. (18)

a+(hq)=s11(hq;hN+)s11(h1;hN+)a+(h1)a(hq)=s21(hq;hN+)s11(h1;hN+)a+(h1),
or Eq. (15) can be used in conjunction with Eq. (13) to yield

Eq. (19)

[a+(hq)a(hq)]=S1(h1;hq)[1s21(h1;hN+)s11(h1;hN+)]a+(h1).

2.2.4.

Calculating fields between interfaces

There are many instances where it is necessary to evaluate the field away from an interface, either internal or external. One approach to this would be to first calculate the field due to the plane wave in question at an interface adjacent to the plane of constant z in which the field must be evaluated. The field at the desired location can then be propagated from the adjacent interface using the plane wave’s propagation vector obtained using Snell’s law. There is a certain numerical efficiency to this approach, however, it requires keeping track of propagation vectors and plane wave phase terms. An alternative approach that is consistent with the matrix formalism discussed above is to introduce a “phantom” interface, which is an interface between two layers with identical refractive indices. In particular, suppose we wish to evaluate the field at some location z=hobs that lies in the layer having refractive index nq. The field at z=hobs can be found by following a similar procedure to that outlined in Sec. 2.2.3, except that we introduce an additional interface at hobs. Then, we evaluate S(hobs;hN+) and S(h1;hobs) instead of S(hq;hN+) and S(h1;hq), respectively, and apply either Eq. (18) or Eq. (19) to evaluate a+(hobs) and a(hobs).

2.2.5.

Conducting media and total internal reflection

In the case where media have a complex-valued refractive index, or when total internal reflection occurs, the wave vector of a plane wave becomes complex valued. This does not pose any difficulties as Snell’s law and the Fresnel reflection and transmission equations remain formally valid (see, for example1921). Snell’s law can be generalized in the notation of this work by defining s,i as the transverse component of a plane wave’s propagation vector in the i’th medium. The propagation vector in this medium will thus be given by (s,i,1|s,i|2). We can then write Snell’s law as

Eq. (20)

n0s=nis,i,
where it is assumed that, as in Sec. 2.2.2, s is the transverse component of the incident plane wave in a homogeneous medium of refractive index n0. The axial component of the plane wave in the i’th medium may then be found as 1|s,i|2.

2.3.

Stratified Medium Matrix Formalism for Focused Illumination

We have shown in Sec. 2.1 how focusing of electromagnetic waves can be described as a superposition of plane waves using the Debye–Wolf integral. We have also shown in Sec. 2.2 how the field at any location within or outside a stratified medium can be evaluated for an arbitrary incident plane wave. Thus, by virtue of the linearity of the system, the forward or back propagating field anywhere within, or outside, the stratified medium may be found according to

Eq. (21)

E±(r,z)=in0k0f2π|s|<NA/n0[aTE±(z)a^TE±+aTM±(z)a^TM±]exp(in0k0s·r)d2s,
where aTE±(z) and aTM±(z) are evaluated using Eqs. (13), (18), or (19) as appropriate, depending on the value of z. We find a^TE± and a^TM± using Eq. (8) combined with knowledge of s in the appropriate medium as obtained from Eq. (20).

2.4.

Numerical Implementation

Equation (21) must be evaluated numerically. This is done by changing the variables of integration from s to (θ,ϕ) as defined by Eq. (4). Evaluation of the Jacobian of this transformation results in d2s=sinθcosθdθdϕ. The resulting integral is then evaluated using Gauss–Legendre numerical integration. The Gauss–Legendre technique allows for a definite integral to be evaluated as22

Eq. (22)

abf(x)dxi=1Nwif(xi),
where xi are the zeros of the N’th-order Legendre polynomial, wi are the associated weights, and N is the number of points used to evaluate the integration. A detailed explanation regarding the calculation of both the xi and wi is given by Davis et al.22 The 2-D integration of Eq. (21) is evaluated by integrating into θ and ϕ separately.

3.

Examples

3.1.

A-Scan Formation in Optical Coherence Tomography

The inverse scattering problem for OCT has been addressed previously in a 1-D setting using the theory of stratified media for plane waves.23 Although we do not consider the inverse scattering problem here, rather the forward problem, we demonstrate how selected problems with practical significance can be addressed using the model developed in this paper. We consider two stratified media as shown in Fig. 3. The stratified medium in Fig. 3(a) is used to demonstrate how the so-called confocal function of an OCT system can be calculated and that of Fig. 3(b) shows how dispersion of water changes the axial resolution of an OCT system. The simulation presented here is idealized and valid for both spectral domain and swept source OCT systems in the sense that a range of wavenumbers spanning λ=1180  nm to λ=1420  nm is uniformly sampled. We consider a fiber-based system with an objective having NA=0.097. The resulting spot size of the beam is determined by the combination of collimator, objective lens and mode field diameter of the single-mode fiber. This is integrated into our model by specifying the functions Exfp(fs) and Eyfp(fs) as

Eq. (23)

ϵxfp(fs)=exp[(F|s|/NA)2]ϵyfp(fs)=0,
where F=1.556 is derived to be consistent with a previously characterized OCT system.24 The A-scan formation model is based upon a recently published full-wave model of OCT image formation.24 In brief, for both stratified media in Fig. 3, the incident beam that would be observed in the absence of a stratified medium is calculated just in front of the first refractive index discontinuity yielding E+(r,0). The field scattered back by the stratified medium is also calculated yielding E(r,0). The modal coefficient describing the scattered field coupled back in to the fiber is given as25

Eq. (24)

αsc=[E+(r,0)·i^][E(r,0)·i^]d2r,
a similar quantity, αref representing the field coupled back into the fiber from the reference arm is evaluated by repeating the above calculation for a single refractive index interface between air and a nonphysical material having refractive index 104. The OCT A-scan is then evaluated as

Eq. (25)

I(zzref)=0S(k)|αsc(k)+αref(k)|2exp[ik2(zzref)]d(1/λ),
where S(k) represents the effective shape of the OCT system’s spectrum. As already mentioned, full details of this model are available in a previous publication.24

Fig. 3

Two stratified media used for calculating example A-scans.

JBO_23_9_090801_f003.png

The first stratified medium considered shows how the confocal function of an OCT system can be calculated using a single A-scan. A refractive index distribution was modeled as shown in Fig. 3(a), where slabs, 100-nm wide, of material with refractive index of n1=1.01 are spaced every 25  μm in air. The narrow slabs act as weak reflectors that do not significantly perturb the illumination beam. The magnitude of the A-scan is plotted in Fig. 4(a), which shows a peak (plotted in blue) for each refractive index discontinuity in the stratified medium. These peaks lie on an envelope, plotted in red, which is evaluated as

Eq. (26)

Iconf(z)=[Efs+(r,z)·i^]2d2r,
where the subscript fs denotes that Efs+(r,z) is evaluated in free space. This example does not produce significant multiple reflections due to the low refractive index contrast between the reflecting slabs (n1=1.01) and the background refractive index (n0=1).

Fig. 4

(a) The magnitude of an A-scan obtained from the refractive index distribution in Fig. 3(a) (blue) and the directly evaluated confocal function (red). (b) The magnitude of an A-scan, on a logarithmic scale, for the refractive index distribution shown in Fig. 3(b). (c) The three peaks in Fig. 3(b) aligned and normalized to the same peak value, demonstrating the effect of dispersion in water.

JBO_23_9_090801_f004.png

The second example of A-scan simulation involves a 500-μm layer of water sandwiched between layers of glass as shown in Fig. 3(b). For this example, we assume that the glass layers are infinitely thick. The refractive index of the glass was assumed to be 1.525 at all wavelengths within the spectrum of the OCT system. The refractive index of water was found by interpolating the data published by Segelstein.26 An A-scan is plotted in Fig. 4(b) on a logarithmic vertical axis. The horizontal axis is scaled by the group refractive index of water obtained by numerically evaluating ng=npλdnp/dλ at the central wavelength of the OCT system’s spectrum, where np is the phase refractive index. The A-scan of Fig. 4(b) shows three peaks separated by a physical distance of 500  μm. The peak at z=0 corresponds to the reflection from the interface between glass and water at z=0 in Fig. 3(b). Similarly, the peak at z=500  μm is from the interface between water and glass at z=500  μm in Fig. 3(b). The third peak is due to light that has been multiply reflected and has undergone two complete round-trips of the water layer. The impact of dispersion due to water on the axial response of an OCT system can be studied by aligning each of these three peaks upon the same axial location as is shown in Fig. 4(c). These plots demonstrate that the peak corresponding to the first reflection does not exhibit dispersion as the A-scan is symmetric about z=0. The two subsequent peaks that correspond to light propagating greater distances in water show how dispersion broadens axial response of an OCT system.

3.2.

Cylindrical Vector Beams Focusing Through a Dielectric Interface

There are many examples in the literature of where the properties of unconventional beams focused through dielectric interfaces are studied (e.g., Refs. 4,7, and 2729). In these, and many other related papers, much effort is devoted to evaluating Eq. (2) partially analytically, for a nonstandard form of ϵ(s) such as may result when focusing, for example, a Bessel–Gauss beam. In such papers, ϵ(s) is evaluated for a specific beam in the aperture of the focusing lens, ϵfp(fs), using the generalized Jones matrix formalism18 outlined in Sec. 2.1. In the absence of a stratified medium, the resulting value of ϵ(s) is substituted into Eq. (2) before changing the variables of integration from s to (θ,ϕ). The integration over ϕ is performed analytically using the identities 02π[sin|cos](nϕ)exp[ixcos(ϕφ)]dϕ=2π(in)Jn(x)[sin|cos](nφ) for some integer n, where Jn is the n’th order Bessel function of the first kind. This results in an expression for the complex amplitude of the electric field in the vicinity of the focus in terms of the so-called “I”-functions, which are functions of (r,z) that are evaluated by performing a numerical integration for each unique combination of |r| and z. The significance of this procedure is that when a new beam type is to be modeled, it is generally necessary to derive the equations that determine the field from scratch. This process is made more difficult when a stratified medium is considered.

In this example, we demonstrate first how a previously published result can be obtained, using the formalism presented in this paper, without the need to perform any analytic derivations. In particular, we consider the focusing of a radially polarized Bessel–Gauss beam originally considered by Biss and Brown27 who calculated how a beam of the form

Eq. (27)

ϵfp(fs)=E0[s^0]exp(|s|2max(|s|2))J1(2|s|max(|s|)),
incident upon the aperture of a focusing lens, interacts with an air interface located at the focus. The lens was assumed to have a numerical aperture of 1.4 when focused into a medium of refractive index n0=1.518.

We consider three different scenarios: focusing into a homogeneous medium with refractive index n0, through an interface at the focus between refractive indices n0 and n1=1 and, finally, through a λ/4 wide layer of refractive index n1 embedded in a refractive index n0. Scenario 1, depicted in Figs. 5(a) and 5(d), shows the profile of the beam in a homogeneous medium. The second scenario [Figs. 5(b) and 5(e)] replicates what was published by van de Nes et al.30 The final result [Figs. 5(c) and 5(f)] that we show as an example is where, instead of having an interface between n0 and n1, there is a λ/4 wide layer of refractive index n1. The subwavelength width of the layer allows for the evanescent waves to couple into the n0 region causing the field after the final interface to more closely resemble that of the homogeneous case. This example is intended to demonstrate how complex focusing problems can be studied using this formalism in a much simpler way than was previously possible.

Fig. 5

(a) and (d) Intensity distributions for the Bessel–Gauss beam described in Eq. (27) focused into free space with refractive index n0=1.518, (b) and (e) through an interface between media with refractive indices n0 and n1=1, and (c) and (f) through a λ/4 wide air gap of refractive index n1 surrounded by a medium of refractive index n0. Subfigures (a)–(c) show |Ez|2 and (d)–(e) show |Ex|2+|Ey|2. Plots have a common colormap, where 0 dB corresponds to the peak value across all six plots.

JBO_23_9_090801_f005.png

3.3.

Simulation of Aberrations Arising in Confocal Microscopy

For the final demonstration of how this formalism can be used, we consider an example from the literature of where a similar model has been applied.12 In particular, we study the illumination point spread function (PSF) that arises when performing confocal microscopy with a dry objective, with a sample embedded in water, covered by a glass cover slip as shown in Fig. 6. This particular example, considered previously by Török and Varga,12 employs a dry objective of numerical aperture 0.9, a 170-μm cover slip of refractive index n1=1.54 and a probe fluorescent molecule situated in water 50  μm below the glass/water interface. The fluorescent molecule is merely a way of depicting the location at which we evaluate the focused field, we do not consider the process of fluorescence.

Fig. 6

Schematic diagrams of the confocal microscopy aberration example (a) and (b) are identical except for relative shifts of the stratified medium and probe fluorescent molecule. Light is focused from air into a glass cover slip of thickness 170  μm and then finally into water. A probe fluorescent molecule, located 50  μm from the glass/water interface, is depicted in green, where the intensity of the incident light is recorded.

JBO_23_9_090801_f006.png

When focusing from air through a glass cover slip, spherical aberration will occur if aberration compensation is not implemented. Although an analytic approach was derived by Török and Varga,12 we will use the stratified medium matrix formalism employed in this paper for this purpose. Suppose that we have an arbitrarily complex stratified medium that we wish to compensate for. We must calculate the phase aberration suffered by each plane wave components, aTE±(z) and aTM±(z), appearing in Eq. (21), due to the stratified medium. This phase aberration can be found, as a function of s, using the matrix S(h1;hN+) that describes the entire stratified medium in Eq. (12). In particular, we know the ratio between the amplitude of the plane wave incident upon the stratified medium [a+(h1)] and that of the transmitted wave [a+(hN+)] is given by Eq. (13) as

Eq. (28)

a+(hN+)a+(h1)=1s11(h1;hN+).
We then use the matrix formalism again to calculate the ratio of amplitudes for the same incident plane wave when all refractive indices are set to that of air, i.e., the aberration free case. We denote this ratio 1/s11fs. We thus define the aberration caused by the stratified medium by exp(iWT) as

Eq. (29)

WT=arg[s11fs/s11(h1;hN+)],
where arg represents the argument of a complex number. Although s11fs/s11(h1;hN+) is often real, it may be a complex value if evanescent waves are produced, which cannot be compensated for.

We used this approach to calculate the rotationally symmetric aberration produced by the glass cover slip as plotted in Fig. 7(b). Another aspect of the model presented in this paper is the ability to model aberration functions with ease. In particular, it is straightforward to incorporate Zernike fringe polynomials into the model. For this work, we consider only the rotationally symmetric Zernike fringe polynomials defined as

Eq. (30)

Z1=1,

Eq. (31)

Z4=1+2ρ2,

Eq. (32)

Z9=16ρ2+6ρ4,

Eq. (33)

Z16=1+12ρ230ρ4+20ρ6,

Eq. (34)

Z25=120ρ2+90ρ4140ρ6+70ρ8,

Eq. (35)

Z36=1+30ρ2210ρ4+560ρ6630ρ8+252ρ10,
where ρ=|s|/max(|s|). As these polynomials are orthogonal on the unit circle, we can decompose the aberration WT onto these functions using numerical integration. The result of this is shown in Fig. 7(a), which shows the contribution due to each individual polynomial. However, since we consider only the first six rotationally symmetric polynomials, aberration compensation described by a sum of these polynomials cannot match, exactly, the actual aberration. This error is plotted in Fig. 7(c) for the case where between three and six polynomials are considered.

Fig. 7

(a) Contributions of each Zernike fringe polynomial to the aberration due to the glass cover slip, (b) the total phase aberration, WT, due to the glass cover slip, and (c) the phase error that results when between three and six Zernike fringe polynomials are used to represent WT.

JBO_23_9_090801_f007.png

Now that the aberration introduced by the glass cover slip has been calculated, it can be compensated for by applying a phase modulation to the beam incident upon the aperture of the focusing lens as

Eq. (36)

ϵxfp(fs)=exp(iΨ)ϵyfp(fs)=0,
where, in this example, Ψ will be evaluated in two different ways. One way is by directly setting Ψ=WT. The alternative way is by setting Ψ=i=1Nci2Zi2, for values of N between 3 and 6, where ci2 is the weight associated with polynomial Zi2 evaluated using numerical integration as ci2=01Zi2(ρ)WT(ρ)ρdρ/01Zi2(ρ)Zi2(ρ)ρdρ.

The unaberrated illumination PSF, evaluated as if the probe fluorescent molecule was in air, is plotted in Fig. 8(f). All other PSFs in Fig. 8 are normalized by the peak of this PSF. The PSFs in Figs. 8(a)8(e) are calculated with the combination of glass cover slip and fluorescent molecule embedded in water. Aberration compensation is applied using an increasing number of Zernike fringe polynomials beginning with three in Fig. 8(a) (i.e., Ψ=i=13ci2Zi2) and going up to six in Fig. 8(d) (i.e., Ψ=i=16ci2Zi2). Significant spherical aberration is present in each of cases (a), (c), and (d) due to the fluorescent molecule being located in water rather than air. It is interesting, however, that case (b), where four Zernike fringe polynomials are considered, is substantially less aberrated. Analysis shows that the combination of aberration due to compensation, combined with propagation through the cover slip and water, results in less overall aberration in the water. Although not considered in this example, we emphasize that noncircularly symmetric aberration polynomials can be just as easily treated within this model.

Fig. 8

(a–e) The intensity of the focused illumination at the location of the fluorescent molecule as a function of shift between the glass/water interface for a variety of aberration compensations and (f) for focusing directly into air without aberration. Intensities are normalized by the peak value in air. Plots (a), (b), (c), and (d) correspond to compensation by three, four, five, and six Zernike fringe polynomials, respectively. Plot (e) corresponds to compensation by the directly evaluated phase aberration.

JBO_23_9_090801_f008.png

4.

Discussion and Conclusions

The model presented in this paper will be useful across many fields of optics, particularly in optical imaging. The strength of the model is that its components that deal with the beam type, stratified medium and focusing are decoupled from one another. This has some important ramifications. One is that this means that incident beams and stratified media of arbitrary complexity can be modeled just as simply as a free-space focusing problem. Another consequence is that the implementation in software closely mirrors the theory itself, making the model amenable to use and modification by others. This approach is quite deliberate as it is intended that this model be available for use by researchers in the field, freeing them from the need to implement their own bespoke focusing code when a new, slightly different, application arises.

Examples including OCT, cylindrical vector beams, and high numerical aperture focusing have been considered. Each of these examples is based upon a result that has already been published, plus some new results. These examples have been chosen to demonstrate the breadth of applications, which this model is able to address.

There are several extensions that we intend to make to this model. Perhaps the most straightforward is the calculation of the magnetic fields, which is trivially evaluated as the magnetic field associated with a plane wave can be directly evaluated from the electric field amplitude and polarization, propagation vector and the refractive index. Consideration of a tilted stratified medium is also possible and of practical significance. Modeling of anisotropic layers within the stratified medium would also be of high practical importance. There are many improvements that can be made to the computational efficiency of the code. Little has been done to date in this regard as the emphasis of the code was on readability. However, there is scope for making the existing code more efficient. Furthermore, as the code is written in MATLAB it should be straightforward to execute it on a graphics processing unit using MATLAB’s parallel computing toolbox.

Disclosures

The author does not declare any financial interests or conflicts of interest.

Acknowledgments

P.M. was supported by a Royal Society University Research Fellowship. P.M. acknowledges fruitful discussions with Professor Peter Török.

References

1. 

V. S. Ignatowsky, “Diffraction by a lens of arbitrary aperture,” Trans. Optical Inst. Petrograd, 1 (4), 1 –36 (1919). Google Scholar

2. 

R. Luneburg, Mathematical Theory of Optics, University of California Press, Berkeley and Los Angeles (1966). Google Scholar

3. 

B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems ii. Structure of the image field in an aplanatic system,” Proc. R. Soc. (London) A, 253 358 –379 (1959). https://doi.org/10.1098/rspa.1959.0200 Google Scholar

4. 

P. Török and P. R. T. Munro, “The use of Gauss–Laguerre vector beams in STED microscopy,” Opt. Express, 12 3605 –3617 (2004). https://doi.org/10.1364/OPEX.12.003605 OPEXFF 1094-4087 Google Scholar

5. 

K. S. Youngworth and T. G. Brown, “Focusing of high numerical aperture cylindrical-vector beams,” Opt. Express, 7 77 –87 (2000). https://doi.org/10.1364/OE.7.000077 OPEXFF 1094-4087 Google Scholar

6. 

P. Török, P. Varga and G. R. Booker, “Electromagnetic diffraction of light focused through a planar interface between materials of mismatched refractive indices—structure of the electromagnetic field I,” J. Opt. Soc. Am. A, 12 (10), 2136 –2144 (1995). https://doi.org/10.1364/JOSAA.12.002136 JOAOD6 0740-3232 Google Scholar

7. 

A. S. van de Nes et al., “Calculation of the vectorial field distribution in a stratified focal region of a high numerical aperture imaging system,” Opt. Express, 12 (7), 1281 –1293 (2004). https://doi.org/10.1364/OPEX.12.001281 OPEXFF 1094-4087 Google Scholar

8. 

R. Kant, “Vector diffraction problem of focusing a category 1 aberrated wavefront through a multilayered lossless medium,” J. Mod. Opt., 51 (3), 343 –366 (2004). https://doi.org/10.1080/09500340408235528 JMOPEW 0950-0340 Google Scholar

9. 

O. Haeberlé, “Focusing of light through a stratified medium: a practical approach for computing microscope point spread functions: Part II: confocal and multiphoton microscopy,” Opt. Commun., 235 (1–3), 1 –10 (2004). https://doi.org/10.1016/j.optcom.2004.02.068 OPCOB8 0030-4018 Google Scholar

10. 

P. R. T. Munro, A. Curatolo and D. D. Sampson, “Full wave model of image formation in optical coherence tomography applicable to general samples,” Opt. Express, 23 2541 –2556 (2015). https://doi.org/10.1364/OE.23.002541 OPEXFF 1094-4087 Google Scholar

11. 

E. Z. Y. Zhang and P. C. Beard, “Ultrahigh-sensitivity wideband Fabry–Perot ultrasound sensors as an alternative to piezoelectric PVDF transducers for biomedical photoacoustic detection,” Proc SPIE, 5320 222 –229 (2004). https://doi.org/10.1117/12.531160 Google Scholar

12. 

P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt., 36 (11), 2305 –2312 (1997). https://doi.org/10.1364/AO.36.002305 APOPAI 0003-6935 Google Scholar

13. 

I. Gregor and J. Enderlein, “Focusing astigmatic Gaussian beams through optical systems with a high numerical aperture,” Opt. Lett., 30 2527 –2529 (2005). https://doi.org/10.1364/OL.30.002527 OPLEDP 0146-9592 Google Scholar

14. 

D. G. Flagello, T. Milster and A. E. Rosenbluth, “Theory of high-NA imaging in homogeneous thin films,” J. Opt. Soc. Am. A, 13 53 –64 (1996). https://doi.org/10.1364/JOSAA.13.000053 JOAOD6 0740-3232 Google Scholar

15. 

D. G. Flagello and T. Milster, “3D modelling of high numerical aperture imaging in thin films,” Proc. SPIE, 1625 246 –261 (1992). https://doi.org/10.1117/12.58952 PSISDG 0277-786X Google Scholar

16. 

P. R. T. Munro, “Multi-layer vectorial focusing,” (2018) http://prtmunro.net September ). 2018). Google Scholar

17. 

P. Török, “An imaging theory for advanced, high numerical aperture optical microscopes,” (2003). Google Scholar

18. 

P. Higdon, P. Török and T. Wilson, “Imaging properties of high aperture multiphoton fluorescence scanning microscopes,” J. Microsc., 193 (2), 127 –141 (1999). https://doi.org/10.1046/j.1365-2818.1999.00448.x JMICAR 0022-2720 Google Scholar

19. 

R. Azzam and N. Bashara, Ellipsometry and Polarized Light, 3rd ed.Elsevier, North Holland, Amsterdam (1992). Google Scholar

20. 

C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles, Wiley Interscience, New York (1983). Google Scholar

21. 

J. Stratton, Electromagnetic Theory, McGraw-Hill, New York (1941). Google Scholar

22. 

P. R. Philip, J. Davis and W. Rheinbolt, Methods of Numerical Integration, Academic Press Limited(1984). Google Scholar

23. 

O. Bruno and J. Chaubell, “One-dimensional inverse scattering problem for optical coherence tomography,” Inverse Prob., 21 499 –524 (2005). https://doi.org/10.1088/0266-5611/21/2/006 INPEEY 0266-5611 Google Scholar

24. 

P. R. T. Munro, “Three-dimensional full wave model of image formation in optical coherence tomography,” Opt. Express, 24 27016 –27031 (2016). https://doi.org/10.1364/OE.24.027016 OPEXFF 1094-4087 Google Scholar

25. 

P. R. T. Munro, “Exploiting data redundancy in computational optical imaging,” Opt. Express, 23 30603 –30617 (2015). https://doi.org/10.1364/OE.23.030603 OPEXFF 1094-4087 Google Scholar

26. 

D. J. Segelstein, “The complex refractive index of water,” University of Missouri–Kansas City, (1981). Google Scholar

27. 

D. P. Biss and T. G. Brown, “Cylindrical vector beam focusing through a dielectric interface,” Opt. Express, 9 (10), 490 –497 (2001). https://doi.org/10.1364/OE.9.000490 OPEXFF 1094-4087 Google Scholar

28. 

K. Hu, Z. Chen and J. Pu, “Generation of super-length optical needle by focusing hybridly polarized vector beams through a dielectric interface,” Opt. Lett., 37 3303 –3305 (2012). https://doi.org/10.1364/OL.37.003303 OPLEDP 0146-9592 Google Scholar

29. 

D. Khoptyar et al., “Tight focusing of laser beams in a λ/2-microcavity,” Opt. Express, 16 9907 –9917 (2008). https://doi.org/10.1364/OE.16.009907 OPEXFF 1094-4087 Google Scholar

30. 

A. van de Nes et al., “Cylindrical vector beam focusing through a dielectric interface: comment,” Opt. Express, 12 967 –969 (2004). https://doi.org/10.1364/OPEX.12.000967 OPEXFF 1094-4087 Google Scholar

Biography

Peter R. T. Munro received his BSc (computer science) and BEng (electrical engineering) degrees from the University of Western Australia in 2000. He completed his PhD in the Department of Physics at Imperial College in 2006. Currently, he holds a Royal Society University research fellowship focused on performing quantitative imaging using a combination of x-ray, optical, and photoacoustic techniques.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE)
Peter R. T. Munro "Tool for simulating the focusing of arbitrary vector beams in free-space and stratified media," Journal of Biomedical Optics 23(9), 090801 (24 September 2018). https://doi.org/10.1117/1.JBO.23.9.090801
Received: 17 April 2018; Accepted: 15 August 2018; Published: 24 September 2018
Lens.org Logo
CITATIONS
Cited by 5 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Refractive index

Interfaces

Optical coherence tomography

Optical simulations

Glasses

Water

Wave propagation

RELATED CONTENT

Tuning Bessel beam propagation properties with liquid media
Proceedings of SPIE (September 09 2019)
Surface enhanced biodetection on a CMOS biosensor chip
Proceedings of SPIE (February 02 2012)
Modelization of a semi leaky waveguide application to a...
Proceedings of SPIE (December 01 1991)
Confocal TIRF microscopy of single molecules
Proceedings of SPIE (June 19 2003)

Back to Top