Inverse-designed nanophotonic neural network accelerators for ultra-compact optical computing

Domain: chemistry
Published: 2026-03-04
Editorial outcome: minor_revision

Abstract
Inverse-designed nanophotonic devices offer promising solutions for analog optical computation, where high-density photonic integration is critical for scaling computational complexity. Here, we present an inverse-designed photonic neural network (PNN) accelerator, enabling ultra-compact and energy-efficient optical computing. Using a wave-based inverse-design method based on three-dimensional finite-difference time-domain simulations, we exploit the linearity of Maxwell’s equations to reconstruct arbitrary spatial fields through optical coherence. Each subwavelength voxel serves as a trainable degree of freedom, yielding a computational density of approximately 400 million parameters per mm². By decoupling the forward-pass process into linearly separable simulations, our approach is highly amenable to computational parallelism. We experimentally demonstrate two inverse-designed PNN accelerators, achieving on-chip MNIST and MedNIST classification accuracies of 89% and 90% respectively, with footprints of 20×20 µm² and 30×20 µm². Our results establish a scalable, energy-efficient platform for photonic computing, bridging inverse nanophotonic design with high-performance optical information processing.

Introduction
Machine learning has transcended conventional signal analysis techniques in numerous research fields ranging from computer vision to natural language processing1,2. Central to this progress is the ability of neural networks to learn complex data representations, enabling robust and adaptive inference. However, as neural models continue to grow in complexity, conventional electronic hardware faces increasing challenges in meeting the demands of computing speed and energy consumption3–6. In response, photonic neural networks (PNNs) have emerged as a promising platform for analog neural computation, offering key advantages such as ultrafast processing and low-power consumption4–10. PNN accelerator hardware based on programmable photonic integrated circuits such as Mach-Zehnder interferometer (MZI) meshes7, microrings8, and waveguide attenuators9, enables in situ training10 and dynamic weight adjustment through optical backpropagation10,11. By contrast, task-optimized architectures are trained offline12 and then realized in hardware for stable and efficiency-critical inference6. This paradigm, analogous to application-specific integrated circuits (ASICs) in electronics, eliminates static tuning power and runtime weight calibration13, leading to compact footprints with low-latency, stable inference14. Moreover, physically embedded weights enable in-memory optical computing4, significantly alleviating memory-bandwidth constraints present at run-time in conventional computing architectures at scale15.

To unlock the full potential of PNNs, it is essential to achieve high-density integration capable of supporting complex models and large datasets on-chip3,6,16. Our approach emphasizes in-memory optical single-shot run-time operation, with highly efficient, stable and low-latency inference. Traditional photonic design approaches, while effective for standard components, often fall short when optimizing for compactness, multifunctionality, and performance simultaneously.

In this context, topology optimization-based inverse-design methods have emerged as a powerful tool17–20. By computationally exploring a vast design space unconstrained by human intuition, it enables the discovery of non-intuitive geometries that maximize light-matter interaction within a minimal footprint21. Unlike diffractive models4–6,22, which rely on sequential, layer-by-layer architectures mimicking feedforward deep neural networks (DNNs), topology optimization allows arbitrary control over scattering, interference, and phase evolution throughout the device’s entire volume, enabling more compact and functionally dense implementations23–25.

However, the practical deployment of topology-optimized PNNs to accommodate increased dataset sizes remains challenging due to the absence of a scalable, physics-accurate training framework. Moreover, increasing the refractive index contrast enhances signal expressivity, which refers to the capacity for a rich set of field transformations within a fixed geometry21, through strengthening confinement, supporting internal resonances26, and enabling stronger interference and engineered scattering effects27,28 - all of which are critical for high-fidelity optical computation.

Consequently, three dimensional solvers such as three dimensional finite-difference time-domain (3D-FDTD) are required, since approximation methods such as the effective index method often fail to accurately capture the resulting out-of-plane scattering dynamics28–30. Moreover, as neural architectures scale to support increasingly complex tasks and larger datasets, computational efficiency and scalability become essential. Thus, achieving ultra-compact PNNs suitable for dense on-chip integration necessitates a design approach that captures the full-wave physics of high-index-contrast nanophotonic structures, while remaining computationally scalable and fabrication feasible.

Here, we demonstrate ultra-compact, inverse-designed PNN accelerators fabricated on a high-index-contrast silicon-on-insulator (SOI) platform operating at a single wavelength, where higher index contrast strengthens light confinement and supports rich interference and mode-mixing dynamics27,28, thereby expanding the space of attainable optical responses per unit volume. To enhance the inverse-design efficiency, we develop a wave-based inverse-design approach grounded in 3D-FDTD simulations, which exploits the linearity of Maxwell’s equations to reconstruct arbitrary spatial fields through optical coherence. By decoupling the forward-pass process into linearly independent simulations, our approach enables efficient parallelization and scales naturally across high-throughput computing hardware, such as graphics processing units (GPUs). Each subwavelength voxel acts as a trainable degree of freedom, yielding an ultra-compact representation of optical transformations with a computational density of approximately 400 million parameters per mm². The fabricated devices, with footprints as small as 20×20 µm² and 30×20 µm², perform on-chip image classification, demonstrating their efficacy on benchmarking datasets (MNIST31 and MedNIST32–34), and achieving experimental accuracies of 89% and 90%, respectively. Our inverse-designed PNN demonstrates the practical feasibility of ultra-compact nanophotonic devices for accelerating on-chip data processing. Owing to the method’s scalability and inherent compatibility with parallel computing architectures, this work is extensible to higher-dimensional processing tasks, offering a scalable pathway to efficiently handle increasingly demanding computational workloads.

Results
PNN accelerator configuration and inverse-design

Figure 1a presents a schematic overview of the on-chip PNN accelerator, wherein compressed and flattened input features from the dataset are encoded onto coherent optical amplitudes at a single wavelength (λ). Within the topology-optimized scattering region, the encoded optical fields undergo complex interference and scattering interactions throughout the device volume. The scattering region itself is designed to support non-uniform optical transport21,28, with routing behavior that adapts to the spatial distribution of optical inputs, efficiently redistributing optical power from the N input encoding ports toward the C class output ports. We denote C as the number of classes in the dataset and B as the output dimension of the PNN, where B=C in the one-hot setting. The PNN outputs represent class-evidence scores. A probability distribution across classification categories is obtained after photodetection, analogous to the output layer of a conventional multi-layer perceptron network (Fig. 1b).

To optimize the device material distribution across the scattering region, we directly perform the optimization within the binarized regime, while incorporating fabrication process-specific design constraints every epoch (see Methods). This approach removes the necessity for hyperparameter scheduling and post-processing discretization35. As illustrated in Fig. 1c, the latent design feature space, which varies continuously between 0 and 1, undergoes low-pass filtering to smooth parameter spatial distributions and is subsequently binarized36 into fabrication-compatible materials: silicon (Si) and silicon dioxide (SiO) as representative examples. A custom material distribution, characterized by its permittivity matrix ϵ(λ)=ϵr(λ)+iϵi(λ), where ϵr(λ) and ϵi(λ) are its real and imaginary components respectively, can model material loss and dispersive effects. We apply morphological operations to remove small features37, and B-spline contour approximations38–40 to further refine and shape the final geometry under minimal feature size and radius of curvature constraints, ensuring manufacturability41.

To accurately capture the performance of our design applied to multi-class image classification, it is suitable to minimize the cross-entropy (CE)1 loss between the one-hot encoded ground truth labels and the output optical power vector. This objective function indicates degraded performance when the optical power at the target output port is low. The corresponding output probability distribution across classes is computed by taking the surface integral of the optical power density at each port normalized by their total sum24. As a result of this normalization, the CE loss function implicitly penalizes higher distributed power at non-target ports, which is important for limiting channel crosstalk. The gradient of the loss function with respect to the full design space is computed using the adjoint variable method (AVM)42,43 and is proportional to the overlap of the forward mode fields and the adjoint fields. For each classification category, an adjoint source is constructed from the one-hot encoded target vector, which is used to derive the reverse-mode field (Eadj).

To simulate the forward-pass optical fields for a dataset containing L samples, each comprising N features, we exploit the linearity of Maxwell’s equations to significantly reduce computational complexity. Rather than performing L separate 3D-FDTD simulations, we reduce the simulation cost by executing one independent simulation for each of the N-input encoding ports, as illustrated in columns 1 and 2 of the inset in Fig. 1c. Leveraging the coherent nature of optical fields44, one can derive sample-wise fields Efwd‘ as a linear combination of pre-computed mode source fields E=[E1, E2, ..., EN], such that Efwd‘ = Σi=1N ai‘Ei, where ai‘ are the coefficients mapped to dataset features of sample ‘ and Ei is the full 3D field from the i-th mode source. More generally, the entire dataset of reconstructed optical fields can be compactly expressed as Efwd = DTE, where D ∈ RN×L is the dataset feature matrix encoding all input samples, and Efwd is the tensor containing the corresponding forward fields. The information is reshaped into a complex tensor with dimensions [L, nx, ny, nz], where nx, ny and nz denote the matrix dimensions in the x, y and z-axes respectively. From this tensor, a subset T⊂L is selected as training fields, which directly inform the iterative updates during inverse-design optimization. Representative training fields corresponding to sample images from the MedNIST dataset are shown in Fig. 1c, inset column 3. The exact gradient update, as illustrated in Fig. 1c inset, column 4 is obtained by solving the AVM equation across the training samples. This gradient acquisition differs from straight-through estimators, which approximate the physical gradient through binary activation functions using an identity operator - a process which is a commonly used in training quantized DNNs45. Using stochastic gradient descent with momentum (Fig. 1c), we iteratively update the latent parameter distribution to minimize the CE loss, thus progressively increasing the classification accuracy during optimization (see Methods).

The total number of required 3D-FDTD simulations per epoch is only N+C, composing N simulations (forward mode) to precompute input-related fields, plus C simulations (reverse mode) to compute gradients for each class category. This approach substantially lowers computational costs compared to directly simulating the entire dataset, especially beneficial when dealing with large datasets, where the total dataset size L is significantly greater than N+C2,46. Furthermore, by virtue of the linear separability of the forward mode spatial fields, this inverse-design approach is highly amenable to computational parallelism, making it particularly well suited for acceleration on GPU and other parallel computing architectures, thereby enhancing scalability across large problem domains.

On-chip image classification tasks

To demonstrate the scalability of our inverse-design scheme applied to image classification, we perform numerical experiments using two benchmark datasets with progressively increased complexity: the MNIST dataset, containing 70,000 grayscale digit images31 with 10 input features (N=10), 10 classes (C=10) and a design region footprint 20×20 µm², and the MedNIST containing 58,954 medical images32–34 with 15 features (N=15) across six classes (C=6) and a design region footprint 30×20 µm². For both datasets, the input data features are spatially encoded onto the amplitude of light at a wavelength of 1550 nm. The MNIST dataset is partitioned into 50,000 training images and 20,000 test images, while the MedNIST dataset includes 48,000 training images and 10,954 test images. Given that these datasets are in grayscale, it is only necessary to restrict the problem domain to the time-harmonic regime in numerical simulations21. The MNIST classification results are shown in Fig. 2. Figure 2a illustrates the normalized power distribution across all output ports and the corresponding electric field (E-field) propagation at epochs 0, 60, and 120 for a randomly selected example from Class ‘5’.

Initially, the randomized high-contrast material topology does not yield meaningful classification outcomes. However, as topology optimization progresses, optical energy becomes increasingly concentrated at the correct output port corresponding to Class ‘5’, while power at neighboring ports is effectively suppressed. The corresponding continuous wave (CW) spatial field response of the design at the final epoch is presented in Supplementary Movie 1, showing the evolution of the optical field distribution over multiple optical cycles.

Figure 2b depicts test loss and accuracy as a function of training epochs, with range intervals capturing variability over three separate training sequences where the datasets are randomly permutated. The final confusion matrix (Fig. 2c) demonstrates a high overall classification accuracy of 97.8%. In addition, the energy density matrix (Fig. 2d), as a measure of the output power distribution averaged across the test set, confirms effective energy localization, showing the optical energy is well concentrated at the correct output ports corresponding to each class, thereby indicating minimal cross-class optical crosstalk. Moreover, Fig. 2e illustrates the classification accuracy with respect to training data size, reinforcing the need for high volumes of training data to perform accurate signal classification. Notably, despite the increase in training data volume, our method maintains a constant number (N+C=20) of required 3D-FDTD simulations per epoch in the inverse-design. As a result, the increase in epoch runtime is just 6.7% when using the training full dataset compared to only 10% of it.

Next, we perform classification on the MedNIST dataset to evaluate the scalability of the proposed method to a higher-dimensional dataset with a larger volume of data. Figure 3a shows the inference fields for a single randomly selected sample from each of the six MedNIST classes, along with their corresponding normalized power distributions, demonstrating correct signal localization at the respective class output port. The CW spatial field response for the ‘BreastMRI’ sample at a wavelength of 1550 nm is presented in Supplementary Movie 2, illustrating how the optimized PNN structure effectively concentrates optical power at the correct output port. As shown in Fig. 3b, the optimizer reached an equilibrium after 150 epochs, yielding a peak numerical accuracy of 99.1%. The optimized material topology at the final training epoch is shown in the 3D-rendered schematic (Fig. 3b Inset), highlighting the 220 nm fully-etched Si waveguide designed to enable efficient wave propagation through its high refractive-index contrast. The confusion matrix (Fig. 3c) quantitatively confirms this high classification accuracy, reaching 100% accuracies for classes ‘ChestCT’, ‘AbdomenCT’, ‘HeadCT’, and ‘BreastMRI’. Furthermore, the energy density matrix (Fig. 3d) demonstrates effective energy localization, with approximately half of the total output power concentrated at the correct output port for five classes. Although slight coupling to neighboring ports is observed, the device maintains high performance even within its ultra-compact footprint.

Figure 3e further validates the performance of our 3D-based optimization approach when the volume of training data increases. The observed classification accuracy improves consistently with increasing dataset size, while the computational overhead remains minimal due to the fixed number of FDTD simulations per epoch. The MNIST and MedNIST PNNs comprise approximately 1.6×10⁵ and 2.4×10⁵ trainable parameters, respectively (see Supplementary Note 1). During each training epoch, the corresponding multiplication operations are parallelized on GPUs, enabling fast and scalable gradient updates. The computational density reaches approximately 400 million trainable parameters per mm² (Supplementary Note 1), highlighting our PNNs’ capability to achieve high computational density within an ultra-compact footprint.

Hardware performance and inverse-design computational scalability

Figure 4 shows the benchmark results for a single simulation job, including key FDTD simulation parameters: mesh time, FDTD time and total wall time for MNIST and MedNIST datasets. The inverse-design process is conducted across three GPU platforms: Nvidia RTX 5090, RTX 4090, and V100. It is observed that the FDTD time for V100 is approximately 64 seconds to complete a single simulation for MNIST, compared to 48 seconds for RTX 4090 and 30 seconds for RTX 5090. The wall-clock time for a single simulation also accounts for factors such as mesh time - typically CPU-bound. Moreover, given the size of the MedNIST chip is 1.5 times larger than MNIST, this increased footprint is reflected in the wall-clock time, which is 1.46 times longer than MNIST.

The cumulative wall-clock time for the inverse-designed PNN is defined as the product of the single wall-clock time, the fixed number (N+C) 3D-FDTD simulations per epoch, and the total number of epochs. The total simulation time accounts for both the cumulative wall-clock time and the total time required for each epoch to derive spatial fields across the entire dataset. Our approach enables the linear separability of each simulation, allowing individual simulations to be independently executed on a clustered computing network47. This capability supports scheduling via makespan minimization, a strategy that optimally balances computational loads across available resources48, significantly reducing overall compute time through parallelization. The inset of Fig. 4 illustrates the total simulation time under three distributed computing schemes. On a single GPU Linux node equipped with an Intel Core i7-processor, 128 GB of RAM, and an Nvidia RTX 5090 (32 GB) graphics card, the inverse-design optimization for MNIST PNN runs for approximately 29.7 hours. Using the same hardware configuration, the MedNIST PNN optimization requires approximately 56.3 hours. Combining RTX 4090 and RTX 5090 nodes reduces these times to 19.7 hours (MNIST) and 37.9 hours (MedNIST). Incorporating an additional V100 node further reduces the total simulation time to 17.1 hours for MNIST and 33.3 hours for MedNIST. These benchmarking results clearly demonstrate the effectiveness of our approach in leveraging parallel computing.

Fabrication and measurement

To experimentally demonstrate the feasibility of our designs, we fabricate PNN accelerators for MNIST and MedNIST datasets on the SOI platform (see Methods). Figures 5a and b show the scanning electron microscope (SEM) images of the fabricated 20×20 µm² MNIST and 30×20 µm² MedNIST PNN accelerators, respectively, highlighting their compact footprints. Both input and output waveguides are designed to support the fundamental transverse-electric (TE) mode at 1550 nm, with a width of 500 nm and height of 220 nm. The minimal feature size across the device is 80 nm.

Optical signals from a CW laser source at 1550 nm are vertically coupled into the chip via vertical grating couplers (VGCs), distributed through a power splitter, and encoded with amplitude and phase information for input data preparation via thermo-optic tuning (see Methods), and then processed via the inverse-designed PNN. To illustrate full device integration, Fig. 5c shows the MedNIST PNN chip wire-bonded to a printed circuit board (PCB) (see Methods) as the electrical interface through which optical amplitude and phase are adjusted for input data preparation. A detailed microscope image of the MedNIST chip is shown in Fig. 5d, clearly illustrating the VGC (Inset 1) for light coupling in and out of the chip, the multimode-interferometer (MMI) power splitter (Inset 2), MZIs (Inset 3), and phase-combiners (Inset 4) for input data preparation, along with the compact PNN accelerator region (Inset 5) and monitoring ports. These monitoring ports enable direct measurement of optical power at critical locations within the device, facilitating real-time calibration (see Supplementary Fig. 1). Finally, the optical power distribution from six output ports of the PNN accelerator is simultaneously measured using a multiport optical power meter to obtain classification results (see Supplementary Fig. 2). The MNIST PNN chip follows a similar layout and packaging (see Supplementary Fig. 3).

For experimental data preparation, we measure analog optical input power to the PNN accelerators across 100 MNIST and 60 MedNIST test samples, where 10 samples are randomly selected from each class (see Supplementary Fig. 4). To verify the accuracy of these inputs, we measure a global mean average error (MAE) of 0.0277 for MNIST and 0.0310 for MedNIST with respect to the simulated inputs (see Supplementary Fig. 5). These low MAE values confirm a strong correlation with the simulated amplitude-encoded features, validating the accuracy of our optical data preparation process. To evaluate how well the experimental dataset represents the test distribution, we compute several distributional benchmarks. A 3D principal component analysis shows that the experimental and test set features largely overlap within their 90% confidence ellipses (see Supplementary Fig. 6a, b). The mean feature distributions yield low per-feature MAE values of 0.0273 for MNIST and 0.0336 for MedNIST, and small average 1D Wasserstein distances49 of 0.0212 and 0.0279, respectively, indicating high distributional agreement (Supplementary Fig. 6c, d). Bootstrap resampling50 over 10,000 iterations (Supplementary Fig. 7a, b) produces global root mean square error precisions of ±2.02% for MNIST and ±3.63% for MedNIST relative to the test set mean at 90% confidence (Supplementary Fig. 7c, d). These results confirm that the stratified subset adequately represents the overall test distribution and that our experimental data remain statistically robust under resampling.

For MNIST, an overall classification accuracy of 89% is achieved, with Class ‘1’ exhibits the highest classification accuracy of 100%. For the MedNIST dataset, we achieve an overall classification accuracy of 90%, with the ‘BreastMRI’ and ‘CXR’ classes reaching the highest classification accuracy of 100%. Moreover, the MedNIST energy density matrix reveals an average of 50.5% of total output power concentrated at the correct output port. Among the correctly classified output ports, our results indicate only a 16.1% variation in maximum output power, demonstrating consistent classification performance across a higher-dimensional sequence of data. The classification score can be further improved by minimizing the power nonuniformities among the VGCs through enhanced fiber-to-chip alignment tolerance51. The deviation between simulated and experimental results mainly arises from fabrication-induced geometry errors. When the fabrication deviation (see Supplementary Fig. 8a) caused by lithographic over- or under-exposure41,52 reaches approximately ±20 nm, changes in classification accuracy become evident for both MNIST and MedNIST devices (see Supplementary Fig. 8b, c). The tolerance to fabrication deviations can be enhanced by independently simulating geometric perturbations53 and updating designs informed by these results, while training with augmented datasets6 can further improve device robustness.

Robustness verification

To evaluate the robustness of the PNN devices experimentally, we explore the effect of input data phase deviations on classification performance (see Methods). Figures 6e and f show scatterplots of the maximum phase variation while still maintaining a correct classification for MNIST and MedNIST datasets, respectively. Despite large input phase swings, the classification accuracy remains robust and consistent across both datasets. This robustness is primarily attributed to the amplitude-dominated encoding scheme, where classification decisions are governed predominantly by the strongest optical intensity features, reducing sensitivity to phase fluctuations. Figure 6g shows representative examples from the MNIST and MedNIST datasets under different mean absolute phase deviation conditions of 0.89 and 1.18 radians, respectively, clearly illustrating the dominant energy distribution at the correct output ports despite these phase deviations.

PNN model scalability to larger neural networks

We scale the model to larger neural networks by stacking PNNs with interleaved photodetector nonlinearity54, where each PNN is inverse-designed for task-optimized optical-field transformations. Figure 7 illustrates the schematic architecture for the PNN model scalable to larger neural networks. Input images are segmented into patches, encoded as feature vectors, and processed by pairs of consecutive PNN cores. The first PNN maps input features to an optical-field distribution across class channels (Fig. 7 Inset a), which are converted after photodetection into a series of patch class-activation vectors55. The second PNN re-encodes this vector into refined patch features (Fig. 7 Inset b), which become increasingly aligned with class prototypes as they propagate in the next stage. Each pair is analogous to a layer block in multilayer DNNs56, where blocks are stacked to extract and refine feature representations57. Successive PNN blocks progressively refine these class-level embeddings, filtering out patches with low confidence e.g., background, and sharpening those with high-class correlation, generating richer interpretability of feature-to-class associations58. This depth-wise refinement expands the effective receptive field59, while a sliding window with a fixed stride integrates broader spatial context. At the network output, patch-level predictions are confidence-weighted via monitored optical power and aggregated to yield the final image-level classification result (Fig. 7 Inset c). Different output encoding schemes, such as one-hot, binary, or error-correcting output codes (ECOC)60 are viable for classification, where multi-class decisions are recovered by code-word decoding at inference54. Weight-sharing across network width enables patches to be evaluated in parallel61. Moreover, optical multiplexing62 such as wavelength or polarization, increases inference throughput and enables multi-task operation to be executed on the same chip (see Supplementary Note 2). By using dual wavelengths, we demonstrate a 30×20 µm² single-chip design which simultaneously classifies MNIST and MedNIST datasets, achieving 95.1% and 98.0% test accuracy, respectively (see Supplementary Note 2).

A patch-efficient adjoint method (see Supplementary Note 3) is incorporated into the inverse-design of each PNN core, ensuring that the training cost is independent of the number of image patches. For a PNN core with input dimension N and output dimension B, the number of FDTD solves per epoch is fixed at N+B under single-wavelength operation (see Methods). This scaling keeps training tractable even for large-class datasets. For instance, with a 4×4 patch kernel (N=16) on ImageNet (1000 classes), where B=16 is determined by ECOC coding54 with a 10-bit minimum code length plus 6 redundancy bits for error correction (see Methods), each PNN core design update requires only 32 FDTD-solves.

Methods
Fabrication-informed inverse-design parametrization

At each inverse-design optimization step, we consider the parametrization of a latent feature space into a fabrication-feasible, binarized material index distribution. These latent parameters vary continuously between 0 and 1, which correlates to SiO2 (ϵmin) and Si (ϵmax), respectively for SOI waveguide. The feature space is randomly sampled from a uniform distribution then low-pass filtered using a Gaussian blur kernel71. To embed fabrication constraints into the design, we first employ a low-pass filter via a conic kernel41. This filtered latent space is projected onto a binary distribution using the modified hyperbolic tan function36 and is fully-binarized to enforce a real material permittivity matrix. The material is biased more toward Si (η < 0.5), where η is the material biasing parameter. A B-spline approximation is utilized to further enforce minimum feature size and radial constraints in the graphic data system (GDS) design file, while facilitating a spatially adaptive resolution and anisotropic topology-subject to considered placement of control points38. A discrete contour which defines the boundary between material and void regions is first constructed using the Suzuki-Abe algorithm72, followed by small feature removal using the skimage-morphology package in Python. The B-spline curve itself is composed of a linear combination of basis functions in 2D, expressed using the Cox-de Boor recursion formula39,40. By tailoring both the degree of the curve k and the number of knots, the density of approximated curves is effectively relaxed, thereby informing minimum feature and curvature constraints. In our case, a cubic (k=3) B-spline approximation maps the discrete material voxels to a topology with a minimum feature size of 80 nm. Open-source Python library gdspy is used to translate these control points to create a GDS file at each design epoch.

Numerical simulation setup

To simulate the optical field distribution, a full-wave 3D-FDTD method (ANSYS Lumerical) is used. A uniform mesh is chosen with a step size less than λ/30 to be equivalent in dimension to the parameter space, and a non-uniform mesh is used outside the design region. A lightweight, semi-supervised autoencoder is used to learn a low-dimensional latent representation of the input data, which is then normalized and encoded onto the amplitude of light at 1550 nm. Open-source Python framework cupy is used to compute and store the 128-bit complex fields, each of size [n, n, n]) in the GPU’s temporary memory. For both designs, a mini-batch size of 250 training samples is chosen to smooth the gradient distribution and mitigate oscillatory behavior in optimization updates. The latent weights are randomly initialized, and Adam73 is used as the stochastic gradient optimizer, with a learning rate of 0.005 and exponential decay rates of 0.667 and 0.9 for the first and second-moment estimates, respectively.

Device fabrication and prototyping

Both the MNIST and MedNIST PNN accelerators are fabricated via electron beam lithography with a foundry-compatible minimum feature size (80 nm)74 on a standard SOI wafer, which has a 220 nm thick top silicon layer sitting on top of a 2 µm thick buried oxide layer, above a 725 µm-thick silicon substrate. VGCs43,75 are fabricated with an etch depth of 70 nm. Plasma-enhanced chemical vapor deposition process is used to deposit a 1 µm layer of SiO2 onto the fabricated devices to act as an insulation layer between the electrodes and the device, which is sufficiently thick to minimize excess absorption due to the metal electrodes76. To facilitate thermo-optic tuning, thin-film titanium microheaters and gold traces are fabricated atop the SiO2 cladding layer via a maskless lithography system (Heidelberg MLA100). For packaging and electrical interfacing, the devices are mounted onto a custom-designed PCB, which utilizes electroless nickel electroless palladium immersion gold contacts to provide the wire bonding material interface and electrical connectivity to a central MCU. The photonic chip is affixed to a thermoelectric cooler to ensure thermal stabilization. Temperature regulation is achieved via a thermistor embedded near the chip and connected through plated vias to the PCB.

Data preparation, calibration, and phase robustness testing

Input light at 1550 nm from a tunable laser (Keysight N7778C) with a relative intensity noise less than −150 dB/Hz77 is fiber-coupled into the chip. A set of 10 and 15 MZIs is used for MNIST and MedNIST PNN devices, respectively to facilitate amplitude control as data preparation. Each MZI is composed of a pair of 50:50 MMI splitters, integrated tapers, and S-bends with 30 µm arc radius. Downstream of MZIs, asymmetric 10:90 MMI are employed to tap a small portion of optical power for monitoring the attenuation in each channel, measured using multiport optical power meters (Keysight N7744C) with dynamic range of −65 dB (exceeding 80 dB with range switching)78 that incorporate built-in averaging to suppress random thermal and shot noises at the photodiode readout79. To monitor optical phase, adjacent channels are similarly tapped using 10:90 MMIs, and the extracted signals are recombined through 50:50 MMI combiners. The phase difference between channels is inferred from the combined optical power, where maximized combined optical power corresponds to phase equalization. For both MNIST and MedNIST chips, a stratified random sample of 10 examples per class is taken from the simulated test set to form the experimental validation set, demonstrating diverse feature representation (see Supplementary Fig. 6). The corresponding input channels are attenuated using MZI-based amplitude control, with feedback from the optical power meter, while simultaneously monitoring and equalizing the relative phase between channels. For both amplitude and phase calibration, the voltage generated by a digital-to-analog converter is swept between 1.1 V and 5 V with increments of 3 mV. The optimal voltage is identified by minimizing the error between the ground-truth and experimental optical transmission at 1550 nm. For phase robustness testing, the phase monitoring ports are used to identify the optimal phase corresponding to the highest classification result - this phase is then served as a reference. The applied phase variation relative to this reference is inferred from the electrical power delivered to the thermo-optic phase shifters, based on their known tuning characteristics80.

Training framework and scalability

Each PNN core is trained to implement task-specific optical-field transformations by minimizing a CE loss1. Class information is embedded through one-hot, binary, or ECOC coding60, with ECOC advantageous for large-class datasets by reducing the number of output channels required54. A patch-efficient adjoint method (see Supplementary Note 3) is incorporated into the optimization, whereby patch contributions are aggregated algebraically before gradient computation. This ensures that the total number of required FDTD solves per epoch for one per PNN core depends only on the input and output dimensions (N+B) and is independent of the number of image patches (see Supplementary Note 3). The optimization incorporates nonlinear photodetection, which converts optical logits into class-activation vectors and provides interpretable measures of patch-to-class association55. To stabilize training and ensure global consistency, auxiliary supervision can be introduced by aggregating patch-level predictions with confidence weighting at intermediate layers, enforcing agreement with image-level classification. Training is proceeded block by block, where each consecutive pair of PNN cores is optimized jointly and then frozen after convergence. Overlapping patches are re-extracted and passed to the next stage, enabling depth scalability through successive refinement. Reuse of pre-trained shallow blocks as priors can further improve efficiency by reducing the number of gradient updates needed for deeper layers81. The final PNN core maps refined patch features into class-activations, and confidence-weighted pooling across patches produces the image-level prediction. All optimizations are performed in simulation, and the resulting geometries are used directly for inference in the fabricated device.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present an interesting paradigm shift, where machine learning-based design does not aim at generating typical\nphotonic components (filters, splitters etc.), but on the contrary, the task at hand is the generation of a scattering-based\nsilicon-on-insulator chip able to do image classification, with optimized performance. The manuscript contains in-depth\nnumerical simulations alongside a small validation using two SOI chips. The tasks used so as to demonstrate the proof of\nconcept are the well-established MNIST and MedNIST, offering relatively high classification performance, with minimum\ntraining-designing downtime.\nOverall, the proposed idea is interesting, the manuscript is well written, well organized,\ncontains both numerical and experimental solid results, and contains the necessary references so as to support its claim. On the other hand, there are some issues that need to be addressed prior to publication.\n• The most critical question that arises concerns the concept itself. In particular, the authors propose a photonic integrated\nchip that should be designed using machine learning to solve a specific machine learning task. In standard machine\nlearning, pre-trained models can tackle multiple problems without adapting their structure or even their training. In the\nintegrated photonic domain, multiple configurable structures (based on MZI meshes, reconfigurable PICs, filters etc.) have\nbeen proposed that are capable of adjusting (weights or even the architecture dynamically) to tackle multiple tasks. On the\nother hand, the proposed scheme suggests a static implementation that is designed for a single task and cannot be\nreconfigured to adapt to new tasks. It would be beneficial for the manuscript if the authors included in their introductory part a\nthorough comparison of the potential benefits of their system compared to other adaptable PIC implementations.\n• Similar to the above, the authors simulate/design and fabricate two different devices with different dimensions and different\nstructure to tackle a similar classification task (MNIST and MedNIST). It would demonstrate the potential flexibility of the\napproach if the larger (MedMNIST) was used in order to classify the MNIST task or vice versa. Such a small-scale\nexperimental or numerical investigation would prove the adaptability of this approach; otherwise, I am not confident about its\nuse.\n• In lines 38–39, the authors claim that dense and large-scale chips should be developed in order to address complex\nproblems. On the other hand, in multiple cases, PNNs exploit time multiplexing and increased bandwidth capabilities so as\nto fit larger models on small-scale chips. In this context, it would enhance the quality of the manuscript if the authors included\nsimilar works and highlighted the merits of their approach compared to such schemes.\n• Furthermore, the authors do not mention or analyze the bandwidth of the amplitude-modulated signals used in\nexperimental validation. This is a critical issue that governs the throughput of such PNNs. In addition, the authors should\ncomment on whether modulation speed affects performance and system training.\n• In line 51, the authors mention that high contrast materials and tight confinement increase signal expressivity. The\nstatement is not clear—what do the authors define as signal expressivity, and how exactly is it affected by light confinement?\n• Another fundamental question is that in most photonic-based computing, inference comes at low latency, marginal power\nconsumption, and increased throughput, while training in most cases entails elaborate simulations using a digital twin\napproach. In these cases (e.g. MZI meshes), the models used can be downgraded to simple transfer function combinations\n(scattering matrix). On the other hand, in this work, inference is equally straightforward, but training-designing entails the use\nof costly physically accurate models. It would be beneficial if the authors could comment on this and provide a table\nincluding the number of parameters (multiplications) needed to train (offline) different PNNs compared to their approach.\n• In line 96, the authors mention that their design allows spatially adaptive transport. What do they mean by this? How is\nadaptability achieved in such a structure?\n• In line 98, the authors state that the output of a Perceptron using PNN is a probability function. That is not correct, as in\norder to transform optical outputs to probabilities, photodetection and digital post-processing (e.g. softmax) should be\napplied.\n• In lines 120–123, the authors describe briefly their strategy for minimizing computational resources. It is my belief that this\nsection is not clear and needs improvement. If I am not mistaken, each input, instead of being injected at each port for the\nsimulation, is injected only to one. If this is correct, does it impact the final accuracy? If the authors used a full-scale training\napproach, would performance be affected, and if yes, how much?\n• The authors generate two different chips to tackle two different tasks; the chips have different areas and offer different\nperformance. Should the size of the chip be scaled with the task’s complexity? And if yes, what type of chip would be\nneeded for an even more complicated task (e.g. ImageNet)? It is my opinion that the authors should provide at least\nnumerical results where they tackle MNIST with chips of different area/complexity so as to highlight this dependence.\n• Accuracy for the MNIST dataset is in the order of 97.8% when simulations are considered. There are effects that are not\ndiscussed within the manuscript regarding these simulations, e.g. the inclusion of thermal/shot noise at the photodiodes,\nlimited bit-resolution at the ADCs, nonlinearity or imperfections at the amplitude modulator, bandwidth limitations of\nequipment etc. All these effects are probably absent in the numerical simulations, thus the really high accuracy of 97.8%.\nOn the other hand, during experimental validation, even when using a small number of samples, accuracy radically drops to\n89%. This last value is not really impressive, based on the fact that 96%–97% accuracy has been demonstrated with PNNs\nin the past. It would be beneficial if the authors could comment on the sensitivity of their design to actual experimental\nconditions (explain this performance degradation) and compare their approach to state-of-the-art PNNs targeting MNIST.\n• In the same context, one key problem in PNN training using digital twins (offline training) is the fact that minor deviations\nbetween the physical chip and the model are translated to significant errors during training (designing here). In this case,\nwhere even nanometer location/size errors would significantly impact the scattering properties of the chip, imperfections\nbetween the physical system and its digital twin can result in major differences between expected performance (97.8%) and\nactual performance (89%). The authors should comment on this in their revised manuscript.\n• During experimental validation, the authors use an extremely limited number of samples from the dataset (100 out of\n60,000). What is the reason for this choice, and how can this affect the statistical properties of the solution?", "reviewer_2": "The results in this study is characterized by both physical modeling (characterization) and computational strategies. The\nauthors employ three-dimensional finite-difference time-domain (3D-FDTD) simulations as the backbone of their design\napproach. The methodology can reconstruct arbitrary optical field patterns through interference, providing a foundation for\nthe photonic neural network’s operation. Results capture the complex wave dynamics within the nanophotonic circuits rather\nthan relying on simplified or purely theoretical models.\nA particularly innovative aspect is how the linearity of optical propagation is used to improve efficiency and scalability.\nInstead of simulating the entire neural network in one monolithic 3D-FDTD run, the authors decouple the forward-pass\ncomputations into multiple linearly separable simulations. In practice, this means different portions or layers of the photonic\nnetwork are simulated independently and can be run in parallel, taking advantage of GPU computing. Running many small\nFDTD simulations at the same time speeds up the design process and makes it easier to scale to larger systems. This\nadvanced method uses the physics of light and powerful computing to design very compact optical devices.\nThe validity is reinforced by experimental verification. After designing the nanophotonic neural network in simulation, two\nprototype photonic chips on a silicon-on-insulator platform were fabricated and tested on standard image-classification\ntasks. The correspondence between simulation and experiment is good: the devices (~20×20 μm² in size) achieved high\nclassification accuracies, approximately 89% on MNIST and 90% on MedNIST. This fabrication and testing step validate the\nmodeling assumptions; the optical network behaved in practice as designed. The strength of this study is to build computer\nsimulation and physical implementation together. It provides a scalable and energy-efficient photonic computing platform in\ntheory and confirms it in experiment.\nDespite these strengths, several limitations and challenges should be mentioned.\nOne concern is the heavy computational demand of 3D-FDTD simulations. Achieving such a detailed level of modeling\naccuracy requires considerable computing resources and time, even with parallelization. The authors reduce this through\nGPU acceleration and by splitting the problem into smaller parts but scaling to truly large or deep neural networks may still\nface practical difficulties. As network size increases, the number of separate simulations and the volume of data could\nincrease dramatically, potentially offsetting the gains from parallelism. While the approach is mentioned scalable to large\nproblem domains, its efficiency for much more complex networks remains to be demonstrated.\nTwo, it was unclear on how these models can be scaled to larger neural networks with several layers. How will the training\noccur? Again with simulation? It was not clear on how to apply the inverse design to large CNN and DNN models.\nThird, is the need for careful incorporation of real-world constraints in the simulation. Factors like fabrication tolerances and\nmaterial losses should be accounted in design to ensure the simulated performance translates into actual devices. If the\ncurrent simulations assume ideal conditions, that could be a limitation. Slight manufacturing imperfections, for example, will\nneed to be solved for larger-scale or commercial applications. It appears to be an ideal case rather than any fabrication\nimperfections.\nIn summary, while the approach appears well suited for small MNIST datasets, it is unclear on the actual scalability of the\nproposed design. Device imperfections are not considered in the design. There is significant computational requirements\nand it was unclear on how the approach will scale with model sizes."}, "author_response": "We thank the reviewers for their comments and recommendations. We have carefully\nconsidered all feedback and revised the manuscript accordingly. Our detailed responses to each\nreviewer are provided below. The reviewer’s comments are shown in italics, and changes in\nthe manuscript are indicated in red.\nReviewer #1\nComments: The authors present an interesting paradigm shift, where machine learning-based\ndesign does not aim at generating typical photonic components (filters, splitters etc.), but on\nthe contrary, the task at hand is the generation of a scattering-based silicon-on-insulator chip\nable to do image classification, with optimized performance. The manuscript contains in-depth\nnumerical simulations alongside a small validation using two SOI chips. The tasks used so as\nto demonstrate the proof of concept are the well-established MNIST and MedNIST, offering\nrelatively high classification performance, with minimum training-designing downtime.\nOverall, the proposed idea is interesting, the manuscript is well written, well organized,\ncontains both numerical and experimental solid results, and contains the necessary references\nso as to support its claim. On the other hand, there are some issues that need to be addressed\nprior to publication.\nReply:\nThank you for reviewing our manuscript and for your positive and encouraging feedback.\nBelow, we carefully address each of your detailed questions, providing clarifications,\nadditional analyses and new results, as well as expanded discussions in the revised manuscript.\nComments: The most critical question that arises concerns the concept itself. In particular, the\nauthors propose a photonic integrated chip that should be designed using machine learning to\nsolve a specific machine learning task. In standard machine learning, pre-trained models can\ntackle multiple problems without adapting their structure or even their training. In the\nintegrated photonic domain, multiple configurable structures (based on MZI meshes,\nreconfigurable PICs, filters etc.) have been proposed that are capable of adjusting (weights or\neven the architecture dynamically) to tackle multiple tasks. On the other hand, the proposed\nscheme suggests a static implementation that is designed for a single task and cannot be\nreconfigured to adapt to new tasks. It would be beneficial for the manuscript if the authors\nincluded in their introductory part a thorough comparison of the potential benefits of their\nsystem compared to other adaptable PIC implementations.\nReply:\nWe thank the reviewer for their comment. In the integrated photonic domain, fully\nreconfigurable PNN architectures based on MZI meshes7, microrings8 and waveguide\nattenuators9 enable dynamic adjustment of weights and structure to support broad task coverage.\nThese approaches emphasize generality and online adaptability. While highly versatile, they\ntypically involve larger footprints and require complex electronic calibration11, which can limit\ninference throughput and stability13,14. Task-optimized architectures, by contrast, are trained\noffline and then realized in hardware for efficiency-critical inference6. This paradigm is\nanalogous to application-specific integrated circuit (ASIC) accelerators in electronics, where\nefficiency for a given task is prioritized over full flexibility. Our approach belongs to this\ncategory and has advantages of energy efficiency, ultra-compactness and single-shot stable\ninference. Moreover, given the weights are embedded physically post-optimization, inference\nproceeds as in-memory optical computing4 without continual weight fetches as characteristic\nof conventional computing architectures, which imposes a memory-bandwidth bottleneck at\nscale15.\nWhile a fabricated PNN core has fixed weights, the overall technique is not restricted to a static\nimplementation that is designed for a single task. Our design remains scalable through cascaded\ncomposition and multiplexing strategies, providing extensible functionality and multi-task\ncapabilities. To demonstrate this, we show that the same 30 × 20 μm² chip can perform two\ndistinct classification tasks in simulation. Wavelength multiplexing further allows task\nselection via input wavelength, while stacking PNN cores supports partial reconfigurability\nthrough evidence extraction, feature embedding, and pooling adjustments. Thus, although each\ncore is weight-fixed, the overall architecture remains scalable and extensible.\nWe have revised the manuscript to clarify these distinctions and justify the scalability and\nflexibility of our design by adding a new section “PNN model scalability to larger neural\nnetworks”.\nOn Page 2:\nIn response, photonic neural networks (PNNs) have emerged as a promising platform for\nanalog neural computation, offering key advantages such as ultrafast processing and low-power\nconsumption4-10. PNN accelerator hardware based on programmable photonic integrated\ncircuits such as Mach-Zehnder interferometer (MZI) meshes7, microrings8, and waveguide\nattenuators9, enables in situ training10 and dynamic weight adjustment through optical\nbackpropagation10,11. By contrast, task-optimized architectures are trained offline12 and then\nrealized in hardware for stable and efficiency-critical inference6. This paradigm, analogous to\napplication-specific integrated circuits (ASICs) in electronics, eliminates static tuning power\nand runtime weight calibration13, leading to compact footprints with low-latency, stable\ninference14. Moreover, physically embedded weights enable in-memory optical computing4,\nsignificantly alleviating memory-bandwidth constraints present at run-time in conventional\ncomputing architectures at scale15.\nOn Pages 16-18\n:\nPNN model scalability to larger neural networks. We scale the model to larger neural\nnetworks by stacking PNNs with interleaved photodetector nonlinearity54, where each PNN is\ninverse-designed for task-optimized optical-field transformations. Figure 7 illustrates the\nschematic architecture for the PNN model scalable to larger neural networks. Input images are\nsegmented into patches, encoded as feature vectors, and processed by pairs of consecutive PNN\ncores. The first PNN maps input features to an optical-field distribution across class channels\n(Fig.7 Inset a), which are converted after photodetection into a series of patch class-activation\nvectors55. The second PNN re-encodes this vector into refined patch features (Fig.7 Inset b),\nwhich become increasingly aligned with class prototypes as they propagate in the next stage.\nEach pair is analogous to a layer block in multilayer DNNs56, where blocks are stacked to\nextract and refine feature representations57. Successive PNN blocks progressively refine these\nclass-level embeddings, filtering out patches with low confidence e.g. background, and\nsharpening those with high-class correlation, generating richer interpretability of feature-to-\nclass associations58. This depth-wise refinement expands the effective receptive field59, while\na sliding window with a fixed stride integrates broader spatial context. At the network output,\npatch-level predictions are confidence-weighted via monitored optical power and aggregated\nto yield the final image-level classification result (Fig.7 Inset c). Different output encoding\nschemes such as one-hot, binary or error-correcting output codes (ECOC)60 are viable for\nclassification where multi-class decisions are recovered by code-word decoding at inference54.\nWeight-sharing across network width enables patches to be evaluated in parallel61. Moreover,\noptical multiplexing62 such as wavelength or polarization, increases inference throughput and\nenables multi-task operation to be executed on the same chip (see Supplementary Section 2).\nBy using dual wavelengths, we demonstrate a 30 × 20 µm² single-chip design simultaneously\nclassifies MNIST and MedNIST datasets, achieving 95.1% and 98.0% test accuracy,\nrespectively (see Supplementary Section 2).\nA patch-efficient adjoint method (see Supplementary Section 3) is incorporated into the\ninverse-design of each PNN core, ensuring that the training cost is independent of the number\nof image patches. For a PNN core with input dimension 𝑁 and output dimension 𝐵, the number\nof FDTD solves per epoch is fixed at 𝑁+𝐵 under single-wavelength operation (see Methods).\nThis scaling keeps training tractable even for large-class datasets. For instance, with a 4 × 4\npatch kernel (𝑁 = 16) on ImageNet (1000 classes), where 𝐵 = 16 is determined by ECOC\ncoding54 with a 10-bit minimum code length plus 6 redundancy bits for error correction (see\nMethods), each PNN core design update requires only 32 FDTD-solves.\nFig. 7 Schematic diagram of PNN model scalable to larger neural networks. Images are segmented into\npatches (1…𝑘…𝑤) by a sliding window and encoded as feature vectors. Patches are processed by a sequence of\ninverse-designed PNN cores [PNN(\"),PNN($),… , PNN(%), PNN(%&\"),… , PNN(')] with interleaved\nphotodetector nonlinearities (NL). Width scalability is achieved by weight sharing across patches, while depth\nscalability arises from stacking cores with re-extraction of overlapping patches. Within each layer block, the first\nPNN core maps a feature vector to an optical-field distribution |𝐄| over 𝐵 channels for class evidence (Inset a)\nand forms the patch class-activation vector after photodetection. This vector is re-encoded onto the optical input\nports |𝐄′| of the second PNN (Inset b). These inputs are used to refine and embed class-prototype information in\npatch features before propagating to the next PNN core. Repeating this process sharpens per-patch class\nattributions, which are aggregated by confidence-weighted pooling into image-level predictions (Inset c). The\nsame forward operations occur during training. In addition, pooling is applied at intermediate layers for auxiliary\nsupervision.\nNew references (main manuscript):\n8. Filipovich, M. J. et al. Silicon photonic architecture for training deep neural networks\nwith direct feedback alignment. Optica 9, 1323-1332 (2022).\n11. Xu, D., Ma, Y., Jin, G. & Cao, L. Intelligent photonics: a disruptive technology to shape\nthe present and redefine the future. Engineering 46, 186-213 (2025).\n12. Momeni, A. et al. Training of physical neural networks. Nature 645, 53–61 (2025).\n13. Wetzstein, G. et al. Inference in artificial intelligence with deep optics and\nphotonics. Nature 588, 39–47 (2020).\n14. McMahon, P. L. The physics of optical computing. Nat. Rev. Phys. 5, 717–734 (2023).\n15. Gholami, A. et al. AI and Memory Wall. IEEE Micro 44, 33–39 (2024).\n54. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n55. Brendel, W. & Bethge, M. Approximating CNNs with bag-of-local-features models\nworks surprisingly well on ImageNet. In Proc. 7th International Conference on Learning\nRepresentations (ICLR, 2019).\nNL\nNetwork Depth\nSliding window\nPatch kernel\n|1| distribution\nClass 1 Class 2\nClass C\n|k| distribution |k|' distribution\nClass 1 Class 2\nClass C\nInput\nFeatures\n|w| distribution |w|' distribution\nRefined\nFeatures\ntupnI PNN\n(1)\nnoitacfiissalC\ntuptuO\nClass Evidence Feature Embedding\nClass 1 Class 2\nClass C\nPatch confidence\npooling\nImage-level\nclassification\nPatch\nPatch\nw\nClass\nFeature\nNL\nNetwork Width\nClass Decision\n|1|' distribution\nClass 1 Class 2\nClass C\nClass 1 Class 1\nClass 2 Class 2\nClass C Class C\nPNN\n(2)\nNL\nPNN\n(i)\nNL\nPNN\n(i+1)\nNL\nPNN\n(l)\n1\nk\nw\n1 2 C\nClass\nFeature\nClass\nFeature\npamtaeH\nssalC\nLayer Block\n56. Alzubaidi, L. et al. Review of deep learning: concepts, CNN architectures, challenges,\napplications, future directions. J. Big Data 8, 53 (2021).\n57. Hou, L. et al. Patch-based convolutional neural network for whole slide tissue image\nclassification. In Proc. IEEE Conference on Computer Vision and Pattern Recognition\n2424–2433 (IEEE, 2016).\n58. Choromanski, K. et al. Rethinking attention with performers. In Proc. 9th International\nConference on Learning Representations (ICLR, 2021).\n59. Luo, W., Li, Y., Urtasun, R. & Zemel, R. Understanding the effective receptive field in\ndeep convolutional neural networks. In Proc. 30th Conference on Neural Information\nProcessing Systems, 4905–4913 (NIPS, 2016).\n60. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\n61. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n62. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nOn Pages 8-9 (Supplementary information):\nSupplementary Section 2. Multiplexing method\nWeight sharing in the form of hardware/time multiplexing can be leveraged to scale the network\nin width13,14. Moreover, wavelength multiplexing15,16 can be incorporated into our inverse\ndesign process to enable multiple optical transformations of different dimensions to be realized\nwithin a single PNN. To demonstrate this flexibility, we design a single 30 × 20 𝜇m2 chip\ncapable of performing both MNIST and MedNIST classifications via wavelength multiplexing,\nachieving simulated test accuracy of 95.1% and 98.0%, respectively.\nAt each inverse-design update, a fixed number 𝑁′ = max(𝑁 ,𝑁 ) of independent\n! !\"#\nsimulations are run to compute the training forward-pass fields 𝐄 as well as 𝐶′ = max(𝐶 ,\n$%# !\n𝐶 ) independent reverse-mode simulations to compute the back-propagation fields, where\n!\"#\n𝑁 ,𝑁 are the input feature dimensions and 𝐶 , 𝐶 are the output class logit dimensions\n! !\"# ! !\"#\nfor MNIST and MedNIST datasets respectively. Computing the gradient fields requires a single\nmultiplication of the forward 𝐄 and adjoint 𝐄 fields per sample of the training dataset.\n$%# &#'\nFor MNIST, the 𝑁 = 10 input features are distributed across all input ports, with zero-padding\n!\napplied to any unused ports. Likewise, for MedNIST, the class outputs occupy the central six\nports, as shown in Supplementary Fig. 9. Supplementary Fig. 9a, b shows the field magnitude\n|𝐄| for random training samples in classes ‘Hand’ and ‘1’ for MedNIST and MNIST datasets\nmultiplexed at wavelengths 1570 nm and 1530 nm respectively, with the field hot-spot\nconcentration at the correct port. The confusion matrix in Supplementary Fig. 9c, d denotes\na test accuracy of 98.0% and 95.1% for MedNIST and MNIST datasets respectively,\ndemonstrating simultaneous classification performance for both datasets.\n'0'\n'1'\n'ChestCT' '2'\n'CXR' '3'\n'AbdomenCT' '4'\n'Hand' '5'\n'HeadCT' '6'\n'BreastMRI' '7'\n'8'\n'9'\nSupplementary Fig. 9 PNN multi-task classification. Field magnitude |𝐄| results for a sample from MedNIST\na and MNIST b datasets are presented, with a hotspot at the correct classification port corresponding to classes\n‘Hand’ at a wavelength of 1570 nm and ‘1’ at a wavelength of 1530 nm respectively. The confusion matrices in c\nand d shows a test accuracy of 98.0% and 95.1% for MedNIST and MNIST classes respectively, demonstrating\nstrong classification performance.\nNew references (Supplementary information):\n13. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n14. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n15. Cheung, A. K. C. et al. Inverse-designed CWDM demultiplexer operated in O-band.\nIn Optical Fiber Communications Conference, W1A-6 (IEEE, 2024).\n16. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nComments: Similar to the above, the authors simulate/design and fabricate two different\ndevices with different dimensions and different structure to tackle a similar classification task\n(MNIST and MedNIST). It would demonstrate the potential flexibility of the approach if the\nlarger (MedMNIST) was used in order to classify the MNIST task or vice versa. Such a small-\nscale experimental or numerical investigation would prove the adaptability of this approach;\notherwise, I am not confident about its use.\nReply:\nWe thank the reviewer for this comment. Our inverse-design framework is adaptable to multi-\ntask classification by optimizing a loss function that comprises of cross-entropy terms for each\ndataset-specific classifier. We demonstrate the flexibility of our approach for multi-task\nclassification by designing a single 30 × 20 𝜇m2 chip which performs both MNIST and\nMedNIST classification.\nWe have revised our manuscript to include discussion on the flexibility of our design method\nto multi-task classification and provided new simulation results and analysis to clarify this.\nOn Page 17:\nMoreover, optical multiplexing62 such as wavelength or polarization, increases inference\nthroughput and enables multi-task operation to be executed on the same chip (see\nSupplementary Section 2). By using dual wavelengths, we demonstrate a 30 × 20 µm² single-\nchip design which simultaneously classifies MNIST and MedNIST datasets, achieving 95.1%\nand 98.0% test accuracy, respectively (see Supplementary Section 2).\nNew reference (main manuscript):\n62. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nOn Pages 8-9 (Supplementary information):\nSupplementary Section 2. Multiplexing method\nWeight sharing in the form of hardware/time multiplexing can be leveraged to scale the network\nin width13,14. Moreover, wavelength multiplexing15,16 can be incorporated into our inverse\ndesign process to enable multiple optical transformations of different dimensions to be realized\nwithin a single PNN. To demonstrate this flexibility, we design a single 30 × 20 𝜇m2 chip\ncapable of performing both MNIST and MedNIST classifications via wavelength multiplexing,\nachieving simulated test accuracy of 95.1% and 98.0%, respectively.\nAt each inverse-design update, a fixed number 𝑁′ = max(𝑁 ,𝑁 ) of independent\n! !\"#\nsimulations are run to compute the training forward-pass fields 𝐄 as well as 𝐶′ = max(𝐶 ,\n$%# !\n𝐶 ) independent reverse-mode simulations to compute the back-propagation fields, where\n!\"#\n𝑁 ,𝑁 are the input feature dimensions and 𝐶 , 𝐶 are the output class logit dimensions\n! !\"# ! !\"#\nfor MNIST and MedNIST datasets respectively. Computing the gradient fields requires a single\nmultiplication of the forward 𝐄 and adjoint 𝐄 fields per sample of the training dataset.\n$%# &#'\nFor MNIST, the 𝑁 = 10 input features are distributed across all input ports, with zero-padding\n!\napplied to any unused ports. Likewise, for MedNIST, the class outputs occupy the central six\nports, as shown in Supplementary Fig. 9. Supplementary Fig. 9a, b shows the field magnitude\n|𝐄| for random training samples in classes ‘Hand’ and ‘1’ for MedNIST and MNIST datasets\nmultiplexed at wavelengths 1570 nm and 1530 nm respectively, with the field hot-spot\nconcentration at the correct port. The confusion matrix in Supplementary Fig. 9c, d denotes\na test accuracy of 98.0% and 95.1% for MedNIST and MNIST datasets respectively,\ndemonstrating simultaneous classification performance for both datasets.\n'0'\n'1'\n'ChestCT' '2'\n'CXR' '3'\n'AbdomenCT' '4'\n'Hand' '5'\n'HeadCT' '6'\n'BreastMRI' '7'\n'8'\n'9'\nSupplementary Fig. 9 PNN multi-task classification. Field magnitude |𝐄| results for a sample from MedNIST\na and MNIST b datasets are presented, with a hotspot at the correct classification port corresponding to classes\n‘Hand’ at a wavelength of 1570 nm and ‘1’ at a wavelength of 1530 nm respectively. The confusion matrices in c\nand d shows a test accuracy of 98.0% and 95.1% for MedNIST and MNIST classes respectively, demonstrating\nstrong classification performance.\nNew references (Supplementary information):\n13. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n14. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n15. Cheung, A. K. C. et al. Inverse-designed CWDM demultiplexer operated in O-band.\nIn Optical Fiber Communications Conference, W1A-6 (IEEE, 2024).\n16. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nComments: In lines 38–39, the authors claim that dense and large-scale chips should be\ndeveloped in order to address complex problems. On the other hand, in multiple cases, PNNs\nexploit time multiplexing and increased bandwidth capabilities so as to fit larger models on\nsmall-scale chips. In this context, it would enhance the quality of the manuscript if the authors\nincluded similar works and highlighted the merits of their approach compared to such schemes.\nReply:\nWe thank the reviewer for this comment. Our work focuses on a design regime optimized for\nsingle-shot inference. In this configuration, the inference latency is determined by the optical\npath length. There is no static tuning power, run-time weight memory loads for calibration, and\nenergy per inference is set by a one-pass drive of input modulators and readout detectors. These\nproperties are attractive for both ultra-low latency and energy efficiency at inference.\nMoreover, our design methodology is compatible with and extensible to multiplexing schemes\nincluding time-domain multiplexing21, where input data could be injected through the same\nports serialized in time to further scale the model complexity where required.\nWe have updated the manuscript to clarify this.\nOn Page 2:\nThis paradigm, analogous to application-specific integrated circuits (ASICs) in electronics,\neliminates static tuning power and runtime weight calibration13, leading to compact footprints\nwith ultra-low-latency stable inference14. Moreover, physically embedded weights realize in-\nmemory optical computing4, significantly mitigating memory-bandwidth constraints present at\nrun-time in conventional computing architectures at scale15.\nTo unlock the full potential of PNNs, it is essential to achieve high-density integration\ncapable of supporting complex models and large datasets on-chip3,6,16. Our approach\nemphasizes in-memory optical single-shot run-time operation, with highly efficient, stable and\nlow-latency inference.\nOn Pages 18-19:\nMultiplexing methods in the time63-65 or frequency-domain65 can be seamlessly integrated\ninto our inverse-design process to encode additional degrees of freedom, enhancing model\ncapacity within ultra-compact footprints. Moreover, incorporating high-speed modulation\nschemes, such as those based on the plasma dispersion effect66, electro-absorption67 and\nPockels effect68,69, along with fast photodetectors and analog-to-digital converters (ADCs)70\ncan facilitate high-bandwidth data interfacing for the PNNs.\nNew references (main manuscript):\n13. Wetzstein, G. et al. Inference in artificial intelligence with deep optics and\nphotonics. Nature 588, 39–47 (2020).\n14. McMahon, P. L. The physics of optical computing. Nat. Rev. Phys. 5, 717–734 (2023).\n15. Gholami, A. et al. AI and Memory Wall. IEEE Micro 44, 33–39 (2024).\n63. Lin, Z. et al. 120 GOPS photonic tensor core in thin-film lithium niobate for inference\nand in situ training.Nat. Commun.15, 9081 (2024).\n64. Song, Y. et al. Integrated electro-optic digital-to-analogue link for efficient computing and\narbitrary waveform generation. Nat. Photon. 19, 1107–1115 (2025).\n65. Ou, S. et al. Hypermultiplexed integrated photonics-based optical tensor processor. Sci.\nAdv. 11, eadu0228 (2025).\n66. De Marinis, L., Andriolli, N. & Contestabile, G. Analysis of integration technologies for\nhigh-speed analog neuromorphic photonics. IEEE J. Select. Top. Quant. Electron. 29\n(2023).\n67. Zhou, X., Yi, D., Chan, D. W. U. & Tsang, H. K. Silicon photonics for high-speed\ncommunications and photonic signal processing. npj Nanophotonics 1, 27 (2024).\n68. Hu, Y. et al. Integrated lithium niobate photonic computing circuit based on efficient and\nhigh-speed electro-optic conversion. Nat. Commun. 16, 8178 (2025).\n69. Zhang, S. et al. Thin-film lithium niobate photonic circuit for ray tracing acceleration.\nNat. Commun. 16, 5938 (2025).\n70. Shekhar, S. et al. Roadmapping the next generation of silicon photonics. Nat. Commun.\n15, 751 (2024).\nComments: Furthermore, the authors do not mention or analyze the bandwidth of the\namplitude-modulated signals used in experimental validation. This is a critical issue that\ngoverns the throughput of such PNNs. In addition, the authors should comment on whether\nmodulation speed affects performance and system training.\nReply:\nWe thank the reviewer for this comment. In our proof-of-concept experimental demonstration,\namplitude modulation is implemented using an MZI-modulator with a thermally tuned phase\nshifter on one arm. High-speed modulation schemes, such as those based on the plasma\ndispersion effect66, electro-absorption67 and Pockels effect68,69 have been experimentally\ndemonstrated. The optical PNN core itself is entirely passive and broadband. By integrating\nhigh-speed modulators along with fast photodetectors and analog-to-digital converters\n(ADCs)70, the system can be scaled to meet the demands of higher-speed operation.\nMoreover, since training is performed offline in the time-harmonic domain, the modulation\nspeed only governs inference throughput but does not affect the training process.\nWe have updated the manuscript to clarify this.\nOn Page 13:\nOptical signals from a CW laser source at 1550 nm are vertically coupled into the chip via\nvertical grating couplers (VGCs), distributed through a power splitter, and encoded with\namplitude and phase information for input data preparation via thermo-optic tuning (see\nMethods), and then processed via the inverse-designed PNN.\nOn Pages 18-19:\nMoreover, incorporating high-speed modulation schemes, such as those based on the plasma\ndispersion effect66, electro-absorption67 and Pockels effect68,69, along with fast photodetectors\nand analog-to-digital converters (ADCs)70 can facilitate high-bandwidth data interfacing for\nthe PNNs.\nNew references (main manuscript):\n66. De Marinis, L., Andriolli, N. & Contestabile, G. Analysis of integration technologies for\nhigh-speed analog neuromorphic photonics. IEEE J. Select. Top. Quant. Electron. 29\n(2023).\n67. Zhou, X., Yi, D., Chan, D. W. U. & Tsang, H. K. Silicon photonics for high-speed\ncommunications and photonic signal processing. npj Nanophotonics 1, 27 (2024).\n68. Hu, Y. et al. Integrated lithium niobate photonic computing circuit based on efficient and\nhigh-speed electro-optic conversion. Nat. Commun. 16, 8178 (2025).\n69. Zhang, S. et al. Thin-film lithium niobate photonic circuit for ray tracing acceleration.\nNat. Commun. 16, 5938 (2025).\n70. Shekhar, S. et al. Roadmapping the next generation of silicon photonics. Nat. Commun.\n15, 751 (2024).\nComments: In line 51, the authors mention that high contrast materials and tight\nconfinement increase signal expressivity. The statement is not clear—what do the authors\ndefine as signal expressivity, and how exactly is it affected by light confinement?\nReply:\nWe thank the reviewer for this comment. By signal expressivity, we refer to the diversity of\ninput-dependent field transformations that a passive photonic structure can implement within\na bounded geometry21. Index contrast governs optical confinement and therefore the spatial\nfrequency bandwidth accessible inside the device27,28. High-contrast platforms support\nstronger confinement, enabling steeper amplitude and phase variations via enhanced local\nscattering27. As a result, a richer set of optical input-output transformations can be realized\nper unit volume. This enhanced expressivity underpins functionally dense integrated photonic\ncircuits, where light propagation is synthesized via local refractive index variations.\nWe have revised this sentence in the manuscript to clarify this definition.\nOn Page 4:\nMoreover, increasing the refractive index contrast enhances signal expressivity, which refers\nto the capacity for a rich set of field transformations within a fixed geometry21, through\nstrengthening confinement, supporting internal resonances26, and enabling stronger\ninterference and engineered scattering effects27,28 - all of which are critical for high-fidelity\noptical computation.\nHere, we demonstrate ultra-compact, inverse-designed PNN accelerators fabricated on a\nhigh-index-contrast silicon-on-insulator (SOI) platform operating at a single wavelength,\nwhere higher index contrast strengthens light confinement and supports rich interference and\nmode-mixing dynamics27,28, thereby expanding the space of attainable optical responses per\nunit volume.\nNew reference (main manuscript):\n27. Meng, Y. et al. Optical meta-waveguides for integrated photonics and beyond. Light Sci.\nAppl. 10, 235 (2021).\nComments: Another fundamental question is that in most photonic-based computing,\ninference comes at low latency, marginal power consumption, and increased throughput, while\ntraining in most cases entails elaborate simulations using a digital twin approach. In these cases\n(e.g. MZI meshes), the models used can be downgraded to simple transfer function combinations\n(scattering matrix). On the other hand, in this work, inference is equally straightforward, but\ntraining-designing entails the use of costly physically accurate models. It would be beneficial\nif the authors could comment on this and provide a table including the number of parameters\n(multiplications) needed to train (offline) different PNNs compared to their approach.\nReply:\nWe thank the reviewer for the comments. Our proposed inverse-design framework for\nphotonic-based computing differs from both reconfigurable MZI-mesh architectures and on-\nchip diffractive optical neural networks (DONNs). It employs computational optimization to\nexplore an expansive device design space for the optimal structure, where each voxel in the\ndesign region contributes to the global scattering, interference and phase evolution of light\ninside device21. Consequently, the proposed PNN architecture achieves a high trainable\nparameter density within a compact footprint. In the following, we compare representative\nconfigurations of reconfigurable MZI-mesh networks, DONNs, and our inverse-designed PNN\nfor the MNIST classification task with respect to both computational density and\nmultiplications required for training.\nAn 𝑁 -port MZI mesh topology implements an arbitrary 𝑁 ×𝑁 unitary matrix\n( ( (\ntransformation on 𝑁 -dimensional input vectors and typically requires 𝒪(𝑁) ) trainable\n( (\nparameters10. Given each trainable parameter is implemented as a discrete phase shifter within\nan MZI unit3,10, the overall footprint scales rapidly with network size, typically reaching\nmillimeter-scale3, which in turn constrains the integration density. Triangular mesh networks10\nhave been theoretically proposed for MNIST classification, comprising two hidden layers of\nunitary 64 × 64 matrix multiplications. The number of equivalent digital forward pass\nmultiplications 𝑀 can be approximated by the floating-point operation equation7 𝑀 ≈\n$ $\nm×𝑁), where 𝑚 is the number of network layers - including the hidden layers and the output\n(\nlayer, equating to approximately 1.23 × 104 forward pass multiplications. The gradient\ncomputation in the backward pass requires approximately double the multiplications for the\nforward pass, owing to additional matrix–vector operations needed to propagate and update the\ngradients1.\nOn-chip diffractive optical neural networks (DONNs) encode trainable parameters as phase\nvalues in compact subwavelength units (SWUs)3,6. Computation in such systems relies\nexclusively on feedforward light propagation through sequential, layer-by-layer\narchitectures3,6,22. The spacing between consecutive layers is governed by diffraction\nconstraints22, which restricts how densely the architecture can be scaled3,22. For an on-chip\nDONN system comprising 𝑚 layers, each with 𝑁 trainable parameters (number of SWUs)\n#\nperforming an 𝑁 ×𝑁 weighting matrix, the number of multiplication operations can be\n# #\napproximated as 𝑀 ≈ 𝑚×𝑁). For instance, a DONN6 chip experimentally demonstrates the\n$ #\nMNIST classification task, whose network architecture is composed of three HLs and a single\noutput layer. Each HL comprises 70 SWUs, with 250 𝜇m distance between successive HLs\nalong the 𝑥-axis and a network width of 105 𝜇m along the 𝑦-axis6. The DONN also requires a\n250 𝜇m propagation distance from the input waveguides to the first HL, and a further 250 𝜇m\ndistance from the last HL to the photodetectors6. Together, these dimensions yield a total\nfootprint of 0.105 mm², and correspondingly, a computational density of approximately 2000\ntrainable parameters per mm². The number of forward multiplications7 to train the device for\nthis task is approximately 1.96 × 104. Similarly, the gradient computation in the backward pass\nrequires approximately twice the number of matrix-vector multiplications of the forward pass1.\nIn contrast, during our inverse-design process, we perform one such multiplication of the\nforward pass (𝐄 ) and backward pass (𝐄 ) 𝐄-fields per design voxel to update its value via\n$%# &#'\ngradient-based optimization. Based on our 3D-FDTD planar design region discretized into 𝑛\nvoxels, the number of trainable parameters is approximately 1.6 × 105 for the MNIST chip.\nAs a result, the number of complex elementwise multiplications 𝑀 per gradient update is\n∇\ncalculated as 𝑀 = 𝑛×𝑑′, where 𝑑′ = 3 denotes the three components of the electric-field\n∇\nvector per voxel. This corresponds to approximately 4.8 × 105 multiplications per update to\ntrain the MNIST classification task. Given an MNIST design footprint of 20 × 20 𝜇m2, the\ncomputational density of this chip is approximately 400 million trainable parameters per mm2,\nwhich demonstrates this method’s potential for significant integration density in comparison to\nrepresentative DONN works4-6. Importantly, by virtue of the adjoint variable method, which\nobtains a gradient of all design voxels from a single overlap of the 𝐄 and 𝐄 fields21, these\n$%# &#'\nmultiplications parallelize naturally across voxels on GPUs, enabling fast, scalable gradient\nupdates across a high-dimensional parameter space.\nAs an additional demonstration of the computational scale of our inverse-design framework,\nwe calculate the number of parameters in the 30 × 20 𝜇m2 MedNIST chip as 2.4 × 105, which\nequates 7.2 × 105 multiplications per design update during training and the same computational\ndensity as the MNIST classification task, given the same spatial resolution.\nAt inference, in contrast to reconfigurable MZI approaches, by physically embedding the\nsystem weights into the design region itself, the necessity for run-time network parameter\nloading and calibration is removed4, greatly reducing the memory-bandwidth bottleneck\npresent in conventional computing schemes.\nWe have updated the manuscript to clarify the compute capacity of our approach in comparison\nto other PNN methods, as well as a respective comparison of training parameter density and\nmultiplications in a new supplementary section “PNN comparative computational density\ncalculation”. We also argue our approach at inference has desirable advantages compared to\nexisting reconfigurable methods along with support references.\nOn Page 1:\nEach subwavelength voxel in the inverse-designed PNN serves as a trainable degree of freedom,\nyielding a computational density of approximately 400 million parameters per mm².\nOn Pages 2-3:\nPNN accelerator hardware based on programmable photonic integrated circuits such as Mach-\nZehnder interferometer (MZI) meshes7, microrings8, and waveguide attenuators9, enables in\nsitu training10 and dynamic weight adjustment through optical backpropagation10,11. By\ncontrast, task-optimized architectures are trained offline12 and then realized in hardware for\nstable and efficiency-critical inference6. This paradigm, analogous to application-specific\nintegrated circuits (ASICs) in electronics, eliminates static tuning power and runtime weight\ncalibration13, leading to compact footprints with low-latency, stable inference14. Moreover,\nphysically embedded weights enable in-memory optical computing4, significantly alleviating\nmemory-bandwidth constraints present at run-time in conventional computing architectures at\nscale15.\nUnlike diffractive models4-6,22, which rely on sequential, layer-by-layer architectures\nmimicking feedforward deep neural networks (DNNs), topology optimization allows arbitrary\ncontrol over scattering, interference, and phase evolution throughout the device’s entire volume,\nenabling more compact and functionally dense implementations23-25.\nOn Page 4:\nEach subwavelength voxel acts as a trainable degree of freedom, yielding an ultra-compact\nrepresentation of optical transformations with a computational density of approximately 400\nmillion parameters per mm².\nOn Page 11:\nThe MNIST and MedNIST PNNs comprise approximately 1.6 × 105 and 2.4 × 105 trainable\nparameters, respectively (see Supplementary Section 1). During each training epoch, the\ncorresponding multiplication operations are parallelized on GPUs, enabling fast and scalable\ngradient updates. The computational density reaches approximately 400 million trainable\nparameters per mm² (Supplementary Section 1), highlighting our PNNs’ capability to achieve\nhigh computational density within an ultra-compact footprint.\nNew references (main manuscript):\n8. Filipovich, M. J. et al. Silicon photonic architecture for training deep neural networks\nwith direct feedback alignment. Optica 9, 1323-1332 (2022).\n11. Xu, D., Ma, Y., Jin, G. & Cao, L. Intelligent photonics: a disruptive technology to shape\nthe present and redefine the future. Engineering 46, 186-213 (2025).\n12. Momeni, A. et al. Training of physical neural networks. Nature 645, 53–61 (2025).\n13. Wetzstein, G. et al. Inference in artificial intelligence with deep optics and\nphotonics. Nature 588, 39–47 (2020).\n14. McMahon, P. L. The physics of optical computing. Nat. Rev. Phys. 5, 717–734 (2023).\n15. Gholami, A. et al. AI and Memory Wall. IEEE Micro 44, 33–39 (2024).\n22. Wu, Z., Zhou, M., Khoram, E., Liu, B. & Yu, Z. Neuromorphic metasurface. Photon.\nRes. 8, 46–50 (2020).\nOn Page 8 (Supplementary information):\nSupplementary Section 1. PNN comparative estimate of computational density and\nmultiplications\nThe following presents a comparative estimate of the computational density and number of\nmultiplications required to train different on-chip PNN approaches.\nAn 𝑁 -port MZI mesh topology implements an arbitrary 𝑁 ×𝑁 unitary matrix\n( ( (\ntransformation on 𝑁 -dimensional input vectors and typically requires 𝒪(𝑁) ) trainable\n( (\nparameters5,6. For example, triangular mesh networks theoretically proposed for MNIST\nclassification comprises two hidden layers (HLs) of unitary 64 × 64 matrix multiplications5.\nThe number of equivalent digital forward pass multiplications 𝑀 can be approximated by the\n$\nfloating-point operation equation7 𝑀 ≈ 𝑚×𝑁) , where 𝑚 is the number of network layers,\n$ (\nincluding the hidden layers and the output layer8, which equates to approximately 1.23 × 104\nforward pass multiplications. The gradient computation in the backward pass requires\napproximately double the multiplications of the forward pass, owing to additional matrix–\nvector operations needed to propagate and update the gradients9.\nOn-chip diffractive optical neural networks (DONNs) encode trainable parameters as phase\nvalues in subwavelength units (SWUs)10,11. For an on-chip DONN system with 𝑚 layers, each\ncomprising 𝑁 trainable parameters (number of SWUs) performing an 𝑁 ×𝑁 weighting\n# # #\nmatrix, the total number of forward pass multiplications can be approximated as 𝑀 ≈\n$\nm×𝑁) 8. For instance, a DONN11 chip experimentally demonstrates the MNIST classification\n#\ntask, whose network architecture is composed of three HLs and a single output layer. Each HL\ncomprises 70 SWUs, with 250 𝜇m distance between successive HLs along the 𝑥-axis and a\nnetwork width of 105 𝜇m along the 𝑦-axis11. The DONN also requires a 250 𝜇m propagation\ndistance from the input waveguides to the first HL, and a further 250 𝜇m distance from the last\nHL to the photodetectors11. Together, these dimensions yield a total footprint of 0.105 mm²,\nand correspondingly, a computational density of approximately 2000 trainable parameters per\nmm². The number of forward multiplications to train the device for this task is approximately\n1.96 × 104. Similarly, the gradient computation in the backward pass requires approximately\ntwice the number of matrix-vector multiplications of the forward pass9.\nIn contrast, during our inverse-design process, we perform one such multiplication of the\nforward pass (𝐄 ) and backward pass (𝐄 ) 𝐄-fields per design voxel to update its value via\n$%# &#'\ngradient-based optimization. The 3D-FDTD planar design region (20 × 20 μm²) of the MNIST\nPNN is discretized with a spatial resolution of approximately 𝜆/30, where 𝜆 is the operating\nwavelength. Under this configuration, the number of design voxels (𝑛), each representing the\nlocal refractive index and serving as a trainable parameter, is approximately 1.6 × 105.\nAccordingly, the number of complex elementwise multiplications 𝑀 per gradient update is\n∇\ncalculated as 𝑀 = 𝑛×𝑑′, where 𝑑′ = 3 denotes the three components of the electric-field\n∇\nvector per voxel. This corresponds to approximately 4.8 × 105 multiplications per update to\ntrain the MNIST classification task. The resulting computational density of this chip is\napproximately 400 million trainable parameters per mm2, which demonstrates this method’s\ncapability for significant integration density. Importantly, by virtue of the adjoint variable\nmethod, which obtains a gradient of all design voxels from a single overlap of the 𝐄 and\n$%#\n𝐄 fields12, these multiplications parallelize across voxels on graphics processing units\n&#'\n(GPUs), enabling fast, scalable gradient updates. As an additional demonstration, we calculate\nthe number of trainable parameters in the 30 × 20 𝜇m2 MedNIST chip as 2.4 × 105, which\ncorresponds to approximately 7.2 × 105 multiplications per design update during training.\nSince the design region spatial resolution is equal to that of the MNIST chip, this also equates\nto approximately 400 million trainable parameters per mm2.\nNew References (Supplementary information):\n5. Pai, S. et al. Experimentally realized in situ backpropagation for deep learning in\nphotonic neural networks. Science 380, 398–404 (2023).\n6. Wu, B. et al. Real-valued optical matrix computing with simplified MZI mesh. Intell.\nComput. 2, 0047 (2023).\n7. Shen, Y. et al. Deep learning with coherent nanophotonic circuits. Nat. Photon. 11,\n441–446 (2017).\n8. Cheng, J. et al. Multimodal deep learning using on-chip diffractive optics with in situ\ntraining capability. Nat. Commun. 15, 6189 (2024).\n9. LeCun, Y., Bengio, Y. & Hinton, G. Deep learning. Nature 521, 436–444 (2015).\n10. Fu, T. et al. Optical neural networks: progress and challenges. Light Sci. Appl. 13, 263\n(2024).\n11. Fu, T. et al. Photonic machine learning with on-chip diffractive optics. Nat.\nCommun. 14, 70 (2023).\n12. Molesky, S. et al. Inverse design in nanophotonics. Nat. Photon. 12, 659–670 (2018).\nComments: In line 96, the authors mention that their design allows spatially adaptive\ntransport. What do they mean by this? How is adaptability achieved in such a structure?\nReply:\nWe thank the reviewer for this comment. Spatially adaptive light transport refers to the\npassive, input-conditioned routing and transformation of optical power within the inverse-\ndesigned structure21. The adaptivity of light transport to different input mode excitations\narises from the optimized material index distribution, which establishes interference and\nmultiple scattering pathways that vary with the incident spatial distribution28.\nSpecifically, the PNNs are inverse-designed so that combinations of 𝑁-input spatial modes\nexcite distinct superpositions of scattering channels, steering optical power toward output\nwaveguides to encode class information. In this way, the structure achieves single-shot\nclassification across a diverse dataset of class features (Supplementary Fig. 4) within a\ncompact footprint.\nWe have updated the main manuscript to clarify this definition.\nOn Page 5:\nWithin the topology-optimized scattering region, the encoded optical fields undergo complex\ninterference and scattering interactions throughout the device volume. Within the topology-\noptimized scattering region, the encoded optical fields undergo complex interference and\nscattering interactions throughout the device volume. The scattering region itself is designed\nto support non-uniform optical transport21,28, with routing behavior that adapts to the spatial\ndistribution of optical inputs, efficiently redistributing optical power from the 𝑁 input\nencoding ports toward the 𝐶 class output ports. We denote 𝐶 as the number of classes in the\ndataset and 𝐵 as the output dimension of the PNN, where 𝐵 = 𝐶 in the one-hot setting.\nComments: In line 98, the authors state that the output of a Perceptron using PNN is a\nprobability function. That is not correct, as in order to transform optical outputs to probabilities,\nphotodetection and digital post-processing (e.g. softmax) should be applied.\nReply:\nWe thank the reviewer for the comments and have revised this statement in the main manuscript.\nOn Page 5:\nThe PNN outputs represent class-evidence scores. A probability distribution across\nclassification categories is obtained after photodetection, analogous to the output layer of a\nmulti-layer perceptron network (Fig. 1b).\nComments: In lines 120–123, the authors describe briefly their strategy for minimizing\ncomputational resources. It is my belief that this section is not clear and needs improvement.\nIf I am not mistaken, each input, instead of being injected at each port for the simulation, is\ninjected only to one. If this is correct, does it impact the final accuracy? If the authors used a\nfull-scale training approach, would performance be affected, and if yes, how much?\nReply:\nWe thank the reviewer for this comment. To clarify, we perform one independent 3D-FDTD\nsimulation for each of the 𝑁-input ports, where 𝑁 is the feature dimension. These solves yield\na set of basis fields, one per input port as shown in columns 1 and 2 of the inset in Fig.1c. Using\nthese basis fields, we synthesize any arbitrary data sample’s response as a weighted sum in the\nlinear time-harmonic domain as shown in column 3 of the inset in Fig. 1c. This decouples the\ndataset length (𝐿) from the simulation cost as forward-pass fields can be computed\nalgebraically without additional Maxwell solves. This method is numerically equivalent to full-\nscale training requiring 𝐿-independent FDTD simulations, which is computationally intractable\nfor large-volume datasets.\nWe have updated the main manuscript to clarify our method for minimizing simulation cost.\nOn Page 6:\nTo simulate the forward-pass optical fields for a dataset containing 𝐿 samples, each comprising\n𝑁 features, we exploit the linearity of Maxwell's equations to significantly reduce\ncomputational complexity. Rather than performing 𝐿 separate 3D-FDTD simulations, we\nreduce the simulation cost by executing one independent simulation for each of the 𝑁-input\nencoding ports, as illustrated in columns 1 and 2 of the inset in Fig. 1c.\nComments: The authors generate two different chips to tackle two different tasks; the chips\nhave different areas and offer different performance. Should the size of the chip be scaled with\nthe task’s complexity? And if yes, what type of chip would be needed for an even more\ncomplicated task (e.g. ImageNet)? It is my opinion that the authors should provide at least\nnumerical results where they tackle MNIST with chips of different area/complexity so as to\nhighlight this dependence.\nReply:\nWe thank the reviewer for this comment. While task complexity and chip size are inherently\nrelated, the relationship is not strictly linear. For example, in our newly added section, “PNN\nmodel scalability to larger neural networks”, we discuss how stacking modularized PNN cores,\nweight sharing and pooling strategies can expand network capacity without relying solely on\nlinear increases in chip footprint. This new section clarifies why scaling in our architecture is\nmore nuanced and efficient and examines how model complexity evolves when extended to\nlarger datasets and multi-layer configurations.\nOn Pages 16-18:\nPNN model scalability to larger neural networks. We scale the model to larger neural\nnetworks by stacking PNNs with interleaved photodetector nonlinearity54, where each PNN is\ninverse-designed for task-optimized optical-field transformations. Figure 7 illustrates the\nschematic architecture for the PNN model scalable to larger neural networks. Input images are\nsegmented into patches, encoded as feature vectors, and processed by pairs of consecutive PNN\ncores. The first PNN maps input features to an optical-field distribution across class channels\n(Fig.7 Inset a), which are converted after photodetection into a series of patch class-activation\nvectors55. The second PNN re-encodes this vector into refined patch features (Fig.7 Inset b),\nwhich become increasingly aligned with class prototypes as they propagate in the next stage.\nEach pair is analogous to a layer block in multilayer DNNs56, where blocks are stacked to\nextract and refine feature representations57. Successive PNN blocks progressively refine these\nclass-level embeddings, filtering out patches with low confidence e.g. background, and\nsharpening those with high-class correlation, generating richer interpretability of feature-to-\nclass associations58. This depth-wise refinement expands the effective receptive field59, while\na sliding window with a fixed stride integrates broader spatial context. At the network output,\npatch-level predictions are confidence-weighted via monitored optical power and aggregated\nto yield the final image-level classification result (Fig.7 Inset c). Different output encoding\nschemes such as one-hot, binary or error-correcting output codes (ECOC)60 are viable for\nclassification where multi-class decisions are recovered by code-word decoding at inference54.\nWeight-sharing across network width enables patches to be evaluated in parallel61. Moreover,\noptical multiplexing62 such as wavelength or polarization, increases inference throughput and\nenables multi-task operation to be executed on the same chip (see Supplementary Section 2).\nBy using dual wavelengths, we demonstrate a 30 × 20 µm² single-chip design simultaneously\nclassifies MNIST and MedNIST datasets, achieving 95.1% and 98.0% test accuracy,\nrespectively (see Supplementary Section 2).\nA patch-efficient adjoint method (see Supplementary Section 3) is incorporated into the\ninverse-design of each PNN core, ensuring that the training cost is independent of the number\nof image patches. For a PNN core with input dimension 𝑁 and output dimension 𝐵, the number\nof FDTD solves per epoch is fixed at 𝑁+𝐵 under single-wavelength operation (see Methods).\nThis scaling keeps training tractable even for large-class datasets. For instance, with a 4 × 4\npatch kernel (𝑁 = 16) on ImageNet (1000 classes), where 𝐵 = 16 is determined by ECOC\ncoding54 with a 10-bit minimum code length plus 6 redundancy bits for error correction (see\nMethods), each PNN core design update requires only 32 FDTD-solves.\nFig. 7 Schematic diagram of PNN model scalable to larger neural networks. Images are segmented into\npatches (1…𝑘…𝑤) by a sliding window and encoded as feature vectors. Patches are processed by a sequence of\ninverse-designed PNN cores [PNN(\"),PNN($),… , PNN(%), PNN(%&\"),… , PNN(')] with interleaved\nphotodetector nonlinearities (NL). Width scalability is achieved by weight sharing across patches, while depth\nscalability arises from stacking cores with re-extraction of overlapping patches. Within each layer block, the first\nPNN core maps a feature vector to an optical-field distribution |𝐄| over 𝐵 channels for class evidence (Inset a)\nand forms the patch class-activation vector after photodetection. This vector is re-encoded onto the optical input\nports |𝐄′| of the second PNN (Inset b). These inputs are used to refine and embed class-prototype information in\npatch features before propagating to the next PNN core. Repeating this process sharpens per-patch class\nattributions, which are aggregated by confidence-weighted pooling into image-level predictions (Inset c). The\nsame forward operations occur during training. In addition, pooling is applied at intermediate layers for auxiliary\nsupervision.\nNew references (main manuscript):\nNL\nNetwork Depth\nSliding window\nPatch kernel\n|1| distribution\nClass 1 Class 2\nClass C\n|k| distribution |k|' distribution\nClass 1 Class 2\nClass C\nInput\nFeatures\n|w| distribution |w|' distribution\nRefined\nFeatures\ntupnI PNN\n(1)\nnoitacfiissalC\ntuptuO\nClass Evidence Feature Embedding\nClass 1 Class 2\nClass C\nPatch confidence\npooling\nImage-level\nclassification\nPatch\nPatch\nw\nClass\nFeature\nNL\nNetwork Width\nClass Decision\n|1|' distribution\nClass 1 Class 2\nClass C\nClass 1 Class 1\nClass 2 Class 2\nClass C Class C\nPNN\n(2)\nNL\nPNN\n(i)\nNL\nPNN\n(i+1)\nNL\nPNN\n(l)\n1\nk\nw\n1 2 C\nClass\nFeature\nClass\nFeature\npamtaeH\nssalC\nLayer Block\n54. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n55. Brendel, W. & Bethge, M. Approximating CNNs with bag-of-local-features models\nworks surprisingly well on ImageNet. In Proc. 7th International Conference on Learning\nRepresentations (ICLR, 2019).\n56. Alzubaidi, L. et al. Review of deep learning: concepts, CNN architectures, challenges,\napplications, future directions. J. Big Data 8, 53 (2021).\n57. Hou, L. et al. Patch-based convolutional neural network for whole slide tissue image\nclassification. In Proc. IEEE Conference on Computer Vision and Pattern Recognition\n2424–2433 (IEEE, 2016).\n58. Choromanski, K. et al. Rethinking attention with performers. In Proc. 9th International\nConference on Learning Representations (ICLR, 2021).\n59. Luo, W., Li, Y., Urtasun, R. & Zemel, R. Understanding the effective receptive field in\ndeep convolutional neural networks. In Proc. 30th Conference on Neural Information\nProcessing Systems, 4905–4913 (NIPS, 2016).\n60. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\n61. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n62. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nMore generally, we observe that task complexity correlates more with degrees of freedom\npresent in the design. Increasing the footprint is one way to achieve this, another way is\nexpanding capacity via multiplexing. We demonstrate this concept by performing classification\nof two tasks of different dimension on the same chip without sacrificing device footprint. This\nshows that even with a fixed chip area, spectral multiplexing can increase effective task\ncapacity.\nOn Pages 8-9 (Supplementary information):\nSupplementary Section 2. Multiplexing method\nWeight sharing in the form of hardware/time multiplexing can be leveraged to scale the network\nin width13,14. Moreover, wavelength multiplexing15,16 can be incorporated into our inverse\ndesign process to enable multiple optical transformations of different dimensions to be realized\nwithin a single PNN. To demonstrate this flexibility, we design a single 30 × 20 𝜇m2 chip\ncapable of performing both MNIST and MedNIST classifications via wavelength multiplexing,\nachieving simulated test accuracy of 95.1% and 98.0%, respectively.\nAt each inverse-design update, a fixed number 𝑁′ = max(𝑁 ,𝑁 ) of independent\n! !\"#\nsimulations are run to compute the training forward-pass fields 𝐄 as well as 𝐶′ = max(𝐶 ,\n$%# !\n𝐶 ) independent reverse-mode simulations to compute the back-propagation fields, where\n!\"#\n𝑁 ,𝑁 are the input feature dimensions and 𝐶 , 𝐶 are the output class logit dimensions\n! !\"# ! !\"#\nfor MNIST and MedNIST datasets respectively. Computing the gradient fields requires a single\nmultiplication of the forward 𝐄 and adjoint 𝐄 fields per sample of the training dataset.\n$%# &#'\nFor MNIST, the 𝑁 = 10 input features are distributed across all input ports, with zero-padding\n!\napplied to any unused ports. Likewise, for MedNIST, the class outputs occupy the central six\nports, as shown in Supplementary Fig. 9. Supplementary Fig. 9a, b shows the field magnitude\n|𝐄| for random training samples in classes ‘Hand’ and ‘1’ for MedNIST and MNIST datasets\nmultiplexed at wavelengths 1570 nm and 1530 nm respectively, with the field hot-spot\nconcentration at the correct port. The confusion matrix in Supplementary Fig. 9c, d denotes\na test accuracy of 98.0% and 95.1% for MedNIST and MNIST datasets respectively,\ndemonstrating simultaneous classification performance for both datasets.\n'0'\n'1'\n'ChestCT' '2'\n'CXR' '3'\n'AbdomenCT' '4'\n'Hand' '5'\n'HeadCT' '6'\n'BreastMRI' '7'\n'8'\n'9'\nSupplementary Fig. 9 PNN multi-task classification. Field magnitude |𝐄| results for a sample from MedNIST\na and MNIST b datasets are presented, with a hotspot at the correct classification port corresponding to classes\n‘Hand’ at a wavelength of 1570 nm and ‘1’ at a wavelength of 1530 nm respectively. The confusion matrices in c\nand d shows a test accuracy of 98.0% and 95.1% for MedNIST and MNIST classes respectively, demonstrating\nstrong classification performance.\nNew references (Supplementary information):\n13. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n14. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n15. Cheung, A. K. C. et al. Inverse-designed CWDM demultiplexer operated in O-band.\nIn Optical Fiber Communications Conference, W1A-6 (IEEE, 2024).\n16. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nFinally, the use of binary-coded or error correcting output codes (ECOC)61 also provides a\nmethod to expand degrees of freedom to avoid a linear increase in the number of output ports\nfor classification. Recent work55 has experimentally classified the 1623-class Omniglot1\ndataset using a limited number (𝐶=17) output ports with ECOC encoding. Combined with\nmultiplexing, for example in the frequency domain62, this provides a practical pathway for\nscaling to much larger datasets such as ImageNet. We provide details on expected simulation\nand computational costs for high-class capacity classification tasks in the new methods\nsection “Training framework and scalability”.\n1. Lake, B. M., Salakhutdinov, R. & Tenenbaum, J. B. Human-level concept learning\nthrough probabilistic program induction. Science 350, 1332–1338 (2015).\nOn Pages 21-22:\nTraining framework and scalability. Each PNN core is trained to implement task-specific\noptical-field transformations by minimizing a CE loss1. Class information is embedded through\none-hot, binary or ECOC coding60, with ECOC advantageous for large-class datasets by\nreducing the number of output channels required54. A patch-efficient adjoint method (see\nSupplementary Section 3) is incorporated into the optimization, whereby patch contributions\nare aggregated algebraically before gradient computation. This ensures that the total number\nof required FDTD solves per epoch for one per PNN core depends only on the input and output\ndimensions (𝑁 + 𝐵) and is independent of the number of image patches (see Supplementary\nSection 3). The optimization incorporates nonlinear photodetection, which converts optical\nlogits into class-activation vectors and provides interpretable measures of patch-to-class\nassociation55. To stabilize training and ensure global consistency, auxiliary supervision can be\nintroduced by aggregating patch-level predictions with confidence weighting at intermediate\nlayers, enforcing agreement with image-level classification. Training is proceeded block by\nblock, where each consecutive pair of PNN cores are optimized jointly and then frozen after\nconvergence. Overlapping patches are re-extracted and passed to the next stage, enabling depth\nscalability through successive refinement. Reuse of pre-trained shallow blocks as priors can\nfurther improve efficiency by reducing the number of gradient updates needed for deeper\nlayers81. The final PNN core maps refined patch features into class-activations, and confidence-\nweighted pooling across patches produces the image-level prediction. All optimizations are\nperformed in simulation, and the resulting geometries are used directly for inference in the\nfabricated device.\nNew references (main manuscript):\n54. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n55. Brendel, W. & Bethge, M. Approximating CNNs with bag-of-local-features models\nworks surprisingly well on ImageNet. In Proc. 7th International Conference on Learning\nRepresentations (ICLR, 2019).\n60. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\n81. Chen, T., Goodfellow, I. & Shlens, J. Net2net: accelerating learning via knowledge\ntransfer. In Proc. International Conference on Learning Representations (ICLR, 2016).\nComments: Accuracy for the MNIST dataset is in the order of 97.8% when simulations are\nconsidered. There are effects that are not discussed within the manuscript regarding these\nsimulations, e.g. the inclusion of thermal/shot noise at the photodiodes, limited bit-resolution\nat the ADCs, nonlinearity or imperfections at the amplitude modulator, bandwidth limitations\nof equipment etc. All these effects are probably absent in the numerical simulations, thus the\nreally high accuracy of 97.8%. On the other hand, during experimental validation, even when\nusing a small number of samples, accuracy radically drops to 89%. This last value is not really\nimpressive, based on the fact that 96%–97% accuracy has been demonstrated with PNNs in\nthe past. It would be beneficial if the authors could comment on the sensitivity of their design\nto actual experimental conditions (explain this performance degradation) and compare their\napproach to state-of-the-art PNNs targeting MNIST.\nReply:\nWe thank the reviewer for their comments. In numerical simulation, we characterize our design\nusing the energy density matrix, which quantifies the output power distribution across different\ndatasets. As shown in Fig. 2d and 3d strong energy localization at the target class output port\n,\nis achieved, indicating the device’s robustness against potential perturbations such as thermal\nand shot noise at the photodiodes. The experiment results in Fig. 6c and 6d demonstrate mean\neffective localization of 50.5% of the optical energy at the output ports associated with each\ntarget class with only a 16.1% variation in maximum output power for MedNIST classification.\nThis underscores the device’s consistency and resilience to perturbation under experimental\nconditions, as deviations in output port selection remain minimal across datasets.\nOur experimental setup is shown in Supplementary Figure 2. Input light at 1550 nm from a\ntunable laser (Keysight N7778C) with a relative intensity noise (RIN) less than -150 dB/Hz77\nis fiber-coupled onto the chip. The optical power at each output channel of the PNN is\nmeasured using multiport optical power meters (Keysight N7744C) that incorporate built-in\naveraging78. This cumulative average time of signals suppress random thermal and shot noises\nat the photodiode readout79, thereby minimizing their impact on system performance. The\noptical power meter (Keysight N7744C) has a dynamic range of -65 dB (exceeding 80 dB with\nrange switching)79 which implies the ADC resolution is sufficiently high. Hence, we do not\nexpect the ADC bit depth to be the cause of the observed performance degradation.\nFor experimental data preparation, we measure analogue optical input power to the PNN\naccelerators, using a random stratified subset of 10 samples per class. Supplementary Figures\n4 and 5 illustrate a global mean average error (MAE) of 0.0277 for MNIST and 0.0310 for\nMedNIST with respect to the simulated inputs. These low MAE values confirm a strong\ncorrelation with the simulated amplitude-encoded features, validating the accuracy of our\noptical data preparation process.\nWe identify the primary cause of the accuracy gap between simulation (97.8%) and experiment\n(89%) can be attributed to fabrication imperfections, as has been commonly reported in state-\nof-art studies targeting classification tasks in the passive design regime3,6. Specifically,\nfabrication imperfections arising from lithographic in-plane etch bias52 can potentially alter the\nenergy distribution at the output ports. Moreover, alignment errors arising from non-\nuniformities in the fiber-to-chip VGC arrays may lead to undesired power variations between\nchannels. In our future work, the power discrepancies among the channels can be minimized\nthrough improved VGC uniformity51 and packaging method6. To more deeply examine the\ncause of this performance degradation, we simulated the effects of lithographic over and under\nexposure on the classification performance of our PNNs.\nOn Page 7 (Supplementary information):\nSupplementary Fig. 8 PNN robustness to fabrication variation. a Zoomed material index distribution of MedNIST chip\nwith overlaid contour for near uniform erosion (+10 nm, + 20 nm) and dilation (-10 nm, -20 nm) fabrication deviation around\nthe nominal design, approximating variation that may be experienced in electron beam lithography from over (under) exposure\nduring fabrication4. The simulated test accuracy (blue dots, left axis) and cross-entropy loss (red crosses, right axis) at each\nfabrication error condition is shown in b and c for MNIST and MedNIST datasets respectively\nNew references (Supplementary information):\n4. Christiansen, R. E. & Sigmund, O. Inverse design in photonics by topology optimization:\ntutorial. J. Opt. Soc. Am. B 38, 496–509 (2021).\nMoreover, our inverse-design framework is readily extensible to incorporate a variety of\nfabrication constraints to reduce the sensitivity of our device to experimental conditions, such\nas consideration of material loss via a custom defined complex material distribution as well as\nlithographic in-plane etch bias53 and out-of-plane thickness52 deviations to approximate the\neffect of process deviations. In addition, applying data augmentation6 during training can\npotentially improve the PNN’s robustness to experimental input permutations and help improve\ngeneralization.\nWe have updated the paper to provide further experimental details, a discussion on fabrication\nfeasibility and supplementary information quantifying the effect of common fabrication\ndeviations on device performance, along with new references.\nOn Pages 5-6:\nTo optimize the device material distribution across the scattering region, we directly perform\nthe optimization within the binarized regime, while incorporating fabrication process-specific\ndesign constraints every epoch (see Methods). This approach removes the necessity for\nhyperparameter scheduling and post-processing discretization35. As illustrated in Fig. 1c, the\nlatent design feature space, which varies continuously between 0 and 1, undergoes low-pass\nfiltering to smooth parameter spatial distributions and is subsequently binarized36 into\nfabrication-compatible materials: silicon (Si) and silicon dioxide (SiO ) as representative\n2\nexamples. A custom material distribution, characterized by its permittivity matrix 𝜖(𝜆) =\n𝜖 (𝜆)+𝑖𝜖 (𝜆), where 𝜖 (𝜆) and 𝜖 (𝜆) are its real and imaginary components respectively,\n+ , + ,\ncan model material loss and dispersive effects. We apply morphological operations to remove\nsmall features37, and B-spline contour approximations38-40 to further refine and shape the final\ngeometry under minimal feature size and radius of curvature constraints, ensuring\nmanufacturability41.\nOn Pages 15-16:\nThe classification score can be further improved by minimizing the power nonuniformities\namong the VGCs through enhanced fiber-to-chip alignment tolerance52. The deviation between\nsimulated and experimental results mainly arises from fabrication-induced geometry errors.\nWhen the fabrication variation (see Supplementary Fig. 8a) caused by lithographic over- or\nunder-exposure41,52 reaches approximately ±20 nm, changes in classification accuracy become\nevident for both MNIST and MedNIST devices (see Supplementary Fig. 8b, c). The tolerance\nto fabrication variations can be enhanced by independently simulating geometric perturbations\nand updating designs informed by these results53, while training with augmented datasets6 can\nfurther improve device robustness.\nOn Page 21:\nData preparation, calibration and phase robustness testing. Input light at 1550 nm from a\ntunable laser (Keysight N7778C) with a relative intensity noise less than -150 dB/Hz77 is fiber-\ncoupled into the chip.\nDownstream of MZIs, asymmetric 10:90 MMI are employed to tap a small portion of optical\npower for monitoring the attenuation in each channel, measured using multiport optical power\nmeters (Keysight N7744C) with dynamic range of -65 dB (exceeding 80 dB with range\nswitching)78 that incorporate built-in averaging to suppress random thermal and shot noises at the\nphotodiode readout79.\nNew references (main manuscript):\n52. Probst, M. J., Khurana, A, Slaby, J. B., Hammond, A. M. & Ralph, S. E. Fabrication\ntolerant multi-layer integrated photonic topology optimization. Opt. Express 32, 31448-\n31462 (2024).\n53. Wang, K. et al. Robust inverse design of digital photonic devices for photonic integrated\ncircuits. Opt. Express 33, 130-143 (2025).\n77. Keysight Technologies. Keysight N777-C family of tunable laser sources. Available at\nhttps://www.keysight.com/it/en/assets/3119-1067/data-sheets/5992-4217.pdf (2024).\n78. Keysight Technologies. N7744C and N7745C optical multiport power meters. Available\nat https://www.keysight.com/it/en/assets/3119-1066/data-sheets/5992-4218.pdf (2023).\n79. Keysight Technologies. Tips for optimizing power meter / sensor measurement speed.\nAvailable at https://www.keysight.com/it/en/assets/7018-03043/application-notes/5990-\n8471.pdf (2017).\nComments: In the same context, one key problem in PNN training using digital twins (offline\ntraining) is the fact that minor deviations between the physical chip and the model are\ntranslated to significant errors during training (designing here). In this case, where even\nnanometer location/size errors would significantly impact the scattering properties of the chip,\nimperfections between the physical system and its digital twin can result in major differences\nbetween expected performance (97.8%) and actual performance (89%). The authors should\ncomment on this in their revised manuscript.\nReply:\nWe thank the reviewer for their comment. We have embedded two fabrication-aware\nconstraints, minimal feature size and radius of curvature, into our design process. Methods to\nfurther reduce translation errors between simulation and experiment include applying\ngeometric perturbations at the length-scale of expected deviations53 and informing design\nupdates based on these field solutions or by applying data augmentation to make the device\ninherently more robust to input variations6 and improve generalization. Moreover, our inverse-\ndesign framework is compatible with fabrication-specific process design kits (PDKs)35, as a\nmethod of design parametrization (Fig. 1) for future large-scale applications.\nWe have updated the manuscript by providing new results to characterize the effect of\nfabrication deviations on the scattering properties of the PNN, and methods to improve its\nperformance.\nOn Pages 5-6:\nA custom material distribution, characterized by its permittivity matrix 𝜖(𝜆) = 𝜖 (𝜆)+𝑖𝜖 (𝜆),\n+ ,\nwhere 𝜖 (𝜆) and 𝜖 (𝜆) are its real and imaginary components respectively,\n+ ,\ncan model material loss and dispersive effects.\nOn Pages 15-16:\nThe classification score can be further improved by minimizing the power nonuniformities\namong the VGCs through enhanced fiber-to-chip alignment tolerance51. The deviation between\nsimulated and experimental results mainly arises from fabrication-induced geometry errors.\nWhen the fabrication variation (see Supplementary Fig. 8a) caused by lithographic over- or\nunder-exposure41,52 reaches approximately ±20 nm, changes in classification accuracy become\nevident for both MNIST and MedNIST devices (see Supplementary Fig. 8b, c). The tolerance\nto fabrication variations can be enhanced by independently simulating geometric perturbations\nand updating designs informed by these results53, while training with augmented datasets6 can\nfurther improve device robustness.\nNew references (main manuscript):\n52. Probst, M. J., Khurana, A, Slaby, J. B., Hammond, A. M. & Ralph, S. E. Fabrication\ntolerant multi-layer integrated photonic topology optimization. Opt. Express 32, 31448-\n31462 (2024).\n53. Wang, K. et al. Robust inverse design of digital photonic devices for photonic integrated\ncircuits. Opt. Express 33, 130-143 (2025).\nOn Page 7 (Supplementary information):\nSupplementary Fig. 8 PNN robustness to fabrication variation. a Zoomed material index distribution of MedNIST chip\nwith overlaid contour for near uniform erosion (+10 nm, + 20 nm) and dilation (-10 nm, -20 nm) fabrication deviation around\nthe nominal design, approximating variation that may be experienced in electron beam lithography from over (under) exposure\nduring fabrication4. The simulated test accuracy (blue dots, left axis) and cross-entropy loss (red crosses, right axis) at each\nfabrication error condition is shown in b and c for MNIST and MedNIST datasets respectively\nNew references (Supplementary information):\n4. Christiansen, R. E. & Sigmund, O. Inverse design in photonics by topology optimization:\ntutorial. J. Opt. Soc. Am. B 38, 496–509 (2021).\nComments: During experimental validation, the authors use an extremely limited number of\nsamples from the dataset (100 out of 60,000). What is the reason for this choice, and how can\nthis affect the statistical properties of the solution?\nReply:\nWe thank the reviewer for their comment. For our proof-of-concept experimental\ndemonstration, we selected a stratified random subset from the test set containing 10 samples\nfrom each class as our experimental measurement set. This provides an adequate performance\nbenchmark consistent in scale with the state-of-art passive PNN results3. To verify that this\nsubset provides sufficient representation of the test dataset, we quantified the statistical\nsimilarity between the test and experimental sets of data using three complementary analyses.\nWe first perform principal component analysis (PCA) of the first three components and show\nthere is significant overlap in the latent space for both datasets. Furthermore, we show that the\nmean average error per-feature distributions are low for MNIST (MAE = 0.0273) and\nMedNIST (MAE = 0.0336) datasets respectively, and the 1D Wasserstein distance49 for MNIST\n(𝑊 = 0.0212) and MedNIST (𝑊 = 0.0279) is small, which indicates strong alignment of the\nunderlying feature distributions. Finally, we confirm a low root mean square error (RMSE)\nprecision of ±2.02% and ±3.63% about the mean distribution via bootstrap resampling50\n(10,000 samples) of the experimental set for MNIST and MedNIST datasets respectively.\nThese results demonstrate both strong sample diversity and statistical robustness of the solution under\nresampling, validating that a limited, stratified subset provides adequate statistical\nrepresentation of the test set.\nWe have updated the main manuscript to clarify our experimental set acquisition and have\nprovided new supplementary figures to visualise both the sample diversity and statistical\nsimilarity between these two sets of data.\nOn Page 14:\nFor experimental data preparation, we measure analogue optical input power to the PNN\naccelerators across 100 MNIST and 60 MedNIST test samples, where 10 samples are randomly\nselected from each class (see Supplementary Fig. 4). To verify the accuracy of these inputs,\nwe measure a global mean average error (MAE) of 0.0277 for MNIST and 0.0310 for MedNIST\nwith respect to the simulated inputs (see Supplementary Fig. 5). These low MAE values\nconfirm a strong correlation with the simulated amplitude-encoded features, validating the\naccuracy of our optical data preparation process. To evaluate how well the experimental dataset\nrepresents the test distribution, we compute several distributional benchmarks. A 3D principal\ncomponent analysis shows that the experimental and test set features largely overlap within\ntheir 90% confidence ellipses (see Supplementary Fig. 6a, b). The mean feature distributions\nyield low per-feature MAE values of 0.0273 for MNIST and 0.0336 for MedNIST, and small\naverage 1D Wasserstein distances49 of 0.0212 and 0.0279, respectively, indicating high\ndistributional agreement (Supplementary Fig. 6c, d). Bootstrap resampling50 over 10,000\niterations (Supplementary Fig. 7a, b) produces global root-mean square error precisions of\n±2.02% for MNIST and ±3.63% for MedNIST relative to the test set mean at 90% confidence\n(Supplementary Fig. 7c, d). These results confirm that the stratified subset adequately\nrepresents the overall test distribution and that our experimental data remain statistically robust\nunder resampling.\nOn Page 21:\nFor both MNIST and MedNIST chips, a stratified random sample of 10 examples per class is\ntaken from the simulated test set to form the experimental validation set, demonstrating\ndiverse feature representation (see Supplementary Fig. 6).\nNew references (main manuscript):\n49. Villani, C. Optimal Transport: Old and New Vol. 338 (Springer, 2009).\n50. Efron, B. & Tibshirani, R. J. An Introduction to the Bootstrap (CRC Press, 1994).\nOn Page 6 (Supplementary information):\n15\n14\n13\n12\n11\n10\n9\n8\n7\n6\n5\n4\n3\n2\n1\n1 0.75 0.5 0.25 0 0.25 0.5 0.75 1\nTransmittance (a.u)\nSupplementary Fig. 6 Statistical representativeness of experimental PNN subset. 3D Principal component analysis (PCA)\nprojection of a MNIST and b MedNIST test dataset features, overlaid with the measured (experimental) set. Colored ellipses\nrepresent 90% covariance regions of each class, showing substantial overlap in the feature space. Bar plots with 90%\nconfidence intervals show mean transmittance with low MAE and Wasserstein (1D) distance2 (𝑊1) computed per-channel.\nThe globally averaged results for c MNIST (MAE = 0.0273, 𝑊1 = 0.0212) and d MedNIST (MAE = 0.0336, 𝑊1 = 0.0279)\nbetween the test (left-axis) and experimental (right-axis) distributions reinforce quantitative similarity.\nrebmuN\nlennahC\nChestCT CXR AbdomenCT Hand\nHeadCT BreastMRI Mean (Test set) Mean (Exp. set)\nMAE=0.034\nW=0.019\nW=0.027\nW=0.025\nMAE=0.034\nW=0.029\nW=0.075\nW=0.032\nW=0.032\nMAE=0.028\nW=0.019\nW=0.027\nW=0.023\nMAE=0.036\nW=0.028\nW=0.018\nMAE=0.030\nW=0.032\nW=0.016\n1 0.75 0.5 0.25 0 0.25 0.5 0.75 1\nTransmittance (a.u)\nrebmuN\nlennahC\nClass 0 Class 1 Class2 Class 3 Class 4\nClass 5 Class 6 Class7 Class 8 Class 9\nMean (Test set) Mean (Exp. set)\nTest Set Experimental Set Test Set Experimental Set\nMAE=0.015\nW=0.019\nMAE=0.035 W=0.018\nMAE=0.023\nW=0.020\nMAE=0.014\nW=0.024\nMAE=0.032\nW=0.022\nMAE=0.032\nW=0.016\nMAE=0.039\nW=0.017\nMAE=0.030\nW=0.033\nMAE=0.013\nW=0.023\nMAE=0.040 W=0.021\nOn Page 7 (Supplementary information):\n300\n200\n100\n0\n0.27 0.275 0.28 0.285 0.29 0.295 0.3 0.305\nRMSE\nSupplementary Fig. 7 Bootstrap analysis of experimental sample precision and stability. Bootstrap resampling3 (10,000\ncounts) of the experimental set are shown in the distribution plots for a MNIST and b MedNIST respectively, with 90%\nconfidence intervals (shaded) and standardized z-scores displayed for statistical reference. The convergence of precision\nestimates yield a full-run root mean square error (RMSE) precision (dashed line) of approximately ±2.02% and ±3.63% for c\nMNIST and d MedNIST datasets about the test-set mean. Together, these results indicate strong sampling stability of the\nmeasured experimental set.\nNew references (Supplementary information):\n2. Villani, C. Optimal Transport: Old and New Vol. 338 (Springer, 2009).\n3. Efron, B. & Tibshirani, R. J. An Introduction to the Bootstrap (CRC Press, 1994).\nReviewer #2\nComments: The results in this study is characterized by both physical modeling\n(characterization) and computational strategies. The authors employ three-dimensional\nfinite-difference time-domain (3D-FDTD) simulations as the backbone of their design\napproach. The methodology can reconstruct arbitrary optical field patterns through\ninterference, providing a foundation for the photonic neural network’s operation. Results\ncapture the complex wave dynamics within the nanophotonic circuits rather than relying on\nsimplified or purely theoretical models.\nA particularly innovative aspect is how the linearity of optical propagation is used to improve\nefficiency and scalability. Instead of simulating the entire neural network in one monolithic\n3D-FDTD run, the authors decouple the forward-pass computations into multiple linearly\nseparable simulations. In practice, this means different portions or layers of the photonic\nnetwork are simulated independently and can be run in parallel, taking advantage of GPU\ncomputing. Running many small FDTD simulations at the same time speeds up the design\nprocess and makes it easier to scale to larger systems. This advanced method uses the physics\nof light and powerful computing to design very compact optical devices.\nThe validity is reinforced by experimental verification. After designing the nanophotonic\nneural network in simulation, two prototype photonic chips on a silicon-on-insulator platform\nwere fabricated and tested on standard image-classification tasks. The correspondence\nbetween simulation and experiment is good: the devices (~20×20 μm² in size) achieved high\nclassification accuracies, approximately 89% on MNIST and 90% on MedNIST. This\nfabrication and testing step validate the modeling assumptions; the optical network behaved\nin practice as designed. The strength of this study is to build computer simulation and\nphysical implementation together. It provides a scalable and energy-efficient photonic\ncomputing platform in theory and confirms it in experiment\nReply:\nThank you for reviewing our manuscript and for your positive and encouraging feedback.\nBelow, we carefully address each of your detailed comments, providing clarifications,\nadditional analyses and new results, as well as expanded discussions in the revised manuscript.\nComments: Despite these strengths, several limitations and challenges should be mentioned.\nOne concern is the heavy computational demand of 3D-FDTD simulations. Achieving such a\ndetailed level of modeling accuracy requires considerable computing resources and time,\neven with parallelization. The authors reduce this through GPU acceleration and by splitting\nthe problem into smaller parts but scaling to truly large or deep neural networks may still\nface practical difficulties. As network size increases, the number of separate simulations and\nthe volume of data could increase dramatically, potentially offsetting the gains from\nparallelism. While the approach is mentioned scalable to large problem domains, its\nefficiency for much more complex networks remains to be demonstrated.\nReply:\nWe thank the reviewer for their comment. We propose a modularized network structure based\non the current design approach to achieve efficient network expansion in both depth (number\nof physical layers) and width (number of ‘neurons’ per layer) as requisite for more complex\nneural networks. The proposed method ensures the number of FDTD-solves scales linearly\nwith network depth and is decoupled from network width. Moreover, we propose methods for\nscaling our approach to datasets with large volume and number of classes, which are readily\nextensible by using our current inverse-design framework.\nWe have updated the paper to clarify this by adding a new section, “PNN model scalability to\nlarger neural networks”, which details the extensibility of our approach to larger network sizes\nand provide the training procedure and simulation cost analysis for scalable PNNs in the\nMethods section, with new results and additional supporting materials presented in the\nSupplementary Information.\nOn Pages 16-18\n:\nPNN model scalability to larger neural networks. We scale the model to larger neural\nnetworks by stacking PNNs with interleaved photodetector nonlinearity54, where each PNN is\ninverse-designed for task-optimized optical-field transformations. Figure 7 illustrates the\nschematic architecture for the PNN model scalable to larger neural networks. Input images are\nsegmented into patches, encoded as feature vectors, and processed by pairs of consecutive PNN\ncores. The first PNN maps input features to an optical-field distribution across class channels\n(Fig.7 Inset a), which are converted after photodetection into a series of patch class-activation\nvectors55. The second PNN re-encodes this vector into refined patch features (Fig.7 Inset b),\nwhich become increasingly aligned with class prototypes as they propagate in the next stage.\nEach pair is analogous to a layer block in multilayer DNNs56, where blocks are stacked to\nextract and refine feature representations57. Successive PNN blocks progressively refine these\nclass-level embeddings, filtering out patches with low confidence e.g. background, and\nsharpening those with high-class correlation, generating richer interpretability of feature-to-\nclass associations58. This depth-wise refinement expands the effective receptive field59, while\na sliding window with a fixed stride integrates broader spatial context. At the network output,\npatch-level predictions are confidence-weighted via monitored optical power and aggregated\nto yield the final image-level classification result (Fig.7 Inset c). Different output encoding\nschemes such as one-hot, binary or error-correcting output codes (ECOC)60 are viable for\nclassification where multi-class decisions are recovered by code-word decoding at inference54.\nWeight-sharing across network width enables patches to be evaluated in parallel61. Moreover,\noptical multiplexing62 such as wavelength or polarization, increases inference throughput and\nenables multi-task operation to be executed on the same chip (see Supplementary Section 2).\nBy using dual wavelengths, we demonstrate a 30 × 20 µm² single-chip design simultaneously\nclassifies MNIST and MedNIST datasets, achieving 95.1% and 98.0% test accuracy,\nrespectively (see Supplementary Section 2).\nA patch-efficient adjoint method (see Supplementary Section 3) is incorporated into the\ninverse-design of each PNN core, ensuring that the training cost is independent of the number\nof image patches. For a PNN core with input dimension 𝑁 and output dimension 𝐵, the number\nof FDTD solves per epoch is fixed at 𝑁+𝐵 under single-wavelength operation (see Methods).\nThis scaling keeps training tractable even for large-class datasets. For instance, with a 4 × 4\npatch kernel (𝑁 = 16) on ImageNet (1000 classes), where 𝐵 = 16 is determined by ECOC\ncoding54 with a 10-bit minimum code length plus 6 redundancy bits for error correction (see\nMethods), each PNN core design update requires only 32 FDTD-solves.\nFig. 7 Schematic diagram of PNN model scalable to larger neural networks. Images are segmented into\npatches (1…𝑘…𝑤) by a sliding window and encoded as feature vectors. Patches are processed by a sequence of\ninverse-designed PNN cores [PNN(\"),PNN($),… , PNN(%), PNN(%&\"),… , PNN(')] with interleaved\nphotodetector nonlinearities (NL). Width scalability is achieved by weight sharing across patches, while depth\nscalability arises from stacking cores with re-extraction of overlapping patches. Within each layer block, the first\nPNN core maps a feature vector to an optical-field distribution |𝐄| over 𝐵 channels for class evidence (Inset a)\nand forms the patch class-activation vector after photodetection. This vector is re-encoded onto the optical input\nports |𝐄′| of the second PNN (Inset b). These inputs are used to refine and embed class-prototype information in\npatch features before propagating to the next PNN core. Repeating this process sharpens per-patch class\nattributions, which are aggregated by confidence-weighted pooling into image-level predictions (Inset c). The\nsame forward operations occur during training. In addition, pooling is applied at intermediate layers for auxiliary\nsupervision.\nNew references (main manuscript):\n54. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n55. Brendel, W. & Bethge, M. Approximating CNNs with bag-of-local-features models\nworks surprisingly well on ImageNet. In Proc. 7th International Conference on Learning\nRepresentations (ICLR, 2019).\n56. Alzubaidi, L. et al. Review of deep learning: concepts, CNN architectures, challenges,\napplications, future directions. J. Big Data 8, 53 (2021).\n57. Hou, L. et al. Patch-based convolutional neural network for whole slide tissue image\nclassification. In Proc. IEEE Conference on Computer Vision and Pattern Recognition\n2424–2433 (IEEE, 2016).\n58. Choromanski, K. et al. Rethinking attention with performers. In Proc. 9th International\nConference on Learning Representations (ICLR, 2021).\n59. Luo, W., Li, Y., Urtasun, R. & Zemel, R. Understanding the effective receptive field in\ndeep convolutional neural networks. In Proc. 30th Conference on Neural Information\nProcessing Systems, 4905–4913 (NIPS, 2016).\nNL\nNetwork Depth\nSliding window\nPatch kernel\n|1| distribution\nClass 1 Class 2\nClass C\n|k| distribution |k|' distribution\nClass 1 Class 2\nClass C\nInput\nFeatures\n|w| distribution |w|' distribution\nRefined\nFeatures\ntupnI PNN\n(1)\nnoitacfiissalC\ntuptuO\nClass Evidence Feature Embedding\nClass 1 Class 2\nClass C\nPatch confidence\npooling\nImage-level\nclassification\nPatch\nPatch\nw\nClass\nFeature\nNL\nNetwork Width\nClass Decision\n|1|' distribution\nClass 1 Class 2\nClass C\nClass 1 Class 1\nClass 2 Class 2\nClass C Class C\nPNN\n(2)\nNL\nPNN\n(i)\nNL\nPNN\n(i+1)\nNL\nPNN\n(l)\n1\nk\nw\n1 2 C\nClass\nFeature\nClass\nFeature\npamtaeH\nssalC\nLayer Block\n60. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\n61. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n62. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nWe also discuss methods to further scale the PNN structure via multiplexing and demonstrate\nthe extensibility of our approach to the multi-task domain by designing a single 30 × 20 𝜇m2\nchip which performs both MNIST and MedNIST classification, multiplexed by wavelength.\nOn Page 17:\nMoreover, optical multiplexing62 such as wavelength or polarization, increases inference\nthroughput and enables multi-task operation to be executed on the same chip (see\nSupplementary Section 2). By using dual wavelengths, we demonstrate a 30 × 20 µm²\nsingle-chip design simultaneously classifies MNIST and MedNIST datasets, achieving 95.1%\nand 98.0% test accuracy, respectively (see Supplementary Section 2).\nNew references (main manuscript):\n62. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nOn Pages 8-9 (Supplementary information):\nSupplementary Section 2. Multiplexing method\nWeight sharing in the form of hardware/time multiplexing can be leveraged to scale the network\nin width13,14. Moreover, wavelength multiplexing15,16 can be incorporated into our inverse\ndesign process to enable multiple optical transformations of different dimensions to be realized\nwithin a single PNN. To demonstrate this flexibility, we design a single 30 × 20 𝜇m2 chip\ncapable of performing both MNIST and MedNIST classifications via wavelength multiplexing,\nachieving simulated test accuracy of 95.1% and 98.0%, respectively.\nAt each inverse-design update, a fixed number 𝑁′ = max(𝑁 ,𝑁 ) of independent\n! !\"#\nsimulations are run to compute the training forward-pass fields 𝐄 as well as 𝐶′ = max(𝐶 ,\n$%# !\n𝐶 ) independent reverse-mode simulations to compute the back-propagation fields, where\n!\"#\n𝑁 ,𝑁 are the input feature dimensions and 𝐶 , 𝐶 are the output class logit dimensions\n! !\"# ! !\"#\nfor MNIST and MedNIST datasets respectively. Computing the gradient fields requires a single\nmultiplication of the forward 𝐄 and adjoint 𝐄 fields per sample of the training dataset.\n$%# &#'\nFor MNIST, the 𝑁 = 10 input features are distributed across all input ports, with zero-padding\n!\napplied to any unused ports. Likewise, for MedNIST, the class outputs occupy the central six\nports, as shown in Supplementary Fig. 9. Supplementary Fig. 9a, b shows the field magnitude\n|𝐄| for random training samples in classes ‘Hand’ and ‘1’ for MedNIST and MNIST datasets\nmultiplexed at wavelengths 1570 nm and 1530 nm respectively, with the field hot-spot\nconcentration at the correct port. The confusion matrix in Supplementary Fig. 9c, d denotes\na test accuracy of 98.0% and 95.1% for MedNIST and MNIST datasets respectively,\ndemonstrating simultaneous classification performance for both datasets.\n'0'\n'1'\n'ChestCT' '2'\n'CXR' '3'\n'AbdomenCT' '4'\n'Hand' '5'\n'HeadCT' '6'\n'BreastMRI' '7'\n'8'\n'9'\nSupplementary Fig. 9 PNN multi-task classification. Field magnitude |𝐄| results for a sample from MedNIST\na and MNIST b datasets are presented, with a hotspot at the correct classification port corresponding to classes\n‘Hand’ at a wavelength of 1570 nm and ‘1’ at a wavelength of 1530 nm respectively. The confusion matrices in c\nand d shows a test accuracy of 98.0% and 95.1% for MedNIST and MNIST classes respectively, demonstrating\nstrong classification performance.\nNew references (Supplementary information):\n13. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n14. Wang, X. et al. Integrated photonic encoder for low power and high-speed image\nprocessing. Nat. Commun. 15, 4510 (2024).\n15. Cheung, A. K. C. et al. Inverse-designed CWDM demultiplexer operated in O-band.\nIn Optical Fiber Communications Conference, W1A-6 (IEEE, 2024).\n16. Zhang, Y., Zhang, K., Hu, P., Li, D. & Feng, S. Multi-wavelength diffractive optical\nneural network integrated with 2D photonic crystals for joint optical\nclassification. Nanophotonics 14, 2891-2899 (2025).\nComments: Two, it was unclear on how these models can be scaled to larger neural networks\nwith several layers. How will the training occur? Again with simulation? It was not clear on\nhow to apply the inverse design to large CNN and DNN models.\nReply:\nWe thank the reviewer for their comment. We scale these models by modularizing PNN\nmodules and stacking them sequentially to extract patch-wise class evidence, which is then\nused to embed class prototype features into the global image information for interpretable\nclassification55. Each PNN core is optimized sequentially using our 3D-FDTD inverse-design\nmethod and reusing earlier-trained blocks as priors can improve training efficiency by\naccelerating optimization in deeper layers78. The simulation cost of these models is linear in\nnetwork depth and decoupled from network width.\nWe have updated the paper to clarify how these PNN models can be scaled by adding a new\nMethods section “Training framework and scalability”, which details the training method\nextended to larger, multi-layer networks, including a patch-efficient adjoint method and\nestimates of the computational cost detailed in the Supplementary section “Patch-efficient\nadjoint gradient method”, along with supportive references.\nOn Pages 21-22:\nTraining framework and scalability. Each PNN core is trained to implement task-specific\noptical-field transformations by minimizing a CE loss1. Class information is embedded through\none-hot, binary or ECOC coding60, with ECOC advantageous for large-class datasets by\nreducing the number of output channels required54. A patch-efficient adjoint method (see\nSupplementary Section 3) is incorporated into the optimization, whereby patch contributions\nare aggregated algebraically before gradient computation. This ensures that the total number\nof required FDTD solves per epoch for one per PNN core depends only on the input and output\ndimensions (𝑁 + 𝐵) and is independent of the number of image patches (see Supplementary\nSection 3). The optimization incorporates nonlinear photodetection, which converts optical\nlogits into class-activation vectors and provides interpretable measures of patch-to-class\nassociation55. To stabilize training and ensure global consistency, auxiliary supervision can be\nintroduced by aggregating patch-level predictions with confidence weighting at intermediate\nlayers, enforcing agreement with image-level classification. Training is proceeded block by\nblock, where each consecutive pair of PNN cores are optimized jointly and then frozen after\nconvergence. Overlapping patches are re-extracted and passed to the next stage, enabling depth\nscalability through successive refinement. Reuse of pre-trained shallow blocks as priors can\nfurther improve efficiency by reducing the number of gradient updates needed for deeper\nlayers81. The final PNN core maps refined patch features into class-activations, and confidence-\nweighted pooling across patches produces the image-level prediction. All optimizations are\nperformed in simulation, and the resulting geometries are used directly for inference in the\nfabricated device.\nNew references (main manuscript):\n54. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n55. Brendel, W. & Bethge, M. Approximating CNNs with bag-of-local-features models\nworks surprisingly well on ImageNet. In Proc. 7th International Conference on Learning\nRepresentations (ICLR, 2019).\n60. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\n81. Chen, T., Goodfellow, I. & Shlens, J. Net2net: accelerating learning via knowledge\ntransfer. In Proc. International Conference on Learning Representations (ICLR, 2016).\nOn Page 10 (Supplementary information):\nSupplementary Section 3. Patch-efficient adjoint gradient method\nThis section provides the adjoint gradient method reformulated for patch-efficiency to\nexplicitly show that the simulation cost is independent of the number of image patches. For\neach PNN core, only a fixed number of independent FDTD simulations are required per epoch:\n𝑁 forward simulations corresponding to the flattened, input patch dimension, and 𝐵 adjoint\nsimulations corresponding to the output PNN dimension. All remaining patch-level gradients\nare resolved algebraically through linear superposition in the time-harmonic domain, without\nperforming additional FDTD solves.\nLet 𝑥 denote the input feature vector of patch 𝑡, and 𝛼 the corresponding output class-\n- -\nevidence vector. The outer product of these quantities is accumulated across all patches to form\na transform matrix 𝐻 given by\n𝐻 = @𝑥 𝛼. ∈ ℝ/×1 (1)\n- -\n-\nwhere ℝ denotes the set of real numbers.\nForward and adjoint field bases are then simulated once for each dimension, represented as\n𝐄 (𝒓) for 𝑖 ∈ {1,…,𝑁} and 𝝍 (𝒓) for 𝑏 ∈ {1,…,𝐵} respectively, where 𝒓 = {𝑥,𝑦,𝑧} is the\n, 2\nposition vector in 3-D.\nOwing to the linearity of Maxwell’s equations, the overlap 𝑀 between forward and adjoint\nfields can be expressed as\n𝑀 (𝒓) = ℛ{𝐄 (𝒓)⊙𝝍 (𝒓)} (2)\n,,2 , 2\nwhere ℛ{⋅} denotes the real-part operator and ⊙ denotes the Hadamard product. This forms a\nprecomputed tensor of fixed dimension [𝑁×𝐵,𝑛 ,𝑛 ,𝑛 ], where 𝑛 ,𝑛 and 𝑛 denote the\n4 5 6 4 5 6\nfield matrix dimensions in the 𝑥, 𝑦 and 𝑧-axes respectively.\n/ 1\n𝜕ℒ\n(𝒓) = @@𝐻 ×𝑀 (𝒓) (3)\n𝜕𝜖 ,,2 ,,2\n,78278\nwhere 𝐻 denotes the (𝑖,𝑏)-th element of 𝐻.\n,,2\nTherefore, through incorporation of the patch-\nefficient adjoint gradient method, arbitrarily large images require simulating a fixed (𝑁+𝐵)\nFDTD fields and algebraically resolving a constant (𝑁×𝐵) superposition-fields via GPU-\naccelerated matrix multiplication. In the one-hot case, 𝐵 = 𝐶, where 𝐶 is the number of output\nclasses. For the ECOC17 case, 𝐵 = ⌈log 𝐶⌉ + 𝑘, where ⌈⋅⌉ denotes the ceiling function and 𝑘 is\n)\nthe number of redundancy bits added for error tolerance13. Consequently, the total number of\nsimulations per epoch per PNN core is fixed at 𝑁+𝐵, decoupled from both the number of\npatches and the dataset size.\nNew references (Supplementary information):\n13. Xu, Z. et al. Large-scale photonic chiplet Taichi empowers 160-TOPS/W artificial\ngeneral intelligence. Science 384, 202-209 (2024).\n17. Dietterich, T. G. & Bakiri, G. Solving multiclass learning problems via error-correcting\noutput codes. J. Artif. Intell. Res. 2, 263–286 (1995).\nComments: Third, is the need for careful incorporation of real-world constraints in the\nsimulation. Factors like fabrication tolerances and material losses should be accounted in\ndesign to ensure the simulated performance translates into actual devices. If the current\nsimulations assume ideal conditions, that could be a limitation. Slight manufacturing\nimperfections, for example, will need to be solved for larger-scale or commercial\napplications. It appears to be an ideal case rather than any fabrication imperfections.\nReply:\nWe thank the reviewer for their comment. We have embedded two fabrication-aware\nconstraints, minimal feature size and radius of curvature, into our design process. Our inverse-\ndesign framework is readily extensible to incorporating other constraints, such as consideration\nof material loss, lithographic in-plane etch bias and out-of-plane thickness deviations as ways\nto approximate process imperfections. Applying data augmentation during training is also\nexpected to improve the PNN’s robustness to experimental input permutations and help\nimprove generalization6. Moreover, our inverse-design framework is compatible with existing\nfabrication-specific process design kits (PDKs)36 as a method of design parametrization (Fig.\n1) for future large-scale applications.\nWe have updated the manuscript by providing new results to characterize the effect of common\nprocess deviations on the scattering properties of each PNN, and methods to improve continuity\nbetween the physical system and its digital twin.\nOn Pages 5-6:\nTo optimize the device material distribution across the scattering region, we directly perform\nthe optimization within the binarized regime, while incorporating fabrication process-specific\ndesign constraints every epoch (see Methods). This approach removes the necessity for\nhyperparameter scheduling and post-processing discretization35. As illustrated in Fig. 1c, the\nlatent design feature space, which varies continuously between 0 and 1, undergoes low-pass\nfiltering to smooth parameter spatial distributions and is subsequently binarized36 into\nfabrication-compatible materials: silicon (Si) and silicon dioxide (SiO ) as representative\n2\nexamples. A custom material distribution, characterized by its permittivity matrix 𝜖(𝜆) =\n𝜖 (𝜆)+𝑖𝜖 (𝜆), where 𝜖 (𝜆) and 𝜖 (𝜆) are its real and imaginary components respectively,\n+ , + ,\ncan model material loss and dispersive effects. We apply morphological operations to remove\nsmall features37, and B-spline contour approximations38-40 to further refine and shape the final\ngeometry under minimal feature size and radius of curvature constraints, ensuring\nmanufacturability41.\nOn Pages 15-16:\nThe deviation between simulated and experimental results mainly arises from fabrication-\ninduced geometry errors. When the fabrication variation (see Supplementary Fig. 8a) caused\nby lithographic over- or under-exposure41,52 reaches approximately ±20 nm, changes in\nclassification accuracy become evident for both MNIST and MedNIST devices (see\nSupplementary Fig. 8b, c). The tolerance to fabrication variations can be enhanced by\nindependently simulating geometric perturbations and updating designs informed by these\nresults53, while training with augmented datasets6 can further improve device robustness.\nNew references (main manuscript):\n52. Probst, M. J., Khurana, A, Slaby, J. B., Hammond, A. M. & Ralph, S. E. Fabrication\ntolerant multi-layer integrated photonic topology optimization. Opt. Express 32, 31448-\n31462 (2024).\n53. Wang, K. et al. Robust inverse design of digital photonic devices for photonic integrated\ncircuits. Opt. Express 33, 130-143 (2025).\nOn Page 7 (Supplementary information):\nSupplementary Fig. 8 PNN robustness to fabrication variation. a Zoomed material index distribution of MedNIST chip\nwith overlaid contour for near uniform erosion (+10 nm, + 20 nm) and dilation (-10 nm, -20 nm) fabrication deviation around\nthe nominal design, approximating variation that may be experienced in electron beam lithography from over (under) exposure\nduring fabrication4. The simulated test accuracy (blue dots, left axis) and cross-entropy loss (red crosses, right axis) at each\nfabrication error condition is shown in b and c for MNIST and MedNIST datasets respectively\nNew references (Supplementary information):\n4. Christiansen, R. E. & Sigmund, O. Inverse design in photonics by topology optimization:\ntutorial. J. Opt. Soc. Am. B 38, 496–509 (2021).\nResponse to Reviewers’ Comments\nWe thank the reviewers for their comments and recommendations. Our response to each\nreviewer is given below.\nReviewer #1\nComments:\nThe authors addressed all the points, clarifying all the obscure points and significantly\nincreasing the overall quality of the manuscript. I believe that the manuscript is fit for\npublication.\nReply: Thank you for your positive feedback and recommendations which helped us improve\nthe quality of the manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors addressed all the points, clarifying all the obscure points and significantly\nincreasing the overall quality of the manuscript. I believe that the manuscript is fit for\npublication."}, "author_response": "Reply: Thank you for your positive feedback and recommendations which helped us improve\nthe quality of the manuscript."}]

Decision letter