The recent resurgence of analog optical information processing has been spurred by advancements in artificial intelligence (AI), especially deep-learning-based inference methods.1^–9 These advances in data-driven learning methods have also benefited optical hardware engineering, giving rise to new computing architectures such as diffractive deep neural networks (${\mathrm{D}}^2\mathrm{NN}$), which exploit the passive interaction of light with spatially engineered surfaces to perform visual information processing. ${\mathrm{D}}^2\mathrm{NNs}$, also referred to as diffractive optical networks, diffractive networks, or diffractive processors, have emerged as powerful all-optical processors9^,10 capable of completing various visual computing tasks at the speed of light propagation through thin passive optical devices; examples of such tasks include image classification,11^–13 information encryption,14^–17 and quantitative phase imaging (QPI),18^,19 among others.20^–24 Diffractive optical networks comprise a set of spatially engineered surfaces, the transmission (and/or reflection) profiles of which are optimized using machine-learning techniques. After their digital optimization (a one-time effort), these diffractive surfaces are fabricated and assembled in 3D to form an all-optical visual processor, which axially extends at most a few hundred wavelengths ($\lambda$).

Our earlier work10^,25 demonstrated that a spatially coherent ${\mathrm{D}}^2\mathrm{NN}$ can perform arbitrary complex-valued linear transformations between a pair of arbitrary input and output apertures if its design has a sufficient number ($N$) of diffractive features that are optimized, i.e., $N\ge N_i{N}_o$, where $N_i$ and $N_o$ represent the space-bandwidth product of the input and output apertures, respectively. In other words, $N_i$ and $N_o$ represent the size of the desired complex-valued linear transformation $\mathit{A}\in {\mathbb{C}}^{{\mathit{N}}_{\mathit{o}}\times {\mathit{N}}_{\mathit{i}}}$ that can be all-optically performed by an optimized ${\mathrm{D}}^2\mathrm{NN}$. For a phase-only diffractive network, i.e., only the phase profile of each diffractive layer is trainable, the sufficient condition becomes $N\ge 2N_i{N}_o$ due to the reduced degrees of freedom within the diffractive volume. Similar conclusions can be reached for a diffractive network that operates under spatially incoherent illumination: Rahman et al.26 demonstrated that a diffractive network can be optimized to perform an arbitrary nonnegative linear transformation of optical intensity through phase-only diffractive processors with $N\ge 2N_i{N}_o$. However, encoding information with spatially incoherent light inherently confines both the input and output to nonnegative values, as they are represented by intensity patterns at the input and output apertures of a ${\mathrm{D}}^2\mathrm{NN}$. To process complex-valued data with spatially incoherent light, other optical approaches were also developed;1^,27^–29 however, these earlier systems are limited to one-dimensional (1D) optical inputs and do not cover arbitrary input and output apertures, limiting their functionality and processing throughput. An extension of these earlier 1D input approaches introduced the processing of 2D incoherent source arrays using relatively bulky and demanding optical projection systems that are hard to operate at the diffraction limit of light.30^,31

Here, we demonstrate the processing of complex-valued data with compact diffractive optical networks under spatially incoherent illumination. We show that a spatially incoherent diffractive network that axially spans $<100\times \lambda$ can perform any arbitrary complex-valued linear transformation on complex-valued input data with negligible error if the number of optimizable diffractive features is above a threshold dictated by the multiplication of the input and output space-bandwidth products, determined by both the spatial extent and the pixel size of the input and output apertures. To represent complex-valued spatial information using spatially incoherent illumination, we preprocessed the input information by mapping complex-valued data to a real and nonnegative, optical intensity-based representation at the input field of view (FOV) of the diffractive network. We term this mapping the “mosaicking” operation, indicating the utilization of multiple intensity pixels at the input FOV to represent one complex-valued input data point. Similarly, we used a postprocessing step, which involved mapping the output FOV intensity patterns back to the complex number domain, which we termed the “demosaicking” operation. Through these mosaicking/demosaicking operations, we show that a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ can be optimized to perform an arbitrary complex-valued linear transformation between its input and output apertures while providing optical information encryption. The presented spatially incoherent visual information processor, with its universality and thin form factor ($<100\times \lambda$), shows significant promise for image encryption and computational imaging applications under natural light.

Figure 1(a) outlines a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ architecture to synthesize an arbitrary complex-valued linear transformation ($\mathit{A}$) such that $\mathit{o}=\mathit{A}\mathit{i}$, where the input is $\mathit{i}\in {\mathbb{C}}^{N_i}$, the target is $\mathit{o}\in {\mathbb{C}}^{N_o}$ and $\mathit{A}\in {\mathbb{C}}^{N_o\times N_i}$. The mosaicking process involves finding the nonnegative (optical intensity-based) representation of each complex-valued element of $\mathit{i}$ using $E$ nonnegative values; here, $E$ bases, $e_k$, $k=0,\cdots ,E-1$ [see Fig. 1(c)], are used for representing the intensity-based encoding of complex numbers. Based on this representation, the 2D input aperture of a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ will have $E{N}_i$ nonnegative (optical intensity) values, denoted as ${\mathit{i}}_{\mathit{r}}\in {\mathbb{R}}_{+}^{E{N}_i}$, representing the input information under spatially incoherent illumination. The output intensity distribution, denoted with $\widehat{{\mathit{o}}_{\mathit{r}}}\in {\mathbb{R}}_{+}^{E{N}_o}$, undergoes a demosaicking process where a complex number is synthesized from the intensity values of $E$ output pixels, yielding the complex output vector $\hat{\mathit{o}}\in {\mathbb{C}}^{N_o}$ such that $\hat{\mathit{o}}\approx \mathit{A}\mathit{i}$.

In our analyses, we used $E=3$, except in Fig. S5 in the Supplementary Material, where $E=4$ results are shown for comparison. We chose the basis complex numbers as $e_k=\exp (jk\frac{2\pi }{E})$, $k=0,\cdots ,E-1$ such that the set of bases $\mathcal{S}$ is closed under multiplication, and the product of any two of the bases in the set is also a basis; for example, for $E=3$ we have $e_k{e}_l=e_{(k+l\mathrm{mod}3)}$. Based on this representation of information, with $E=3$ and $e_0$, $e_1,e_2$, we can decompose any arbitrarily selected complex-valued transformation matrix $\mathit{A}$ into $E=3$ matrices (${\mathit{A}}_0$, ${\mathit{A}}_1$, ${\mathit{A}}_2$) with real nonnegative entries such that $$\mathit{A}=e_0{\mathit{A}}_0+e_1{\mathit{A}}_1+e_2{\mathit{A}}_2.$$(1)

For a given complex-valued input $\mathit{i}=e_0{\mathit{i}}_0+e_1{\mathit{i}}_1+e_2{\mathit{i}}_2$, where ${\mathit{i}}_{\mathit{k}}\in {\mathbb{R}}_{+}$, the corresponding target output vector can be written as $$\mathit{o}=\mathit{A}\mathit{i}=(e_0{\mathit{A}}_0+e_1{\mathit{A}}_1+e_2{\mathit{A}}_2)(e_0{\mathit{i}}_0+e_1{\mathit{i}}_1+e_2{\mathit{i}}_2),$$(2)$$\mathit{o}=e_0({\mathit{A}}_0{\mathit{i}}_0+{\mathit{A}}_2{\mathit{i}}_1+{\mathit{A}}_1{\mathit{i}}_2)+e_1({\mathit{A}}_1{\mathit{i}}_0+{\mathit{A}}_0{\mathit{i}}_1+{\mathit{A}}_2{\mathit{i}}_2)+e_2({\mathit{A}}_2{\mathit{i}}_0+{\mathit{A}}_1{\mathit{i}}_1+{\mathit{A}}_0{\mathit{i}}_2),$$(3)i.e., we have $${\mathit{o}}_{\mathit{r}}=\left[\begin{array}{c}{\mathit{o}}_0\\ {\mathit{o}}_1\\ {\mathit{o}}_2\end{array}\right]=\left[\begin{array}{ccc}{\mathit{A}}_0& {\mathit{A}}_2& {\mathit{A}}_1\\ {\mathit{A}}_1& {\mathit{A}}_0& {\mathit{A}}_2\\ {\mathit{A}}_2& {\mathit{A}}_1& {\mathit{A}}_0\end{array}\right]\left[\begin{array}{c}{\mathit{i}}_0\\ {\mathit{i}}_1\\ {\mathit{i}}_2\end{array}\right]={\mathit{A}}_{\mathit{r}}{\mathit{i}}_{\mathit{r}},$$(4)with a nonnegative real-valued matrix ${\mathit{A}}_{\mathit{r}}$$${\mathit{A}}_{\mathit{r}}=\left[\begin{array}{ccc}{\mathit{A}}_0& {\mathit{A}}_2& {\mathit{A}}_1\\ {\mathit{A}}_1& {\mathit{A}}_0& {\mathit{A}}_2\\ {\mathit{A}}_2& {\mathit{A}}_1& {\mathit{A}}_0\end{array}\right].$$(5)

For $E=4$, where $e_k{e}_l=e_{(k+l\mathrm{mod}4)}$ and $\mathit{A}=e_0{\mathit{A}}_0+e_1{\mathit{A}}_1+e_2{\mathit{A}}_2+e_3{\mathit{A}}_3$, a similar analysis yields $${\mathit{A}}_{\mathit{r}}=\left[\begin{array}{cccc}{\mathit{A}}_0& {\mathit{A}}_2& {\mathit{A}}_3& {\mathit{A}}_1\\ {\mathit{A}}_2& {\mathit{A}}_0& {\mathit{A}}_1& {\mathit{A}}_3\\ {\mathit{A}}_1& {\mathit{A}}_3& {\mathit{A}}_0& {\mathit{A}}_2\\ {\mathit{A}}_3& {\mathit{A}}_1& {\mathit{A}}_2& {\mathit{A}}_0\end{array}\right].$$(6)

Based on these equations, one can conclude that to all-optically implement an arbitrary complex-valued transformation $\mathit{o}=\mathit{A}\mathit{i}$ using a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$, the layers of the ${\mathrm{D}}^2\mathrm{NN}$ need to be optimized to perform an intensity linear transformation ${\mathit{A}}_{\mathit{r}}\in {\mathbb{R}}_{+}^{E^2{N}_i{N}_o}$ such that ${\mathit{o}}_{\mathit{r}}={\mathit{A}}_{\mathit{r}}{\mathit{i}}_{\mathit{r}}$. The entire system, upon convergence, performs the predefined complex-valued linear transformation $\mathit{A}$ on any given input data using spatially incoherent light, based on Eqs. (2) and (4). In the following sections, we numerically explore the number of optimizable diffractive features ($N$) needed for accurate approximation of $\mathit{A}$ using a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$.

We numerically demonstrated the capabilities of diffractive optical processors to universally perform any arbitrarily chosen complex-valued linear transformation with spatially incoherent light. Throughout the paper, we used $N_i=N_o=16$. To visually represent the data, we rearranged the 16-element vectors into $4\times 4$ arrays of complex numbers, hereafter referred to as the “complex image.” We arbitrarily selected a desired complex-valued transformation, $\mathit{A}\in {\mathbb{C}}^{16\times 16}$, as shown in Fig. 1(b).

To explore the number of diffractive features needed, we trained nine models with varying values of $N$ and evaluated the mean-squared-error (MSE) between the numerically measured ($\widehat{{\mathit{A}}_{\mathit{r}}}$) and the target all-optical linear transformation, ${\mathit{A}}_{\mathit{r}}$ (see Fig. 2). Our results, summarized in Fig. 2, highlight that with a sufficient number of optimizable diffractive features, i.e., $N\ge 2E^2{N}_i{N}_o=2N_{i,r}{N}_{o,r}$, our system achieves a negligible approximation error with respect to the target ${\mathit{A}}_{\mathit{r}}\in {\mathbb{R}}_{+}^{48\times 48}$. In Fig. 2(c), we also visualize the resulting all-optical intensity transformation $\widehat{{\mathit{A}}_{\mathit{r}}}$ compared to the ground truth ${\mathit{A}}_{\mathit{r}}$. In essence, this comparison reveals the spatially varying incoherent point spread functions (PSFs) of our diffractive system optimized using deep learning; a negligible MSE between $\widehat{{\mathit{A}}_{\mathit{r}}}$ and ${\mathit{A}}_{\mathit{r}}$ shows that the resulting spatially varying incoherent PSFs match the target set of PSFs dictated by ${\mathit{A}}_{\mathit{r}}$.

We also evaluated the numerical accuracy of our complex-valued transformation in an end-to-end manner, as illustrated in Fig. 2(d). For this numerical test, we sequentially set each entry of $\mathit{i}$ to $e_0$, evaluated the corresponding complex output $\hat{\mathit{o}}$, and stacked them to form $\widehat{{\mathit{A}}_0}$, where the subscript represents that the measurement was evaluated using the complex impulse along the basis $e_0$ as input. Then, we repeated this process for the other two bases to obtain $\widehat{{\mathit{A}}_1}$ and $\widehat{{\mathit{A}}_2}$, and stacked these matrices as a block matrix $[\widehat{{\mathit{A}}_0}|\widehat{{\mathit{A}}_1}|\widehat{{\mathit{A}}_2}]$, shown in Fig. 2(d). Each row of the images $\mathrm{amp}(\hat{\mathit{o}})$ and $\text{phase}(\hat{\mathit{o}})$ in Fig. 2(d) represents one of these complex output vectors, while the corresponding target vectors are presented in the same figure through $\mathrm{amp}(\mathit{o})$ and $\text{phase}(\mathit{o})$. The small magnitude of the error $\mathit{\epsilon }={|\hat{\mathit{o}}-\mathit{o}|}^2$ shown in Fig. 2(d) illustrates the success of this spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ model in accurately approximating the complex-valued linear transformation $\mathit{o}=\mathit{A}\mathit{i}$, implemented for an arbitrarily selected $\mathit{A}$.

In this section, we demonstrate a complex number-based image encryption–decryption scheme using a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$. In the first scheme, shown in Fig. 1(d), the message is encoded into a complex image, employing either amplitude and phase encoding or real and imaginary part encoding. Then, a digital lock encrypts the image by applying a linear transformation (${\mathit{A}}^{-1}$) to conceal the original message within the image. At the optical receiver, the encrypted message is deciphered by an optimized incoherent ${\mathrm{D}}^2\mathrm{NN}$ that all-optically implements the inverse transformation $\mathit{A}$. In an alternative scheme, as depicted in Fig. 1(e), the key and lock are switched, i.e., the spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ is used to encrypt the message with a complex-valued $\mathit{A}$ while the decryption step involves the digital inversion using ${\mathit{A}}^{-1}$.

For our analysis, we used the letters “U,” “C,” “L,” and “A” as sample messages. “U” and “C” are used in amplitude-phase-based encoding (Fig. 3), whereas “L” and “A” are used for real-imaginary-based encoding of information (Fig. S1 in the Supplementary Material), forming complex-number-based images. To accurately model the spatially incoherent propagation26 of light through the ${\mathrm{D}}^2\mathrm{NN}$, we averaged the output intensities over a large number of $N_{\phi }=\mathrm{20,000}$ of randomly generated 2D phase profiles at the input (see Sec. 4 for details).

In Fig. 3(a), we show the results corresponding to digital encryption and optical diffractive decryption, i.e., the system shown in Fig. 1(d). The digitally encrypted complex information $\mathit{i}={\mathit{A}}^{-1}\mathit{o}$, and its intensity representation ${\mathit{i}}_{\mathit{r}}$, are shown in Fig. 3(a). The optically decrypted output $\hat{\mathit{o}}$ (through the spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$) and its intensity-based representation $\widehat{{\mathit{o}}_{\mathit{r}}}$ are shown in the same Fig. 3(a), together with the resulting error maps, i.e., ${|\hat{\mathit{o}}-\mathit{o}|}^2$ and ${|\widehat{{\mathit{o}}_{\mathit{r}}}-{\mathit{o}}_{\mathit{r}}|}^2$, which reveal a very small degree of error. This agreement of the recovered and the ground-truth messages in both the intensity and complex-valued domains confirms the accuracy of the diffractive decryption process through an optimized spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$. Figure 3(b) shows the successful performance of the sister scheme [Fig. 1(e)], which involves diffractive encryption through a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$ and digital decryption, also revealing a negligible amount of error in both ${|{\mathit{A}}^{-1}\hat{\mathit{o}}-\mathit{i}|}^2$ and ${|\widehat{{\mathit{o}}_{\mathit{r}}}-{\mathit{o}}_{\mathit{r}}|}^2$. As reported in Fig. S1 in the Supplementary Material, we also conducted a numerical experiment using the letters “L” and “A,” encoded using the real and imaginary parts of the message. The visualizations are arranged the same way as in Fig. 3, where for both schemes depicted in Figs. 1(d) and 1(e), the degree of error between the recovered and the original messages is negligible, affirming the success of using the real and imaginary part-based encoding method. For the assessment of the approximation errors when the number of diffractive features is smaller, we compared the decryption performance of three models with different numbers of diffractive features/neurons, i.e., $N=(0.5,0.7,2)\times 2E^2{N}_i{N}_o$, for the same setup outlined in Fig. S1(a) in the Supplementary Material. The results are summarized in Fig. S2 in the Supplementary Material: for models with $N<2E^2{N}_i{N}_o$, the decryption quality is compromised, exhibiting a pixel absolute error of $>0.1$. However, this error reduces to $<0.05$ for $N=4E^2{N}_i{N}_o$ where the decrypted images display significantly enhanced contrast and reduced noise levels.

To further evaluate the efficacy of our encryption method, we analyzed the complex image entropy, examining both the real and imaginary components separately (refer to the Sec. 4 for details). The original image $\mathit{i}$, ${\mathrm{D}}^2\mathrm{NN}$ encrypted output $\hat{\mathit{o}}$, and the digitally encrypted output $\mathit{A}\mathit{i}$, along with the corresponding image entropies, are shown in Fig. S3(a) in the Supplementary Material for two complex image examples. We repeated this analysis for a set of 1000 complex images with the resulting entropy distributions reported in Fig. S3(b) in the Supplementary Material. These results demonstrate that the entropy of the encrypted images is statistically higher than that of the original images. This increase in entropy signifies a heightened level of randomness within the encrypted images, thereby validating the effectiveness of our encryption process. In addition, the entropy distributions of the ${\mathrm{D}}^2\mathrm{NN}$ encrypted images show excellent agreement with the digitally encrypted corresponding images, further demonstrating the success of our spatially incoherent optical encryption scheme.

How we assign each element in the vector ${\mathit{i}}_{\mathit{r}}$ and ${\mathit{o}}_{\mathit{r}}$ to the pixels at the input and output FOVs of the diffractive network does not affect the final accuracy of the image/message reconstruction. For example, we can arrange the FOVs in such a manner that the components ${\mathit{i}}_{r,k}$ corresponding to a basis $e_k$ are assigned to the neighboring pixels, in two adjacent rows, as shown in Fig. S4(a) in the Supplementary Material; in an alternative implementation, the assignment/mapping can be completely arbitrary, which is equivalent to applying a random permutation operation on the input and output vectors (see Sec. 4). When compared to each other, these two approaches of mosaicking and demosaicking schemes show negligible differences in the error of the final reconstruction of the letters “U,” “C,” “L,” and “A” as shown in Fig. S4(b) in the Supplementary Material. These results underscore that the specific arrangement of the mosaicking/demosaicking schemes at the input and output FOVs does not impact the performance of the incoherent ${\mathrm{D}}^2\mathrm{NN}$ system.

In this article, we employed a data-free PSF-based ${\mathrm{D}}^2\mathrm{NN}$ optimization method (see Sec. 4),26 since we can determine the nonnegative intensity transformation ${\mathit{A}}_{\mathit{r}}$ from the target complex-valued transformation $\mathit{A}$ based on the mosaicking and demosaicking schemes; the columns of ${\mathit{A}}_{\mathit{r}}$ represent the desired spatially varying PSFs of the ${\mathrm{D}}^2\mathrm{NN}$. The advantage of this data-free learning-based ${\mathrm{D}}^2\mathrm{NN}$ optimization approach is that computationally demanding simulation of wave propagation with large $N_{\phi }$ is not required during the training. Coherent propagation is appropriate for simulating the spatially varying PSFs, point by point, since a point emitter at the input aperture coherently interferes with itself during optical diffraction within a ${\mathrm{D}}^2\mathrm{NN}$; this approach makes the training time much shorter. On the other hand, this approach necessitates prior knowledge of ${\mathit{A}}_{\mathit{r}}$, which might not always be available, e.g., for tasks such as data classification. An alternative to this data-free PSF-based optimization approach is to train the diffractive network in an end-to-end manner, using a data-driven direct training approach.26 This strategy advances by minimizing the differences between the outputs and the targets on a large number of randomly generated examples, thereby learning the spatially varying PSFs implicitly from numerous input-target intensity patterns corresponding to the desired task – instead of learning from an explicitly predetermined ${\mathit{A}}_{\mathit{r}}$. This direct approach, however, requires a longer training time, necessitating the simulation of incoherent propagation for each training sample on a large data set.

In our presented approach, the choice of $E$ is not restricted to $E=3$, as we have used throughout the main text. As another example of encoding, we show the image encryption results with $E=4$ in Fig. S5 in the Supplementary Material, where the four bases are $\exp (j\frac{\pi }{2}k)(k=\mathrm{0,1},\mathrm{2,3})$. The reconstructed “U,” “C,” “L,” and “A” letters are also reported in the same figure, confirming that given sufficient degrees of freedom (with $N\ge 2E^2{N}_i{N}_o$), the linear transformation performances are similar to each other. However, compared to $E=3$, this choice of $E=4$ necessitates 4/3 times more pixels on both the diffractive network input and output FOVs—reducing the throughput (or spatial density) of complex-valued linear transformations that can be performed using a spatially incoherent ${\mathrm{D}}^2\mathrm{NN}$. Accordingly, more diffractive features and a larger number of independent degrees of freedom (by 16/9-fold) are required within the ${\mathrm{D}}^2\mathrm{NN}$ volume to achieve an output performance level that is comparable to a design with $E=3$. Note that while $E\ge 3$ is sufficient to reconstruct the original complex-valued images regardless of the image complexity, the redundancy provided by larger $E$ values might offer increased resilience against noise at the cost of reducing the image-processing throughput (per input aperture area) with larger $E$.

Our framework offers several flexibilities in implementation, which could be useful for different applications. First, the flexibility to arbitrarily permute the input and the output pixels following different mosaicking and demosaicking schemes (as introduced earlier in Sec. 2) could enhance the security of optical information transmission. A user would not be able to either spam or hack valuable information that is transferred optically without specific knowledge of the mosaicking and demosaicking schemes, thus ensuring the security of this scheme. Note that this enhancement in security is achieved without adding complexity to the system by just permuting the assignment of data elements to the pixels of the input and output devices, e.g., spatial light modulators (SLMs) and complementary metal-oxide-semiconductor (CMOS) detector-arrays. Second, the flexibility in choosing $E$, as discussed above, could be useful in adding an extra layer of security against unauthorized access, albeit with a trade-off in system throughput that comes with larger $E$. Furthermore, we can use different sets of bases for mosaicking and demosaicking operations by applying offset phase angles ${\theta }_i$ and ${\theta }_o$, respectively, to the original bases $e_k=\exp (jk\frac{2\pi }{E})$, $k=0,\cdots ,E-1$. This will result in a set of modified/encrypted bases: $e_{k,i}=\exp [j(k\frac{2\pi }{E}+{\theta }_i)]$ for mosaicking and $e_{k,o}=\exp [j(k\frac{2\pi }{E}+{\theta }_o)]$ for demosaicking. This powerful flexibility in representation further enhances the security of the system.

Regarding image encryption-related applications, we demonstrated two approaches [Figs. 1(d) and 1(e)] to utilize ${\mathrm{D}}^2\mathrm{NNs}$ for encryption or decryption. However, it is also possible to deploy a pair of diffractive systems in tandem, with one undertaking the matrix operation $\mathit{A}$ for encryption and the other undertaking the inverse operation ${\mathit{A}}^{-1}$ for decryption. Furthermore, potential extensions of our work could explore a harmonized integration of polarization state controls32 and wavelength multiplexing33 to build a multifaceted, fortified encryption platform. In addition to increasing the data throughput, these additional degrees of freedom enabled by different illumination wavelengths and polarization states would further enhance the security of a diffractive processor-based system.

In this work, we focused on the numerical analysis of the presented concept. However, we should note that various ${\mathrm{D}}^2\mathrm{NNs}$ designed using deep-learning-based approaches have been experimentally validated over different parts of the electromagnetic spectrum, e.g., from terahertz (THz)9^,14 to near-infrared (NIR)15 and visible wavelengths,24 showing a good agreement between the numerical and experimental results. To address some of the experimental challenges associated with fabrication errors and mechanical misalignments, a “vaccination” strategy34^,35 has been devised. This approach enhances the robustness of the diffractive optical designs by incorporating such aberrations/imperfections as random variables during the training phase, thereby preparing the system to better withstand and adapt to the uncertainties inherent in real-world experimental conditions.

Although spatially coherent light is more suitable for complex-valued information processing in laboratory settings, the use of spatially incoherent light offers various practical advantages. For example, speckle noise, which is inevitable in coherent systems, can be suppressed by using partially or fully incoherent illumination. An additional benefit of spatially incoherent designs is the range of viable illumination sources that can be used: instead of using specialized coherent sources, a spatially incoherent system can work with standard light-emitting diodes (LEDs), or even under natural light, which is important for some applications of diffractive information processing.

To conclude, we demonstrated the capability of spatially incoherent diffractive networks to perform arbitrary complex-valued linear transformations. By incorporating various forms of mosaicking and demosaicking operations, we paved the way for a wider array of applications by leveraging incoherent ${\mathrm{D}}^2\mathrm{NNs}$ for complex-valued data processing. We also showcased potential applications of these spatially incoherent ${\mathrm{D}}^2\mathrm{NNs}$ for complex number-based image encryption or decryption, highlighting the security benefits arising from the system’s flexibility. Our exploration marks a significant stride toward enhanced versatility and robustness in optical information processing with spatially incoherent diffractive systems that can work under natural light.

In this paper, we use $N_i=N_o=16$ so that $\mathit{A}\in {\mathbb{C}}^{16\times 16}$; see Fig. 1(b). To generate $\mathit{A}$, we randomly sample the amplitude of each element from the uniform distribution $Uniform(\mathrm{0,}1)$ and the phases from $Uniform(\mathrm{0,}2\pi )$. For the encryption application, to ensure that the result of inversion is not sensitive to small errors, we performed QR-factorization on $\mathit{A}$ to obtain a condition number of one.36

Following Eq. (4), the complex-valued input and target vectors $\mathit{i}\in {\mathbb{C}}^{N_i}$ and $\mathit{o}\in {\mathbb{C}}^{N_o}$ are represented by the corresponding real and nonnegative intensity vectors ${\mathit{i}}_{\mathit{r}}=[\begin{array}{ccc}{\mathit{i}}_0^T& \cdots & {\mathit{i}}_{\mathit{E}-1}^T\end{array}{]}^T\in {\mathbb{R}}_{+}^{E{N}_i}$ and ${\mathit{o}}_{\mathit{r}}=[\begin{array}{ccc}{\mathit{o}}_0^T& \cdots & {\mathit{o}}_{\mathit{E}-1}^T\end{array}{]}^T\in {\mathbb{R}}_{+}^{E{N}_o}$, where $\mathit{i}=\sum _{k=0}^{E-1}{e}_k{\mathit{i}}_k$ and $\mathit{o}=\sum _{k=0}^{E-1}{e}_k{\mathit{o}}_k$. The desired all-optical intensity transformation ${\mathit{A}}_r$ between ${\mathit{i}}_{\mathit{r}}$ and ${\mathit{o}}_{\mathit{r}}$ is derived from the target complex-valued linear transformation $\mathit{A}$ following Eqs. (1) and (5). We should note that deriving ${\mathit{A}}_r$ from $\mathit{A}$ requires mapping each complex element $a$ to its real and nonnegative representation $(a_0,\cdots ,a_{E-1})$ based on the $E\ge 3$ complex bases $e_k$ such that $a=\sum _{k=0}^{E-1}{e}_k{a}_k$. To define a unique mapping, we follow an algorithm29 by imposing additional constraints: $a_k=0$ if $\frac{2\pi }{E}\le \text{phase}(e_k{a}^{*})\le 2\pi -\frac{2\pi }{E}$, i.e., $a_k=0$ if the angle between $a$ and $e_k$ is greater than $\frac{2\pi }{E}$; here $a^{*}$ represents the complex conjugate of $a$. The same constraints were also used while mapping the complex input vectors $\mathit{i}$ to the real and nonnegative intensity vectors ${\mathit{i}}_r$.

For mosaicking (demosaicking) assignment of each element of ${\mathit{i}}_{\mathit{r}}$ (${\mathit{o}}_{\mathit{r}}$) to one of the $N_{i,r}=E{N}_i$ ($N_{o,r}=E{N}_o$) pixels of the 2D input (output), the arrangement of the FOV can be regular, e.g., in a row-major order as shown in Fig. S4(a) in the Supplementary Material, “Regular mosaicking.” Alternatively, the pixel assignment on the input (output) FOV can follow any arbitrary mapping which can be defined by a permutation matrix ${\mathit{P}}_{\mathit{i}}$ (${\mathit{P}}_{\mathit{o}}$) operating on the input (output) vector; see Fig. S4(a) in the Supplementary Material, “Arbitrary mosaicking.” For such cases, when ordered in a row-major format, intensities on the input (output) FOVs ${\mathit{i}}_{\mathit{r}}$ (${\mathit{o}}_{\mathit{r}}$) can be written as ${\mathit{i}}_{\mathit{r}}={\mathit{P}}_{\mathit{i}}[\begin{array}{ccc}{\mathit{i}}_0^T& \cdots & {\mathit{i}}_{\mathit{E}-1}^T\end{array}{]}^T$ (${\mathit{o}}_{\mathit{r}}={\mathit{P}}_{\mathit{o}}[\begin{array}{ccc}{\mathit{o}}_0^T& \cdots & {\mathit{o}}_{\mathit{E}-1}^T\end{array}{]}^T$). Accordingly, such an arbitrary arrangement of pixels was accounted for by redefining the all-optical intensity transformation as ${\mathit{P}}_{\mathit{o}}{\mathit{A}}_{\mathit{r}}{\mathit{P}}_{\mathit{i}}^{\mathit{T}}$.

The 1D vector ${\mathit{i}}_{\mathit{r}}$ is rearranged as a 2D distribution of intensity $I(x,y)$ at the input FOV of the ${\mathrm{D}}^2\mathrm{NN}$. To numerically model the spatially incoherent propagation of the input intensity distribution $I(x,y)$ through the ${\mathrm{D}}^2\mathrm{NN}$, we coherently propagated the optical field $\sqrt{I}\exp (j\phi )$ through the trainable diffractive surfaces to the output plane, where $\phi$ is a random 2D phase distribution, i.e., $\phi (x,y)\sim Uniform(0,2\pi )$ for all $(x,y)$. If we denote the coherent field propagation operator as $\mathfrak{D}\{·\}$ (see Sec. 4.5), then the instantaneous output intensity is $|\mathfrak{D}\{\sqrt{I(x,y)}\exp [j\phi (x,y)]\}{|}^2$ and the time-averaged output intensity $O(x,y)$ for spatially incoherent light can be written as $$O(x,y)=⟨{\left|\mathfrak{D}\right\{\sqrt{I(x,y)}\exp [j\phi (x,y)]|}^2⟩.$$(7)

The average output intensity can be approximately calculated by repeating the coherent wave propagation $\mathfrak{D}\{·\}$$N_{\phi }$-times, each time with a different random phase distribution ${\phi }_r(x,y)$, and averaging the resulting $N_{\phi }$ output intensities, $$O(x,y)=\underset{N_{\phi }\to \infty }{\lim }\frac{1}{N_{\phi }}\sum _{r=1}^{N_{\phi }}{\left|\mathfrak{D}\right\{\sqrt{I(x,y)}\exp [j{\phi }_r(x,y)]\left\}\right|}^2.$$(8)

We used $N_{\phi }=\mathrm{20,000}$ for estimating the incoherent output intensity $O(x,y)$ corresponding to any arbitrary input intensity $I(x,y)$. Note that when only one pixel at the input aperture is activated, with all other input pixels being inactive with zero intensity, as is the case while evaluating spatially varying PSFs, the application of Eq. (8) becomes redundant, although one could still use it. In this scenario, all the light diffracted from a single point source is mutually coherent. Consequently, for the purposes of evaluating the spatially varying PSFs of the system, as elaborated later in Sec. 4.7, employing a coherent propagation model for each point emitter at the input aperture is accurate and provides a faster solution.

The propagation of spatially coherent light patterns through a diffractive processor, denoted by $\mathfrak{D}\{·\}$, involves a series of interactions with consecutive diffractive surfaces, interleaved by wave propagation through the free space separating these surfaces. We assume that these modulations are introduced by phase-only diffractive surfaces, i.e., the field amplitude remains unchanged during the light–matter interaction. Specifically, we assume that a diffractive surface alters the incident optical field, symbolized as $u(x,y)$, in a localized manner according to the optimized phase values ${\varphi }_M(x,y)$ of the diffractive features, resulting in the phase-modulated field $u(x,y)\exp [j{\phi }_M(x,y)]$. The diffractive surfaces are coupled by free-space propagation, allowing the light to travel from one surface to the next. We used the angular spectrum method to simulate the free-space propagation,37$$u(x,y;z=z_0+d)={\mathcal{F}}^{-1}\{\mathcal{F}\{u(x,y;z=z_0)\}\times H(f_x,f_y;d)\},$$(9)where $\mathcal{F}\{·\}$ is the 2D Fourier transform and ${\mathcal{F}}^{-1}\{·\}$ is its inverse operation. $H(f_x,f_y;d)$ is the free-space transfer function corresponding to propagation distance $d$. For wavelength $\lambda$, $$H(f_x,f_y;d)=\{\begin{array}{ll}\exp \left[j\frac{2\pi }{\lambda }d\sqrt{1-{(\lambda f_x)}^2-{(\lambda f_y)}^2}\right],& f_x^2+f_y^2<1/{\lambda }^2\\ 0,& \text{otherwise}\end{array}.$$(10)

The fields were discretized with a lateral sampling interval of $\delta \approx 0.53\lambda$ to accommodate all the propagating modes and sufficiently zero-padded to remove aliasing artifacts.38

We modeled the diffractive surfaces by their laterally discretized heights $h$, which correspond to phase delays ${\phi }_M=\frac{2\pi }{\lambda }(n-1)h$, where $n$ is the refractive index of the material. The connectivity between consecutive diffractive layers9 was kept equal across the diffractive designs with varying $N$ by setting the separation between the layers as $d=\frac{W\delta }{\lambda }$, where the width of each diffractive layer is $W=\sqrt{\frac{N}{K}}\delta$. Here, $K$ is the number of diffractive layers; we used $K=4$ throughout the paper.

For performing an arbitrary complex-valued linear transformation with a diffractive processor, we used the PSF-based data-free design approach, where the diffractive features were optimized so that the all-optical intensity transformation of the diffractive processor achieves $\widehat{{\mathit{A}}_{\mathit{r}}}\approx {\mathit{A}}_{\mathit{r}}$. To evaluate $\widehat{{\mathit{A}}_{\mathit{r}}}$, we used $E{N}_i$ intensity vectors $\{{\mathit{i}}_{\mathit{r},t}{\}}_{t=1}^{E{N}_i}$, where ${\mathit{i}}_{\mathit{r},t}[l]=1$ if $l=t$ and 0 otherwise. In other words, $\{{\mathit{i}}_{\mathit{r},t}{\}}_{t=1}^{E{N}_i}$ are unit impulses, located at different input pixels. We simulated the all-optical output intensity vectors ${\{{\mathit{o}}_{\mathit{r},t}^{\prime }\}}_{t=1}^{E{N}_i}$ corresponding to these unit impulses and stacked them, i.e., $$\widehat{{\mathit{A}}_{\mathit{r}}^{\prime }}=[{\mathit{o}}_{\mathit{r},1}^{\prime }|{\mathit{o}}_{\mathit{r},2}^{\prime }|\cdots |{\mathit{o}}_{\mathit{r},E{N}_i}^{\prime }].$$(11)

Finally, we compensated for the optical diffraction efficiency-related scale mismatch through multiplication by a scalar, i.e., $$\widehat{{\mathit{A}}_{\mathit{r}}}=\sigma \widehat{{\mathit{A}}_{\mathit{r}}^{\prime }},$$(12)where $\sigma$ was defined as $$\sigma =\frac{{\sum }_{n=1}^{E{N}_i}{\sum }_{m=1}^{E{N}_o}{\mathit{A}}_{\mathit{r}}[m,n]\widehat{{\mathit{A}}_{\mathit{r}}^{\prime }}[m,n]}{{\sum }_{n=1}^{E{N}_i}{\sum }_{m=1}^{E{N}_o}{(\widehat{{\mathit{A}}_{\mathit{r}}^{\prime }}[m,n])}^2}.$$(13)

The MSE loss function to be minimized was defined as $${\mathcal{L}}_{\mathrm{PSF}}=\frac{1}{N_i{N}_o}\sum _{n=1}^{E{N}_i}\sum _{m=1}^{E{N}_o}{({\mathit{A}}_{\mathit{r}}[m,n]-\widehat{{\mathit{A}}_{\mathit{r}}}[m,n])}^2.$$(14)

The height $h$ of the diffractive features at each layer was constrained between zero and a maximum value $h_{\max }$ by employing a latent variable $h_{\text{latent}}$. The relationship between the constrained height $h$ and the latent variable $h_{\max }$ was defined as $h=\frac{h_{\max }}{2}\times [\sin (h_{\text{latent}})+1]$, where we chose $h_{\max }\approx \frac{\lambda }{n-1}$, which corresponds to a differential phase modulation of $2\pi$. The latent variables were initialized randomly from the standard normal distribution $\mathcal{N}(\mathrm{0,1})$.

The optimization of the diffractive layers was carried out using the AdamW optimizer39 for 12,000 iterations, with an initial learning rate of ${10}^{-3}$. The model state corresponding to the minimum of the MSEs evaluated after every 400 iterations was selected for the final evaluation. The ${\mathrm{D}}^2\mathrm{NN}$ models were implemented and trained using PyTorch (v1.12.1)40 with Compute Unified Device Architecture (CUDA) version 12.2. Training and testing were done on GeForce RTX 3090 graphics processing units (GPUs) in workstations with 256 GB of random-access memory (RAM) and Intel Core i9 central processing unit (CPU). The training time of the models varied with the size of the models. For example, the model used in Figs. 2(b) and 2(c) took around 1 h for 12,000 iterations. Inference for each input vector with $N_{\phi }=\mathrm{20,000}$ takes around 30 s.

To visualize the all-optical transformation error in Fig. 2, we used the error matrix ${\mathit{\epsilon }}_{\mathit{r}}={({\mathit{A}}_{\mathit{r}}-\widehat{{\mathit{A}}_{\mathit{r}}})}^2$; here $(·{)}^2$ denotes an element-wise operation. To evaluate the error $\mathit{\epsilon }$ at complex linear transformation, we applied demosaicking to the columns of $\widehat{{\mathit{A}}_{\mathit{r}}}$ to form the block matrix $[\widehat{{\mathit{A}}_0}|\cdots |\widehat{{\mathit{A}}_{\mathit{E}-1}}]\in {\mathbb{C}}^{{\mathit{N}}_{\mathit{o}}\times {\mathit{E}\mathit{N}}_{\mathit{i}}}$. Here, the subscript $k$ represents that $\widehat{{\mathit{A}}_{\mathit{k}}}$ is measured by applying the columns of $e_k\mathit{I}$ as input and stacking the corresponding demosaicked (complex-valued) output vectors. Accordingly, we have $$\mathit{\epsilon }={|[e_0\mathit{A}-\widehat{{\mathit{A}}_0}|\cdots |e_{E-1}\mathit{A}-\widehat{{\mathit{A}}_{\mathit{E}-1}}]|}^2.$$(15)

Here, $|{·|}^2$ represents an element-wise operation.

For the evaluation of the image encryption strength, we computed the entropy separately for the real and imaginary parts of a complex image as follows: $$H^{\mathrm{Re}/\mathrm{Im}}(\mathit{x})=-\sum _i{p}_i^{\mathrm{Re}/\mathrm{Im}}·\log (p_i^{\mathrm{Re}/\mathrm{Im}}),$$(16)where we calculated the distribution of either the real or imaginary part (denoted by the superscript) over the pixels of $x$; here $x$ denotes the complex image. $p_i$ is the probability (normalized histogram count) for a certain pixel value $i$.

For the histograms presented in Fig. S3(b) in the Supplementary Material, the data set is adapted from the Extended MNIST (EMNIST).41 For the creation of the input complex images, we randomly selected two distinct images from the EMNIST data set, using one as the real part and the other as the imaginary part of the complex image. To ensure compatibility with the input dimensionality, these images were bilinearly downsampled to a resolution of $4\times 4\text{pixels}$. We randomly formed a set of 1000 such complex images to compile the histograms presented in Fig. S3(b) in the Supplementary Material.

Biographies of the authors are not available.