Computation Theory of Large-Scale Partially Coherent Imaging by the Modified Modal Expansion Method

Abstract

The numerical calculation of partially coherent imaging involves a fourfold integral, which is numerically complex and impracticable to be calculated directly. The use of coherent-mode decomposition (CMD) can make this problem more manageable but finding the coherent-modes for complicated partial coherent fields (without already-known coherent-mode expansion) are rather computationally intensive. In this letter, a modified modal expansion method is proposed, which significantly reduces the requirement of computational resources. The propagation of partial coherence in imaging systems with extremely large sampling number could be handled by an ordinary computer. A comparison between the new method and the traditional method in terms of memory resource requirements and computational time consumption is also detailed in this article. We will also show that this method could deal with the anisoplanatic imaging cases while maintaining the same computational efficiency.

1. Introduction

Partially coherent imaging was first introduced by Hopkins [1] with the introduction of a four-dimensional transmission factor and has been investigated intensively over the past few decades [2,3,4,5,6]. The simulation of partially coherent imaging is concerned with evaluating the following integrals, which is known as the cross-spectral density (CSD) function propagation law [7] (propagation through the imaging system in partially coherent imaging cases):

$$\begin{array}{l}W\left({\mathit{r}}_1,{\omega }_1,{\mathit{r}}_2,{\omega }_2\right)\approx 〈\left[{\displaystyle \underset{\infty }{\iint }U\left({\mathit{\rho }}_1^{\prime },{\omega }_1\right)G\left({\mathit{r}}_1,{\omega }_1;{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}\right]·{\left[{\displaystyle \underset{\infty }{\iint }U\left({\mathit{\rho }}_2^{\prime },{\omega }_2\right)G\left({\mathit{r}}_2,{\omega }_2;{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }}\right]}^{*}〉\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \iint {\displaystyle {\iint }_{\infty }〈U\left({\mathit{\rho }}_1^{\prime },{\omega }_1\right)U*\left({\mathit{\rho }}_2^{\prime },{\omega }_2\right)〉G\left({\mathit{r}}_1,{\omega }_1;{\mathit{\rho }}_1^{\prime }\right)G*\left({\mathit{r}}_2,{\omega }_2;{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \iint {\displaystyle {\iint }_{\infty }W\left({\mathit{\rho }}_1^{\prime },{\omega }_1,{\mathit{\rho }}_2^{\prime },{\omega }_2\right)G\left({\mathit{r}}_1,{\omega }_1;{\mathit{\rho }}_1^{\prime }\right)G*\left({\mathit{r}}_2,{\omega }_2;{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }}}\end{array}$$

(1)

where $\mathit{\rho }$ and $\mathit{r}$ are the position vectors of the source plane and observation plane, respectively. The CSD is a function of variables in either four or six dimensions, while the field is a function of variables in either two or three dimensions. The numerical complexity of the propagation calculation is exponential in the number of dimensions, so direct calculation of the propagation of the CSD may be impracticable.

The use of coherent-mode [8,9] or pseudo-modes [10] expansions of CSD can potentially make this problem more manageable. In the view of both methods, the partially coherent wide-sense stationary (WSS) field can be generated by the incoherent sum of many statistically independent modes. The probability of appearance of each mode, or the amplitude of each mode, is concerned with the incoherent sum by weighting the modes, which corresponds to two different physical interpretations. Although the basic principle of these two methods is somewhat similar, the realization of the two methods is completely different. The two methods both have advantages and drawbacks. For the method of coherent-mode expansion, we must know the modes of the field. However, the coherent-mode expansion of complicated fields may be difficult to find analytically. Up to now, only a small number of coherent-mode expansions have been found, which were summarized well by references [11,12]. The simulation of partially coherent imaging with modal expansion method was reported in reference [13]. The pseudo-modes for a given CSD function are much easier to find than its coherent-modes expansion and thus, this method can be used to simulate nearly any random light source. The application of pseudo-modes method in partially coherent imaging has been evaluated in reference [2,5]. However, the modes involved in pseudo-modes are not necessarily orthogonal nor unique, which is the reason why they are called pseudo-modes. Moreover, the pseudo-modes are not discrete. Both of these mean that many more pseudo-modes are required to accurately represent a random source than coherent modes. Therefore, for partially coherent fields whose modes are already known, the coherent-mode method is more desired in numerical simulations for the sake of efficiency and accuracy. Moreover, simulations of partially coherent imaging could also be accelerated in some particular methods, such as the transmission cross coefficient (TCC) method widely used in microlithography research [3,4]; however, these methods are only restricted in the applications of some certain specific types of partially coherent light fields.

For complicated fields whose modes are difficult to be found analytically, coherent-mode representations may also be found numerically by diagonalizing the Hermitian symmetric matrix [14] or using an LDL factorization [15]. The main challenge for looking for coherent modes numerically is the computational cost, especially regarding memory [16]. Let us consider 2D partial coherent imaging along spatial dimensions x and y; a sampling of 500 points are carried out along each axis on object plane. In this case, in order to diagonalize the CSD matrix W, the computer must be able to store a ${500}^2 \times {500}^2$ complex matrix, which means 500⁴ = 62.5$\times$10⁹ points of data or 250 gigabytes of fast memory at single precision. As a comparison, only 500² points of data or 1 Mb are needed for the cases of complete coherence or incoherence. Thus, numerically finding the coherent-modes expansion are practical only for small-scale CSD matrices, which correspond to small sampling numbers in each spatial dimension. For partial coherent fields with a large sampling number, it is not practical to find their coherent-modes numerically. The “large-scale” defined in the title of this paper refers to such cases where the storage and manipulation requirements of the CSD matrix could not be fulfilled by the computer. Surely, the specific criteria for the judgment of “large-scale” are related to the computer hardware configuration conditions.

In this letter, we propose a new model of partially coherent imaging. A modified coherent-mode decomposition method is proposed, which could be used for finding coherent-modes numerically with a significant reduction in the requirement of computational resources. And thus, this method could be easily handled numerically, while maintaining the computational efficiency of the classical coherent mode decomposition method over that of other methods such as the pseudo-modes method.

2. Theory

2.1. A Review of the Theory of Partially Coherent Imaging with Coherent-Mode Representation

First, let us assume an imaging system with partially coherent illumination on the object. We specify the object by a transmission function $t\left(\mathit{\rho }\right)$. The CSD (cross-spectral density) of the stationary field incident on the object is $W\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)=〈u({\mathit{\rho }}_1)u^{*}({\mathit{\rho }}_2)〉$, where $\mathit{\rho }=(x,y)$, $u(\mathit{\rho })$ is a single field realization and the brackets denote ensemble averaging. The CSD of the field emerging from the object is

$$W_0\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)=t\left({\mathit{\rho }}_1^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_2^{\prime }\right)W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)$$

(2)

Let $K(\mathit{\rho },\mathit{\rho }\prime )$ be the transfer function characterizing light propagation from a point $\mathit{\rho }$ to an image point $\mathit{\rho }\prime$. The image-plane CSD $W^i\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)$ can be expressed as [17]

$$W^i\left({\mathit{\rho }}_1,{\mathit{\rho }}_2\right)={\displaystyle \int {\displaystyle {\int }_{-\infty }^{\infty }{W}_0\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)}}K\left({\mathit{\rho }}_1,{\mathit{\rho }}_1^{\prime }\right)K^{\ast }\left({\mathit{\rho }}_2,{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }$$

(3)

And the intensity on the image plane could be obtained by

$$I\left(\mathit{\rho }\right)=W^i\left(\mathit{\rho },\mathit{\rho }\right)$$

(4)

Let us consider the situation where $t\left(\mathit{\rho }\right)$ is only non-zero in a finite area A, which is the usual case in real imaging systems. The integral in Equation (3) outside area A is zero. Thus, we are allowed to restrict the integral in Equation (3) within area A, which is equivalent to the original equation. Substituting Equations (2) and (3) into the expression in Equation (4), the intensity on the image plane could be expressed as [18,19]

$$I\left(\mathit{\rho }\right)={\displaystyle \underset{A}{\iint }W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)}K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right)K^{\ast }\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }$$

(5)

here, the integral has been restricted in area A.

It is well known that the coherent modes of a genuine CSD can be found via a Mercer (eigen-value) expansion of the Hermitian and non-negative-defined CSD [14,20,21],

$$W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)={\displaystyle \sum _n{\lambda }_n{\psi }_n\left({\mathit{\rho }}_1^{\prime }\right){\psi }_n^{*}\left({\mathit{\rho }}_2^{\prime }\right)}$$

(6)

where ${\lambda }_n$ denote the Eigenvalues and ${\psi }_n$ denotes the corresponding Eigenfunctions of $W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)$.

With this series representation in Equation (6), the fourfold integral of Equation (5) can be reduced into a summation of twofold integrals:

$$I\left(\mathit{\rho }\right)={\displaystyle \sum _n{\lambda }_n{\phi }_n\left(\mathit{\rho }\right){\phi }_n^{*}\left(\mathit{\rho }\right)}$$

(7)

where ${\phi }_n$ are the image plane modes propagated from the object of the imaging system [22], formed by a two-fold integral in 2D plane coordinates:

$${\phi }_n\left(\mathit{\rho }\right)={\displaystyle {\int }_{A}t\left({\mathit{\rho }}_1^{\prime }\right){\psi }_n\left({\mathit{\rho }}_1^{\prime }\right)K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}$$

(8)

2.2. Computation Theory of Large-Scale Partial Coherent Imaging with Modified Modal Expansion

In many cases, the coherent-mode expansions of $W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_1^{\prime }\right)$ in Equation (6) need to be found numerically, if the coherent-modes are not known analytically. If a large-scale object is to be considered, which means a large number of sampling numbers, the decomposition would require huge storage of the computer, which is often not applicable.

The solution to this problem proposed in this letter is to split the illuminated area A into P subregions ${\mathit{A}}_p$. Since Equation (5) is a double integral over area A, it could be rewritten as the sum of the integrals over each subregion and the integrals across each pair of different subregions, as follows:

$$I\left(\mathit{\rho }\right)=I_1\left(\mathit{\rho }\right)+I_2\left(\mathit{\rho }\right)$$

(9)

where $I_1$ is the summation of the double integral over each subregion ${\mathit{A}}_p$, as follows:

$$I_1\left(\mathit{\rho }\right)={\displaystyle \sum _p{\displaystyle \underset{A_p}{\iint }W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)}}K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right)K^{\ast }\left(\mathit{\rho },{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }$$

(10)

and $I_2$ is the summation of the double integral across different subregions, as follows:

$$I_2\left(\mathit{\rho }\right)={\displaystyle \sum _p{\displaystyle \sum _{q\ne p}{\displaystyle \underset{A_p}{\int }{\displaystyle \underset{A_q}{\int }W\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)}}}}K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right)K^{\ast }\left(\mathit{\rho },{\mathit{\rho }}_2^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }$$

(11)

It should be noticed that, for the special case of completely incoherent imaging, Equation (11) does not need to be considered, since there is no coherence between different subregions.

The above summation could also be understood from the perspective of CSD matrix W. If the illuminated area A is split into P subregions, the CSD matrix W could be rewritten into $P^2$ submatrices, as follows:

$$\mathit{W}=\left[\begin{array}{ccccc}{\mathit{W}}_{11}& {\mathit{W}}_{12}& \dots & \dots & {\mathit{W}}_{1p}\\ {\mathit{W}}_{21}& {\mathit{W}}_{22}& & & ⋮\\ ⋮& & \ddots & & ⋮\\ ⋮& & & \ddots & ⋮\\ {\mathit{W}}_{p1}& \dots & \dots & \dots & {\mathit{W}}_{pp}\end{array}\right]=\left[\begin{array}{ccccc}〈{\left|{\tilde{u}}_1\right|}^2〉& 〈{\tilde{u}}_1{\tilde{u}}_2^{*}〉& \dots & \dots & 〈{\tilde{u}}_1{\tilde{u}}_p^{*}〉\\ 〈{\tilde{u}}_2{\tilde{u}}_1^{*}〉& 〈{\left|{\tilde{u}}_2\right|}^2〉& & & ⋮\\ ⋮& & \ddots & & ⋮\\ ⋮& & & \ddots & ⋮\\ 〈{\tilde{u}}_p{\tilde{u}}_1^{*}〉& \dots & \dots & \dots & 〈{\left|{\tilde{u}}_p\right|}^2〉\end{array}\right]$$

(12)

where $u_p$ represents the random optical field in each subregion, the brackets denote ensemble averaging, and the tildes in equation indicate filling in the number sequence. The diagonal submatrices, ${\mathit{W}}_{pp}(p=1\dots P)$, are the CSD of each subregion; whereas, the other submatrices ${\mathit{W}}_{pq}(p\ne q)$ present the mutual interference between different subregions. We should notice that ${\mathit{W}}_{pq}$ are zero matrices for fully incoherent light. Moreover, ${\mathit{W}}_{pq}$ and ${\mathit{W}}_{qp}$ are transposed conjugated.

If we are considering the case of discrete optical fields, the integral in Equation (5) should be carried out by a numerical integration over the entire matrix W. When dealing with large-scale W matrices, performing integration directly on the entire matrix can be time-consuming and computationally intensive. A common strategy is to divide the large matrix into smaller submatrices ${\mathit{W}}_{pq}$ and then perform integration separately on each submatrix. This is another interpretation of the process in Equations (9)−(11) from the discrete representation.

First, let us consider the integrations on each diagonal submatrix or the summation in Equation (10). It could be easily found that they have the same form as Equation (5) and thus, they could be accelerated numerically with the coherent-modes method, by following the steps of Equations (5)−(8). Compared with performing the integration on the original CSD matrix W, the new calculation method is significantly less computationally intensive and has less a reduced demand for memory space.

The next step is to perform the integration over each non-diagonal submatrices, which correspond to the continuous representation in Equation (11). Since the direct integration is still a daunting task, we need to find a more efficient method. Traditional methods for accelerating the integration on the CSD matrix, such as the coherent mode method, are no longer applicable because the non-diagonal matrices ${\mathit{W}}_{pq}=〈{\tilde{u}}_p{\tilde{u}}_q^{*}〉$ are not standard CSD matrices, as they are not Hermitian symmetric. The physical meaning of these non-diagonal matrices also differs from that of standard CSD matrices, describing the characteristics of mutual coherence between the light fields of different subregions.

In order to accelerate the integrations on these non-diagonal submatrices, we would also like to find a similar expansion form of ${\mathit{W}}_{pq}$. The expansion could be performed with singular-value decomposition (SVD). It should also be noted that the Mercer decomposition in coherent-mode decomposition is a special case of singular-value decomposition [14].

By utilizing the SVD method, we could diagonalize the non-Hermitian symmetric matrix ${\mathit{W}}_{pq}$ as follows:

$${\mathit{W}}_{pq}{=\mathit{U}}_{pq}{\Lambda }_{pq}{\mathit{V}}_{pq}^{†}$$

(13)

where ${\Lambda }_{pq}$ is the diagonal matrix consisting the singular values of matrix ${\mathit{W}}_{pq}$. This expression can be rewritten in a more explicit form, as follows:

$$W_{pq}\left(m,n\right)={\displaystyle \sum _{k=1}^N{\lambda }_{kpq}}{u}_{kpq}\left(m\right)v_{kpq}^{†}\left(n\right)$$

(14)

where $W_{pq}\left(m,n\right)$ denotes the elements of ${\mathit{W}}_{pq}$, ${\lambda }_{kpq}$ denotes the $k$th singular value of ${\mathit{W}}_{pq}$, and $u_{kpq}\left(m\right)$ and $v_{kpq}\left(n\right)$ denote the $m$th and $n$th element of the $k$th column of ${\mathit{U}}_{pq}$ and ${\mathit{V}}_{pq}$. Equation (14) can also be written in the following continuous form:

$$W_{pq}\left({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime }\right)={\displaystyle \sum _n{\lambda }_{npq}{\psi }_{npq}\left({\mathit{\rho }}_1^{\prime }\right){\zeta }_{npq}^{*}\left({\mathit{\rho }}_2^{\prime }\right)}$$

(15)

where ${\psi }_{npq}$ and ${\zeta }_{npq}$ are the coherent modes corresponding to subregion ${\mathit{A}}_p$ and ${\mathit{A}}_q$, decomposed from the mutual coherence between these two subregions.

Comparing Equation (15) with Equation (6), we could find that the decomposed modes ${\psi }_{npq}$ and ${\zeta }_{npq}$ in Equation (15) are not transposed conjugated to each other, which results from the differences in characters between ${\mathit{W}}_{pp}$ and ${\mathit{W}}_{pq}$. Therefore, the propagation of the decomposed modes through the imaging system should be calculated, respectively, as follows:

$$u_{npq}\left(\mathit{\rho }\right)={\displaystyle {\int }_{A_p}t\left({\mathit{\rho }}_1^{\prime }\right){\psi }_{npq}\left({\mathit{\rho }}_1^{\prime }\right)K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}$$

(16)

$$v_{npq}\left(\mathit{\rho }\right)={\displaystyle {\int }_{A_q}t\left({\mathit{\rho }}_1^{\prime }\right){\zeta }_{npq}\left({\mathit{\rho }}_1^{\prime }\right)K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}$$

(17)

here, $u_{npq}$ and $v_{npq}$ are the propagated modes at the imaging plane. Then, the image intensity contributed by the mutual coherence between different subregions could be expressed as follows:

$$\begin{array}{ll}{I}_2(\mathit{\rho })& =\sum _p\sum _{q\ne p}{\int }_{A_p}{\int }_{A_q}W({\mathit{\rho }}_1^{\prime },{\mathit{\rho }}_2^{\prime })t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)K(\mathit{\rho },{\mathit{\rho }}_1^{\prime })K^{\ast }(\mathit{\rho },{\mathit{\rho }}_2^{\prime }){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }\\ & =\sum _p\sum _{q\ne p}\sum _n{\lambda }_{nqq}{\int }_{A_p}{\int }_{A_q}{\psi }_{nqq}\left({\mathit{\rho }}_1^{\prime }\right){\zeta }_{nqq}^{\ast }\left({\mathit{\rho }}_2^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)K(\mathit{\rho },{\mathit{\rho }}_1^{\prime })K^{\ast }(\mathit{\rho },{\mathit{\rho }}_2^{\prime }){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }\\ & =\sum _p\sum _{q\ne p}\sum _n{\lambda }_{nqq}{\int }_{A_p}{\psi }_{rpq}\left({\mathit{\rho }}_1^{\prime }\right)t^{\ast }\left({\mathit{\rho }}_1^{\prime }\right)K(\mathit{\rho },{\mathit{\rho }}_1^{\prime }){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }{\int }_{A_q}{\zeta }_{rpq}^{\ast }\left({\mathit{\rho }}_2^{\prime }\right)t\left({\mathit{\rho }}_2^{\prime }\right)K^{\ast }(\mathit{\rho },{\mathit{\rho }}_2^{\prime }){\mathbf{d}}^2{\mathit{\rho }}_2^{\prime }\\ & =\sum _p\sum _{q\ne p}\sum _n{\lambda }_{npq}{u}_{npq}(\mathit{\rho })v_{rpq}^{\ast }(\mathit{\rho })\end{array}$$

(18)

The SVD-based method helps simplify the integral mathematically. The partially coherent light fields on each pair of subregions, ${\mathit{A}}_p$ and ${\mathit{A}}_q$, whose mutual coherence satisfies the mathematical description of the ${\mathit{W}}_{pq}$ matrix, are found mathematically in the form of two sets of weighted spatial coherent modes ${\psi }_{npq}$ and ${\zeta }_{npq}$. We should notice that the partially coherent light fields summed from these weighted coherent modes ${\psi }_{npq}$ and ${\zeta }_{npq}$ are not equivalent to the real partially coherent light field on each subregion. An extreme example is that when we are considering two subregions that are incoherent with each other, the light field constructed on each subregion would be zero light field, for the fulfillment of a zero ${\mathit{W}}_{pq}$ matrix. And during the simulation process, we would construct several different partially coherent light fields on the same subregion ${\mathit{A}}_p$, to help describe its mutual coherence characters with every other subregion ${\mathit{A}}_q$.

Finally, by following Equation (9), the complete image of the partially coherent system can be computed and the entire calculation is finished.

2.3. Summary of the Computation Procedure

As a summary, the procedure for executing a large-scale partial coherent imaging simulation is as follows:

• Split the original partially coherent illuminated area $\mathit{A}$ into P subregions ${\mathit{A}}_1$~${\mathit{A}}_P$;
• Divide the original CSD matrix W into several sub-matrices according to the division of the lighting area as Equation (12). Through this step, we can obtain P diagonal sub-matrices ${\mathit{W}}_{pp}$ representing the CSD of each subregion and $P·(P-1)$ sub-matrices ${\mathit{W}}_{pq}$ representing the degree of coherence between different subregions;
• Perform CMD decomposition on each diagonal sub-matrix ${\mathit{W}}_{pp}$ and calculate the image plane intensity of each CSD sub-matrix propagated through the imaging system by following a similar procedure as Equations (5)−(8);
• Perform SVD decomposition on each off-diagonal submatrix ${\mathit{W}}_{pq}$, decomposing each submatrix into two sets of weighted spatially coherent light field modes ${\psi }_{npq}$ and ${\zeta }_{npq}$, by following Equation (14). Through this mathematical decomposition method, we can construct two partially coherent light fields in each pair of subregions, whose cross-correlation properties satisfy the mathematical description of the ${\mathit{W}}_{pq}$ matrix;
• Propagate the weighted spatial coherent modes ${\psi }_{npq}$ and ${\zeta }_{npq}$ to the imaging plane, obtaining the propagated light fields $u_{npq}$ and $v_{npq}$ on the entire image plane, by following Equations (16) and (17);
• Calculate the image plane intensity contributed by the mutual coherence between each pair of subregions, by following Equation (18);
• Calculate the total intensity in the image plane, by following Equation (9).

3. Simulation Results

3.1. A Tutorial of Simulation with the Proposed Method (2 × 2 Splitting Case)

In this section, we present an example of imaging with a partially coherent system. Here, we assume Kohler illumination before the imaging system illustrated in Figure 1. It should be noticed that the assumption of Kohler illumination is given as an example and that the method proposed in this letter is not restricted to any specific type of partially coherent source.

The partially coherent imaging of a square object is considered, as a simplest tutorial. The object is split into 2 × 2 subregions (${\mathit{A}}_1$~${\mathit{A}}_4$), as shown in Figure 1.

First, by utilizing Equation (10), we calculate the image formed by each subregion after passing through the optical system. With the traditional coherent mode decomposition method, we can obtain four series of weighted coherent modes equivalent to the partially coherent light field of the four subregions. As an example, we present the first 14 eigenmodes of the subregion ${\mathit{A}}_1$ obtained through the coherent mode expansion in Figure 2. In practical calculations, to achieve more accurate results, it is necessary to compute all the eigenmodes obtained through the expansion.

It is important to note here that, to obtain more accurate results, we need to calculate the intensity distribution projected from each subregion onto the entire image plane, which means different sampling windows on the object plane and the image plane. If we use FFT for the calculations, we would need to perform zero-padding on the coherent light field modes of the object plane to ensure that their matrix sizes are the same as the sampling matrix size of the image plane.

By computing the light field distribution after the zero-padded light field matrices pass through the optical imaging system, we obtain the corresponding intensity distribution results of the partially coherent light fields from each independent subregion, as shown in Figure 3.

Next, let us calculate the coherence decomposition between each pair of subregions. SVD decompositions are performed on each submatrix that is not on the diagonal, decomposing ${\mathit{W}}_{pq}$ into two sets of weighted spatial coherent modes. Through this mathematical decomposition method, we can construct partial coherent light fields in two subregions, whose cross-correlation properties satisfy the mathematical description of the ${\mathit{W}}_{pq}$ matrix. Figure 4 shows the decomposed coherence modes between adjacent subregions ${\mathit{A}}_1$ and ${\mathit{A}}_2$, with the first 14 modes in each subregion. The Kohler illumination satisfies the property that the mutual coherence between different spatial positions decreases with increasing distance between these spatial positions. It could be found in Figure 4 that, in each coherent mode obtained through decomposition, the pixels closer to the adjacent subregions have larger amplitudes. This is consistent with theoretical expectations.

Similarly, Figure 5 illustrates the coherence between diagonal subregions ${\mathit{A}}_1$ and ${\mathit{A}}_4$ in Figure 1, with the first 14 modes of the two sets of partially coherent light fields decomposed from the mutual coherence of subregions ${\mathit{A}}_1$ and ${\mathit{A}}_4$. Similarly, we can also observe in Figure 5 that, among the various coherent modes obtained through decomposition, the pixels closer to the adjacent subregions exhibit larger amplitudes, similar to the former case.

By propagating the coherent modes decomposed from each subregion to the image plane, we can obtain the intensity distribution contributed by the mutual coherence between various subregions by referring to Equation (18). The image plane intensity contributed by each pair of mutual coherence is shown in Figure 6. It is worth noting that before calculating the propagation of the coherent modes of each subregion, we also need to perform zero-padding on these coherent modes.

As a final step, all calculated intensity distributions on the image plane are summed together. By following the description in Equation (9), the complete imaging intensity distribution could be obtained, as shown in Figure 7. In Figure 7, $I_1$ represents the sum of imaging results of the CSD on each subregion itself, while $I_2$ represents the sum of the contribution of the coherence between various subregions to the intensity distribution after propagating to the image plane. It could be found that the imaging result is imprecise if the mutual coherence between different subregions is not taken into account in partially coherent imaging cases. Moreover, three cases of different degrees of partial coherence are given in Figure 7. The first case (R = 1.5r) represents an overfill (1.5$$) of the lens pupil with the image of the incoherent primary source, which appears similar to incoherent cases. The second case (R = 0.7r) presents an apparent partially coherent imaging effect, while the last case (R = 0.1r) is almost indistinguishable from the fully coherent case.

Figure 8 gives a further illustration of the partially coherent effect. The black line in Fig shows the ideal geometrical image of the edge and the other curves are the diffracted images under different coherence conditions. The red curve (R = 1.5r) represents the overfill case and the yellow curve (R = 0.1r) is near the fully coherent case, while the blue curve (R = 0.7r) presents an apparent partially coherent imaging effect.

3.2. Performance in Comparison with Conventional Coherent Mode Decomposition Method

A comparison of performance between the traditional CMD algorithm and the newly proposed algorithm is given in this section. A group of identical partially coherent imaging parameters are set for the two simulation methods, as shown in

Table 1. Key parameters of the partially coherent imaging simulations under different settings.

Simulating Method	Wavelength (nm)	Numerical Aperture	Illumination Filling Factor	Object Scale (um)	Sampling Number of a Single Unit in the Space Domain	Sampling Number of a Single CSD Matrix	Number of Decomposed Coherent Modes
Conventional coherent mode decomposition	1	0.75	1	2 × 2	60$\times$60	3600$\times$3600	3600
Modified method (2$\times$2)					30$\times$30	900$\times$900	14,400
Modified method (3$\times$ 3)					20$\times$20	400$\times$400	32,400
Modified method (4$\times$4)					15$\times$15	225$\times$225	57,600
Modified method (5$\times$5)					12$\times$12	144$\times$144	90,000

. Similar to the previous section, the partially coherent illumination is assumed to be Kohler illumination. The source wavelength, numerical aperture, illumination filling factor, and spatial size of the object plane are present in Table 1. The simulation under several different subregion partitioning strategies, namely 2 × 2 subregions, 3 × 3 subregions, 4 × 4 subregions, and 5 × 5 subregions, are also demonstrated. The performance comparison between different partitioning strategies could help us understand how to decide the optimal strategy in different simulation works. Since a 60 × 60 discrete sampling was performed in the original area, the number of sampling points within each sub-region will decrease proportionally as 30, 20, 15, and 12. Correspondingly, the size of the CSD matrix of each subregion and between different subregions will also decrease rapidly. However, as the number of subregions increases, the number of coherent modes obtained from the decomposition will rapidly increase, as shown in Equations (13)−(15).

The performance comparison of the simulation methods will be presented from the following dimensions: the requirement for computer memory resources, the time consumption, and the calculation precision. The requirement of computational memory resources mainly arises from the coherent mode decomposition process, which decomposes the four-dimensional CSD matrix into a set of weighted spatial coherent light field modes. We will break down the simulation time consumption into the coherent decomposition phase (primarily SVD calculations) and the imaging simulation phase and, finally, provide the total time consumption. For the evaluation of calculation precision, the simulation results of the traditional coherent mode decomposition method were used as the benchmark. The results of the modified method under various partitioning strategies were compared with the benchmark and the proportion of error was given.

The simulation work was carried out using MATLAB R2021a software on a computer hardware platform with an I5 processor and 16GB of RAM. The comparison of simulation results is shown in

Table 2. Comparison of performance between the conventional method and the proposed method.

Simulating Method	Memory Usage of Coherent Mode Decomposition (MB)	Computation Time of SVD (s)	Computation Time of Imaging (s)	Total Computation Time (s)	Precision (Proportion of Error in Comparison with the Conventional Result)
Conventional coherent mode decomposition	648	46.47447	0.579663	47.05413	0
$\mathrm{Modified} \mathrm{method} (2 \times$2)	40.5	9.019392	6.084642	15.10403	$~1 \times {10}^{-16}$
$\mathrm{Modified} \mathrm{method} (3 \times$3)	8	2.5074331	13.66546	16.1729	$~1 \times {10}^{-16}$
$\mathrm{Modified} \mathrm{method} (4 \times$4)	2.53125	2.0448191	23.15931	25.20413	$~1 \times {10}^{-16}$
$\mathrm{Modified} \mathrm{method} (5 \times$5)	1.0368	2.024232	39.94231	41.96654	$~1 \times {10}^{-16}$

. As the number of subregions increases, the consumption of computer memory resources decreases significantly. At the same time, the time consumption in the SVD calculation phase also decreases rapidly with the increase in the number of subregions. However, the time consumption of the imaging simulation calculation phase increases with the increase in the number of subregions. Therefore, the total time consumption does not show a monotonic increase or decrease in trend with the increase in the number of subregions but, rather, exhibits first decreasing and then increasing patterns. In comparison of calculation precision, under different partitioning strategies, the calculation results are basically consistent with the calculation results of the traditional coherence mode decomposition method. The difference in simulation results mainly originates from the numerical calculation accuracy of the MATLAB software.

To explain the trend of simulation time changes in the table above, we have organized the overall simulation process as well as the calculation types and complexity of each step in Figure 9. Assume the original spatial domain sampling number is $N \times N$. In the modified method, the original spatial domain is divided into $M \times M$ subregions, with each subregion having a sampling number of $n \times n$. Obviously, we have $N=M \times n$.

In the coherent mode decomposition (CSD) process, the computational complexity of the singular value decomposition (SVD) operation for a single CSD matrix is $O(N^6)$, where the dimension of the CSD matrix is $N^2 \times N^2$. When the original space domain was divided into $M^2$ subregions, the SVD operation would be performed a total of $M^4$ times, with the time complexity of a single SVD operation being $O(n^6)$. Since the accelerating speed of $M^4$ is far slower than the decelerating of $n^6$, the overall time will exhibit a significant downward trend.

The simulation of imaging process primarily involves two Fourier transforms. We often utilize Fast Fourier Transform (FFT) to accelerate the numerical computation of the Fourier transform. For the traditional method, during the first step of modal decomposition, $N^2$ coherent modes were obtained and each mode would undergo an imaging calculation process. Therefore, the Fourier transforms need to be carried out by $2N^2$ times. However, during the modal decomposition of the modified method, we obtain a greater number of coherent modes ($M^4 \times n^2$, see Table 1). Similarly, each mode needs to go through the propagation process of the imaging system. Moreover, since we ultimately aim to obtain the intensity distribution within the entire spatial domain on the image plane, each Fourier transform must be conducted within the complete spatial domain, which means that the computational complexity of a single Fourier transform is the same as that of the traditional method. Thus, during the imaging calculation, the improved method requires executing several times more numerical calculations than the traditional method, resulting in significantly higher time consumption. Additionally, as the number of subregions increases, the time consumption will also rise rapidly.

Finally, we present a more explicit demonstration of the reduction in storage requirements with the proposed method, with an example shown in Figure 10, where the sampling number of complex amplitude on the object plane is 500$\times$500 (in single precision mode). If the entire plane is to be computed by the traditional coherent-mode method, 465 GB of memory is needed for the storage of the CSD matrix and another 976 GB of memory is needed for executing the matrix decomposition. If the object plane is split into several subregions, the storage requirements would decrease rapidly as the subregions become smaller.

4. Further Discussion

4.1. Discussion on the Potential Limitations of the Modified Modal Expansion Method

Firstly, the applicability of the modified modal expansion method is discussed. The method proposed in this paper is a purely numerical simulation approach; as long as the required input conditions for the numerical calculation method are met, the method is applicable. The input conditions referred to here are the CSD function or matrix of an incident partially coherent optical field. The CSD can be either a mathematical function with a known analytical expression or numerical results without a known analytical form (such as a measurement result). Therefore, the application of this method is not limited to certain specific types of optical fields. All partially coherent optical fields can be simulated and analyzed using this method.

Next, we discuss the performance limitations of the method. According to the demonstration in Section 3.2, as the number of subregions increases, the computational load for simulating the subsequent propagation of the optical field (such as through imaging systems or free space diffraction) gradually rises. For instance, when the number of subregions in the simulation example of Section 3.2 exceeds 5 × 5, the total time consumption surpasses that required by the traditional method. And the time consumption continues to increase with more subregions. This deserves particular attention in simulation tasks where the computation of subsequent optical field propagation is rather extensive. At this point, we should find an optimal balance between computer memory requirements and simulation time consumption. On the other hand, if the computational load of subsequent propagation is very small, further increasing the number of segmented regions might be more beneficial, which could reduce computation time due to the smaller matrix size processed by the SVD procedure. In summary, during the application of this method, we should focus on achieving the best balance between requirements on memory and simulation time consumption, by selecting the optimal number of segmented regions.

4.2. Applicability on Anisoplanatic Imaging Simulations

For a further acceleration of the propagation of modes (as Equation (8) or Equation (16), the imaging systems are always supposed to be isoplanatic [18], such that $K\left(\mathit{\rho },{\mathit{\rho }}_1^{\prime }\right)=K\left(\mathit{\rho }-{\mathit{\rho }}_1^{\prime }\right)$. Then, the integral could be rewritten in a convolution form, which could be accelerated by fast Fourier transformation (FFT). Actually, the assumption of an isoplanatic system is a compromise between computation efficiency and accuracy: in real imaging systems, the pupil function is more likely to be dependent on object coordinates. If we are interested in the imaging of a large object, the assumption of an isoplanatic system could introduce errors in the calculated image.

In our proposed method, the imaging systems are supposed to be isoplanatic in each subregion ${\mathit{A}}_p$, where the computational precision could be ensured if the object plane is split into small enough subregions for their isoplanatic nature. From Equations (10), (16), and (17) we found that the transfer functions corresponding to different subregions are allowed to be different and that the introduction of different transfer functions did not slow down the computation. Thus, if a partial coherent imaging system is to be modeled, we can split its object plane into several small areas, which are small enough to be nearly isoplanatic. Then, each subfield of the object field could be propagated with different transfer functions $K_n\left(\mathit{\rho }-{\mathit{\rho }}_1^{\prime }\right)$, as shown in Figure 11. The propagation of modes in Equations (16) and (17) could be rewritten as the following expression:

$$u_{npq}\left(\mathit{\rho }\right)={\displaystyle {\int }_{A_p}t\left({\mathit{\rho }}_1^{\prime }\right){\psi }_{npq}\left({\mathit{\rho }}_1^{\prime }\right)K_p\left(\mathit{\rho }-{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}$$

(19)

$$v_{npq}\left(\mathit{\rho }\right)={\displaystyle {\int }_{A_q}t\left({\mathit{\rho }}_1^{\prime }\right){\zeta }_{npq}\left({\mathit{\rho }}_1^{\prime }\right)K_q\left(\mathit{\rho }-{\mathit{\rho }}_1^{\prime }\right){\mathbf{d}}^2{\mathit{\rho }}_1^{\prime }}$$

(20)

5. Conclusions

In conclusion, we have shown that the numerical computation of partially coherent imaging could be carried out through the propagation of subfields in the object plane. The computation of subfields under partially coherent illumination, compared with the direct handling of the complete random optical field, significantly reduces the requirement of computational resources. With the newly proposed method, we are able to carry out the partially coherent imaging computation of arbitrary large sampling numbers, with an ordinary computer. The performance between the newly proposed method and the traditional CMD method is also given in detail, in terms of the storage requirement, time consumption, and precision. From the experimental results and the corresponding theoretical analysis, we could find that the increase in the number of subregions mainly affects the computational load in subsequent propagation stages. Therefore, we need to consider both the reduced memory requirements brought by the increase in sub-regions and the increase in simulation time. We should find the optimal balance point in different application scenarios, which might be the future work of this research. Moreover, we have also shown that the application of the newly proposed method is not restricted to some special model of partial coherence, nor the assumption of the isoplanatic nature of the imaging system.