Differential phase contrast (DPC) microscopy, with its simplicity, non-invasiveness, and label-free imaging method, plays a pivotal role in optical microscopy. First introduced in scanning microscopy, DPC technique leverages a split detector to measure the differential intensity influenced by a sample’s phase gradient [1, 2]. Further advancements enabled widefield DPC imaging on a conventional brightfield microscope platform with split illumination [3], fostering real-time imaging vital for live-cell research. The understanding of the relationship between a sample’s phase information and the resulting differential intensity has made DPC imaging invaluable for quantitative phase imaging (QPI). As a result, it is possible to visualize in real-time, intrinsic optical phase changes of biological samples tied to dynamic processes like cell growth, movement, and behavior in their natural state. Recent advancements in widefield DPC microscopy techniques using coded illumination via light-emitting diode (LED) arrays [4–6] and liquid crystal display (LCD) [7–9], as well as coded detection [10, 11], have successfully showcased superior QPI capabilities for live cells and organisms [8–19].

However, despite these advancements in widefield DPC microscopy, the prevailing DPC imaging theory for phase retrievals in QPI still relies on paraxial diffraction theory [1, 2]. While effective for low numerical aperture (NA) DPC microscopy systems, this theory becomes inaccurate for high NA systems, which are critical for high-resolution research at the sub-organelle level.

To address this gap, we introduce a nonparaxial DPC imaging theory for high NA DPC microscopy, based on scalar diffraction theory. We derive a nonparaxial DPC imaging equation with transmission cross-coefficient (TCC), considering factors such as the apodization of the high NA objective lens, the nonparaxial propagation of light, and the angular distribution of illumination or detection strength. We also present nonparaxial TCCs for two practical sample object simplifications and compare these with traditional paraxial TCCs. Our work establishes a theoretical framework with potential to enhance phase retrieval accuracy in high NA DPC microscopy.

2.1. Schematic of DPC Microscopy

Figure 1 presents a scanning DPC microscopy system, a reciprocal equivalent to the widefield DPC microscopy system as per Helmholtz’s reciprocity theory [3]. The existing paraxial DPC imaging theory [1, 2] has its roots in such a scanning DPC microscopy system and has been aptly applied to widefield DPC imaging without compromising on generality. In sample-scanning DPC microscopy, an on-axis point source illuminates a sample that is in motion. The light field transmitted from the sample alters its propagation direction based on the local phase gradient of the sample. The split detector, in turn, detects this laterally shifted field differentially across all scanning positions sequentially. Thus, the phase gradient of the sample correlates with the differential intensity in DPC microscopy. Widefield DPC microscopy operates in reverse; each pixel of a camera (serving as a point detector) captures intensity differences at its conjugate sample locations between two illuminations from the split source. This enables a direct two-dimensional (2D) DPC imaging without scanning. We note that our schematic doesn’t account for a collector lens placed between the sample and the split detector, as many recent DPC microscopy setups employing LED array sources operate without a condenser lens [12–19].

Figure 1. A schematic of a scanning differential phase contrast (DPC) imaging system, with the light’s pathway from a point source to a split detector, labeled as “L” and “R”. This configuration is the reciprocal equivalent of a widefield DPC imaging system, where light proceeds from a split source to a point detector array. t(x): sample transmittance, P(ξ): objective pupil, P_d(x_d): split source/detector pupil. The vectors x = (x, y), ξ = (ξ, η), and x_d = (x_d, y_d) denote 2D coordinates.

Our formulation incorporates three nonparaxial factors compared with the paraxial approach. Firstly, we contemplate a high numerical aperture (NA) aplanatic objective. Its in-focus point spread function (PSF) in the sample plane, denoted as h(x), can be approximated using the scalar Debye integral [20]:

where k is the wavenumber, P(θ, ϕ) is the pupil function of the objective lens in spherical coordinate (θ, ϕ) with an apodization factor of $\sqrt{\cos \theta }$, and N is a unit normal vector to the objective pupil surface (Σ) pointing towards the geometrical focus. The Debye approximation is valid for the Fresnel number much greater than one [21], a condition met by high NA objectives. Here, the light field incident on the objective pupil from an on-axis isotropic point source is assumed to be uniform in amplitude across the aperture, attributed to the very low NA value (<0.03) of a typical tube lens in compound microscopes. Therefore, Eq. (1) can depict the quasi-monochromatic PSF for the point illumination on the sample. Upon transforming the integral geometry to the Cartesian pupil coordinate, Eq. (1) becomes

where f_o signifies the focal length of the objective lens. The in-focus PSF can thus be calculated by a 2D Fourier transform of an effective pupil function:

where P(ξ) denotes a conventional paraxial pupil function that is a circle function, circ(ξ). The nonparaxial apodization factor results in an effective pupil function with more intense weighing at the edge of the pupil in the Cartesian coordinate.

Secondly, our formulation considers nonparaxial light propagation from the sample to the detector (Fig. 2). When their distance (d) is much greater than the light’s wavelength, the complex amplitude at the detector plane could be expressed by the first Rayleigh-Sommerfeld diffraction integral [20] as

Figure 2. The diffraction geometry for light propagation from the sample plane to the detector plane.

where E(x) is the scalar field right after the sample, defined over a planar diffraction geometry Σ_S with the normal vector N = (0, 0, 1). R = (x_d − x, y_d − y, d) is a distance from a sample position to a detector position with $R=\sqrt{d^2+{\left|x_d-x\right|}^2}$. Assuming the sample’s dimension is much smaller than either d or the split detector radius, b, (or x is close to the origin (O) in Fig. 2), a Debye-like approximation [22] can be applied such that R ≈ M(R_d − x∙xd/R_d), where M = (x_d/R_d, y_d/R_d, d/R_d) is a unit vector from the origin to a detector point with $R_{\text{d}}=\sqrt{d^2+{\left|x_d\right|}^2}$. With a further approximation, 1/R ≈ 1/R_d, Eq. (4) simplifies to

It is noteworthy that Eq. (5) reflects the inverse square law (1/R_d) and the inclination factor (d/R_d) in the Huygens-Fresnel principle, both of which are neglected in the paraxial approximation.

Thirdly, we consider an angle-dependent light sensitivity of the split detector, which could be modeled as

where θ stands for a polar angle relative to the detector’s perpendicular axis, and g shapes the angular distribution of detection sensitivity (Fig. 3). In the context of widefield DPC microscopy, g corresponds to the Lambertian order [23] of a light source, controlling the angular distribution of emission intensity; it is isotropic for g = 0 and Lambertian for g = 1 (Fig. 3). For a planar LED array, the Lambertian order is mainly influenced by the lens shape encapsulating the LED chip, and it is crucial to consider this factor because each LED has a distinct polar angle on the sample.

Figure 3. Angular distribution of detector sensitivity or source intensity, depicted in polar coordinates, based on the g factor.

2.2. Nonparaxial DPC Imaging Equation

Given a thin sample at a scanning position x_s = (x_s, y_s) with amplitude transmittance represented as t(x), the electric field right after the sample could be approximated as h(x)t(x_s − x). Following its nonparaxial propagation to the detector plane, as guided by Eq. (5), this transmitted field may evolve to

Assuming a perfectly incoherent detector, mirroring the incoherence of the LED array source in a reciprocal manner, the total detected intensity at the sample position x_s can be calculated as

where D(x_d) = (d/R_d)^g |P_d(x_d)|² is an intensity detection sensitivity, angularly dependent on (d/R_d)^g from Eq. (6). Here, |P_d(x_d)|² corresponds to the shape of the detector, often a circle, and is weighted by a sign function, sgn(x_d), for differential detection in typical DPC imaging. We note that an LED array source in widefield DPC microscopy could be modeled by multiplying D(x_d) with a 2D Dirac comb function, ∑ _j δ(x_d − x_j), where x_j represents the location of the j^th LED.

To derive a TCC [24], also termed as a partially coherent transfer function in spatial frequency domain, we can use the Fourier domain representations of both the sample transmittance and the PSF as follows:

where T(f) denotes the Fourier transform of t(x) in the 2D spatial frequency domain f = (f_x, f_y), and λ is the wavelength of light. After substituting these two equations to Eq. (8) and performing further mathematical manipulations, a DPC imaging equation can be expressed in the frequency domain as

where the asterisk (*) stands for the complex conjugate and C(f₁, f₂) is the four-dimensional (4D) TCC derived as

Thus, DPC imaging is governed by a 4D bilinear partially coherent process. In contrast to the paraxial TCC dictated by the objective and detector pupils, P(ξ) and P_d(x_d), [1, 2], Eq. (12) introduces three additional elements. First is the high NA objective’s apodization factor, (1 − |ξ|²/f_o²)^−1/4, incorporated in the effective objective aperture, P_e(ξ) as seen in Eq. (3). The second addition is the nonparaxial propagation factor, d ²/R_d⁴, tracing the path from the sample to the detector, elucidating the inverse square law of intensity and the inclination factor in light diffraction. The last is the angular dependence factor of detector responses, (d/R_d)^g, incorporated in D(x_d). We note that, in widefield DPC microscopy, the scanning position x_s in Eq. (11) can be directly replaced by the sample coordinate x and D(x_d) can be perceived as the illumination’s intensity distribution across the condenser aperture.

2.3. TCC in Normalized Coordinates

It is often convenient to describe TCC defined in a normalized pupil coordinate as ξ′ = ξ/a, where a denotes the objective lens’s pupil radius, and in a normalized spatial frequency as m = f/(NA/λ) with NA = a/f_o for an aplanatic objective lens. This could be done by relating the detector coordinate with the objective pupil coordinate via

The Jacobian J(ξ′) for this coordinate transformation, i.e., dx_d = J(ξ′)dξ′, is derived as a² R_d⁴/( f_o² d²). Plugging Eq. (13) to $R_{\text{d}}=\sqrt{d^2+{\left|x_d\right|}^2}$ yields $R_{\text{d}}^{-1}=d^{-1}\sqrt{1-{\text{NA}}^2{\left|{\xi }^{\prime }\right|}^2}$. Using these relations, Eq. (12) can be reformed into a 2D integral over the normalized objective pupil coordinate as

If the detector has a radius of b and is split by a sign function (two semicircles with 1 and −1), P_d[x_d(ξ′)] can be simplified to P_d(bξ′/σ) [22]. Here, σ is the partial coherent factor, defined as the ratio of a detection NA to an objective NA: $b/\sqrt{d^2+b^2}/\text{NA}$. For the rotationally symmetric objective pupil, namely, P(ξ′) = P(−ξ′), TCC has a form of

where the ‘n’ subscript in each pupil function denotes the normalized pupil to the unit circle.

The distinction between paraxial and nonparaxial TCCs become evident upon close examination of Eq. (15). In the frequency domain, the paraxial TCC represents the geometric overlap of three pupils with uniform amplitudes [24], whereas the nonparaxial TCC consists of three pupils with radially varying amplitudes. These comprise a detector pupil whose sensitivity radially decreases by (1 − NA² |ξ′|²)^g/2 and the two displaced objective pupils apodized with (1 − NA² |ξ′|²)^−1/4. As expected, the TCC derived from these radially nonuniform pupils exhibits a dependency on NA and converges to the paraxial TCC when NA approaches zero. It is interesting to note that when using an isotropic detector (g = 0), or reciprocally, an isotropic LED source in wide-field DPC imaging, the detector sensitivity remains radially uniform irrespective of NA. If one considers a point array detector as represented by a Dirac comb, the TCC in Eq. (15) could be represented as

where ξ_j′ = x_j /NA∙(d² + |x_j|²)^−1/2 denotes the j^th location of the point detector.

In DPC imaging, the differential detection leads to a TCC value of zero at the zero spatial frequency, i.e., C(0;0) = 0, which precludes its use as a TCC normalization factor. Instead, a brightfield TCC, C_BF(f₁;f₂), when |P_d,n(ξ′/σ)|² = circ(ξ′/σ) in Eq. (15), could be employed as follows:

with the normalization factor at zero spatial frequency derived as

The same factor can be obtained directly from Eq. (12), allowing Eq. (18) to be used for normalizing both C(m₁;m₂) and C(f₁;f₂).

2.4. Two Possible Simplifications

The 4D bilinear process of 2D DPC imaging, represented by Eq. (11), is quite complex for practical purposes. As a workaround, it can be simplified to 2D imaging processes for two idealized objects [1, 2]. Here, we first consider a weak object expressed as t(x) = 1 + t_w(x) where |t_w(x)| ≪ 1. Replacing its Fourier transform, T(f) = δ(f) + T_w(f), into Eq. (11), using C(f₁;f₂) = C^*(f₂;f₁) due to I(x) = I^*(x) (as intensity is real), and neglecting a cross-product term [22], the DPC imaging intensity simplifies to

where C(0;0) is zero in DPC imaging and Re[ ] takes the real part. This adjustment leads to a linear 2D image formation, where the image intensity is directly tied to an inverse Fourier transform of the product of the weak object transfer function (WOTF) C(f;0) and the weak object frequency spectrum T_w(f).

The DPC imaging intensity is often normalized by the brightfield image intensity, such that

where I_L and I_R are the measured intensities on the left and right semicircle detectors, respectively, as shown in Fig. 1. With this normalization, the DPC imaging intensity of the weak object can be approximated [22] as

where C_N,LR(f;0) = C_LR(f;0)/C_BF(0;0) is the normalized 2D WOTF in the ‘LR’ DPC configuration. Here, the ‘LR’ subscript added to Eq. (17) indicates the inclusion of sgn(ξ′/σ) in the detector sensitivity. Hence, the WOTF model allows the object information to be inversely retrieved from the measured DPC image. In practice, DPC images are also gathered from the top/bottom detectors or multiple azimuthal orientations for improved spatial frequency coverage. These multi-orientation WOTFs are integrated into the retrieval algorithm [4].

Another idealized object of interest is a slowly varying phase object, whose spatial phase profile changes more slowly than the imaging resolution. In this case, the object transmittance is given by t(x) = exp(i∆ϕ(x)∙x), where ∆ϕ(x) is a local phase gradient. As the Fourier transform of t(x) is T(f) = δ(f − ∆ϕ(x)/(2π)), the phase gradient at x is directly mapped to a spatial frequency of ∆ϕ(x)/(2π). This allows simplification of Eq. (11) to

where C(f;f) is called a 2D phase gradient transfer function (PGTF). For the DPC intensity definition in Eq. (20), it can be shown that

where C_N,LR(f;f) = C_LR(f;f)/C_BF(f;f) serves the normalized 2D PGTF in the ‘LR’ DPC configuration. Consequently, the DPC intensity measured is directly associated with the magnitude of the PGTF at the spatial frequency that corresponds to the object’s local phase gradient.

We conducted numerical simulations with parameters λ = 0.5 μm and σ = 1 to contrast nonparaxial TCC with paraxial TCC in DPC imaging using the ‘LR’ split geometry. Our first step was calculating the normalized 2D WOTF, C_N,LR(m;0), across varying ranges of NA and g values, as shown in Fig. 4. Each WOTF was rapidly computed in 20–31 ms using MATLAB (on an Intel Core i7) via the convolution theorem [20] of the Fourier transform. Owing to the ‘LR’ split geometry, WOTFs exhibited antisymmetry about the vertical axis (m_x = 0), with the peak WOTF magnitude occurring near m = (±0.9, 0), where the displaced objective pupil has maximum overlap with one of the semicircle sources/detectors. In this scenario, less than half of the shifted objective pupil area resided within the semicircle, yielding a maximum WOTF value closer to 0.4. As NA and g increase, enhancing the nonuniformity of the objective and source/detector pupils, the 2D nonparaxial WOTF shape showed the greatest deviation from the paraxial WOTF, particularly at 0.9 NA with g = 2 (Fig. 4). For a more detailed comparison, we examined the WOTF horizontal cross-sections at m_y = 0 (Fig. 5). While the WOTFs at 0.3 NA closely resembled the paraxial WOTF with negligible deviations (<1.7%) across g = 0–2, those at 0.9 NA displayed substantial deviations of −8.0%, −16.2%, and −25.7% for g = 0, 1, and 2, respectively. Compared with the paraxial WOTF, nonparaxial WOTFs showed reduced intensity values (or lesser pupil overlap) over m_x = 0–0.9, and exhibited peak WOTF values that could fall either below or above the paraxial peak, depending on the NA and g values.

Figure 4. Normalized 2D WOTF for ‘LR’ DPC imaging (λ = 0.5 μm) across different NA and g values. The paraxial WOTF is included for comparison. WOTF, weak object transfer function; DPC, differential phase contrast; NA, numerical aperture.

Figure 5. Cross-sections of 2D weak object transfer function (WOTF) at m_y = 0. Deviations indicate the intensity difference between nonparaxial and paraxial WOTFs, normalized by the paraxial WOTF peak of 0.4097.

Next, we computed the normalized 2D PGTF, C_N,LR(m;m), over a range of NA and g values, as shown in Fig. 6. Similar to WOTFs, PGTFs were also antisymmetric about m_x = 0 (vertical axis) and were bound between −1 and 1 upon normalization with C_BF(m;m). This clamping of PGTF arises when two objective pupils, equally displaced by m, partially overlap with a single semicircle of the source/detector. For example, this clamping initiates when |m_x| ≥ 1 for m_y = 0, and any local phase gradient in a slowly varying phase sample with an absolute value of normalized frequency exceeding one becomes undifferentiable in Eq. (23). Additionally, the normalized PGTF may only be defined within the circle where |m| ≤ 2, as this is the region where C_BF(m;m) is non-zero. Compared with the 2D paraxial PGTF, nonparaxial PGTFs displayed significant differences, especially at higher NAs like 0.9 across all g values (Fig. 6). In the PGTF cross-sections at m_y = 0 (Fig. 7), PGTF values for all NAs were lower than their paraxial PGTF counterparts, with deviations at m_x = ~0.55 for 0.9 NA being as large as −11.6%, −20.1%, −27.2% when g = 0, 1, and 2, respectively. Much like the WOTFs, PGTFs at 0.3 NA closely resembled the paraxial PGTFs, with negligible deviations less than 1.8%.

Figure 6. Normalized 2D PGTF for ‘LR’ DPC imaging (λ = 0.5 μm) across different NA and g values. The paraxial PGTF is included for comparison. PGTF, phase gradient transfer function; DPC, differential phase contrast; NA, numerical aperture.

Figure 7. Cross-sections of 2D phase gradient transfer function (PGTF) at m_y = 0. Deviations indicate the intensity difference between nonparaxial and paraxial PGTFs, normalized by the paraxial PGTF maximum of 1.

Interestingly, we observed that nonparaxial PGTFs did not always show a monotonic increase along the m_x direction. For instance, a non-monotonic behavior was noticed near m = (±0.11, 0) at 0.9 NA and g = 2 in Fig. 6. This behavior was more evident in the PGTF cross-section shown in Fig. 7, where PGTF values near m_x = 0.11 were below zero (which was the value at m_x = 0). This reversal in PGTF originated from the large nonuniformity of pupil functions, evident in Fig. 8. While the geometrical overlap of the three pupils for m = (0.11, 0) was wider at the left semicircle (positive) detector, the dominant influence came from the intensified edges of the two apodized objective pupils on the right semicircle (negative) detector, resulting in a negative differential intensity, and hence, a negative PGTF value. In fact, the overall weighting factor in the net pupil overlap at m = (0.89⁻, 0) in Fig. 8 was −0.82, an absolute magnitude 2.6× larger than 0.32 at m = (−1⁺, 0). We further identified that a similar inversion in intensity can occur in nonparaxial WOTF, such as near m = (±0.05, 0) when NA = 0.95 and g = 2, although the reversal was less stark, given that only one objective pupil undergoes shifting in WOTF.

Figure 8. Detector and objective pupils, along with their multiplication, illustrating the net pupil overlap for PGTF(m; m) when m = (0.11, 0), g = 2, NA = 0.9, and σ = 1. Pupils are represented on normalized objective pupil coordinates, equivalent to normalized spatial frequencies. Cross-sectional profiles at m_y = 0 are additionally shown. PGTF, phase gradient transfer function.

Finally, we simulated high NA DPC imaging for a one-dimensional (1D) sinusoidal phase object (phase amplitude: 2 radians) as shown in Fig. 9. Using the 1D PGTF forward model (NA = 0.9 and g = 0), the high NA DPC intensity profile exhibited reduced variation compared to its paraxial counterpart, due to the lower intensity values of the nonparaxial PGTF (as discussed in Fig. 7). We then back-calculated the phase gradient from the nonparaxial DPC intensity profile using both 0.9 NA PGTF and paraxial PGTF models, respectively. The nonparaxial PGTF model yielded an accurate phase profile (in red), while the paraxial PGTF model resulted in significant deviations, with phase errors surpassing 60%. This highlights the potential risk of using paraxial TCC models in high NA QPI.

Figure 9. Quantitative phase imaging (QPI) simulation of a slowly varying 1D sinusoidal phase object.

We have, for the first time, derived a nonparaxial DPC imaging equation in the scanning DPC imaging geometry, considering the apodization of the objective lens, the nonparaxial propagation of light, and the angular variation in detection sensitivity. The formulation developed herein is also applicable to widefield DPC microscopy systems implemented with a split source with angularly varying emission, such as an LED array. In comparison to the paraxial TCC, our derived nonparaxial TCC introduces two additional parameters, NA and g. These parameters induce radial nonuniformity in the pupil functions across the aperture, differing from the uniform pupils traditionally seen in paraxial DPC imaging theory. Our numerical investigation of simplified TCCs, namely WOTF and PGTF, suggests that paraxial TCC can maintain accuracy at lower NAs (below 0.3) but deviates significantly at higher NAs and g values, reaching more than 25% deviation at 0.9 NA and g = 2. At such high NA and g values, our model also predicts the possible occurrence of TCC reversal, a phenomenon not anticipated by paraxial theory. Our high NA QPI simulation also reveals that misusing paraxial PGTFs could lead to significant phase reconstruction errors. Therefore, the nonparaxial TCCs derived in this work can potentially enhance phase retrieval accuracy in QPI based on high NA DPC imaging, which will propel cell biology research.

Further explorations could involve extensive numerical studies on different source shapes (e.g., annular split sources and LED array sources), source distribution variations (g), ranges of partial coherent sigma (σ), and spectral bandwidths of sources. Such parameter studies, enabled by the developed formulation, could deepen the understanding of nonparaxial DPC imaging and guide the development of optimal high NA DPC imaging systems to maximize imaging performance. Moreover, further theoretical development considering optical defocus, spherical aberrations, and light-specimen interactions that include the polarization states of light, could provide more profound insights into high NA DPC imaging and further enhance TCC accuracy.

This work was in part supported by the Research Institute for Convergence Science, Seoul National University.

Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2022R1A6A1A03063039).

The author declares no conflicts of interest.

Data underlying the results presented in this paper may be obtained from the authors upon reasonable request.