Ultrafast neural sampling with spiking nanolasers

Domain: physics
Published: 2025-12-03
Editorial outcome: major_revision

Abstract
Owing to their significant advantages in terms of bandwidth, power efficiency, and latency, optical neuromorphic systems have arisen as interesting alternatives to digital electronic devices. Recently, photonic crystal nanolasers with excitable behavior were first demonstrated. Depending on the pumping strength, they emit short optical pulses – spikes – at various intervals on a nanosecond timescale. In this theoretical work, we show how networks of such photonic spiking neurons can be used for Bayesian inference through sampling from learned probability distributions. We provide a detailed derivation of translation rules from conventional sampling networks, such as Boltzmann machines, to photonic spiking networks and demonstrate their functionality across a range of generative tasks. Finally, we provide estimates of processing speed and power consumption, for which we expect improvements of several orders of magnitude over current state-of-the-art neuromorphic systems.

Introduction
In the overwhelming majority of artificial neural networks (ANNs) in use today, neuronal outputs are determined by smooth and continuous activation functions, and their values are updated synchronously for multiple neurons within a network. In contrast, biological neurons, especially in the mammalian cortex, communicate asynchronously with short, stereotypical pulses called spikes. While the non-differentiability of such spike-based codes indeed makes them more difficult to train using error backpropagation and therefore less attractive for conventional deep learning (however, see refs. 1–3), they possess several properties that are extremely relevant in physical neuronal networks, biological and artificial alike. First, they can be much more efficient than time-continuous or rate codes: as information becomes implicitly encoded in the timing of the spikes, rather than the value of neuronal outputs, a spike code can save both energy and bandwidth. Thus, spiking neurons have become the de facto standard model for neuromorphic systems, both purely digital 4–6 and mixed-signal 7–10. Second, due to the all-or-nothing nature of spikes, the transmission of information is significantly less prone to disruption by noise. This makes spiking neural networks particularly attractive for mixed-signal devices, whose analog components (neurons and synapses) invariably introduce spatio-temporal noise into the network dynamics.

By associating neuronal spikes with binary states, one can link spiking neural networks to Boltzmann machines (BMs) 11–13 (but see also ref. 14 for an alternative implementation). This machine learning model forms the basis for powerful variants such as deep belief networks 15 and autoencoders 16, which have been used for applications ranging from acoustic modeling 17 to medical diagnostics 18 and quantum tomography 19. By moving to the spiking domain, one can harness the speed and efficiency of neuromorphic substrates for such applications. Indeed, first demonstrations of such spike-based sampling networks on highly accelerated neuromorphic hardware have already shown promising results 20–22.

However, despite the impressive results of electric and electronic networks, electrical interconnects between computing elements limit bandwidth, latency, and energy efficiency 23. In contrast, photonic transmission of information offers significant advantages in terms of loss, bandwidth, and computing speed. Not only is the issue of Joule heating alleviated, but the improved fan-in/-out by wavelength multiplexing also facilitates large-scale connectivity. The initial motivation of improving interconnections between digital electronic computing elements has recently extended to actual photonic computing, in particular, implementing matrix-vector multiplication. Two studies have recently and independently introduced a 64 × 64 photonic matrix demonstrating exceptional computing speed and latencies 24,25. Current research on photonic neural networks implements artificial neurons similar to those in machine learning, which are deterministic and have a continuous activation function 26.

In this work, a photonic equivalent of a neuron is a laser with an excitable response 27–33; interaction between such neurons with highly accelerated intrinsic dynamics is being developed 34. Among a variety of possible implementations, a semiconductor laser with a saturable absorber (SA) is of particular interest here 35 because of the strong analogy with the leaky integrate-and-fire (LIF) neuron model 36. Semiconductor lasers can be miniaturized and integrated in a photonic circuit and need tiny amounts of power (about 10 μW) 37, while, by trading higher levels of current (100 μW), nearly 14% wall-plug efficiency directly into a silicon photonic circuit has been demonstrated 38. These are essential requirements for advancing the neuromorphic state of the art. Very recently, this technology was shown to produce an experimentally observed excitable response by adding a resonator section acting as an SA 39, thus paving the way for efficient integrated photonic spiking neurons (PSNs). The combination of a pulsed, stochastic response followed by a pronounced refractory period allows the implementation of neural sampling as a means of approximating the Glauber dynamics found in BMs. In contrast to most photonic neural networks found in literature, such networks operate on probability distributions rather than signals, which fills an important niche of processing non-deterministic and partial inputs. In this work, we demonstrate the capability of such PSN networks to implement BMs through spike-based sampling.

Results
Preliminaries
A BM is a stochastic neural network composed of binary neurons connected by a symmetric weight matrix. All possible states of the network can be described with a set of vectors z, where z k ∈ {0, 1}, i.e., not active or active. The probability of the BM being in a particular state z is p(z) ∝ exp(−E(z) / k B T), where k B is the Boltzmann constant, T is the ensemble (Boltzmann) temperature, and E(z) is the energy of this state (Fig. 1a).

E(z) = − ∑ k<j W kj z k z j − ∑ k b k z k, (1)

where b k are the neuronal biases, and W is an asymmetric coupling matrix that represents synaptic weights. Since the influence of the temperature k B T can be absorbed into the weights and biases, we can disregard the units, assume k B T = 1, and omit this term for brevity from here on.

In general, a BM can be used for Bayesian inference in previously learned probability spaces. First, biases and weights are adapted such that p(z) follows a particular distribution determined by the training data. Then, an input (evidence or observation) is provided by changing biases, and the output (the result of the inference problem) is simply obtained by measuring the resulting changes to p(z). For example, one can train a BM to follow a distribution over natural images and their labels. Afterwards, when provided with a partial image (by bias-clamping the neurons representing the corresponding pixels to their respective values), the BM will generate states of the unobserved neurons that correspond to plausible completions of the partial image, along with the corresponding labels.

However, the state space of N binary neurons has a size of 2 N, which renders the computation of the probability distribution intractable even for modest networks. Sampling can be a viable alternative: a histogram of samples approximates the true probability distribution, and the accuracy increases with more samples drawn. Fast sampling is therefore preferable, for which the nature of the hardware implementation matters a lot. Here, we consider an implementation with spiking neurons, as their dynamics lend themselves to ultrafast and energy-efficient implementation with nanolasers, as discussed further below.

To start, consider a time-continuous counterpart to a classical BM neuron. Whenever a neuron fires a spike, it enters a refractory period of some duration τ during which it cannot spike again. Thus, a spike can be interpreted as representing the onset of an active state

z k(t) = 1 ⇔ neuron has fired in (t − τ; t], (2)

with z k(t) = 0 otherwise. Conversely, a network of N spiking neurons can then serve to represent the full state z(t) ∈ {0, 1} N.

In spike-based sampling, a spiking neural network approximates p(z) through a time-continuous analogon of Gibbs sampling 11, a sampling algorithm where components z i of a current sample z(n) are updated iteratively using the conditional probability distribution: z i(n+1) ∼ p(z i | z k(n) ∀ k ≠ i). A particularly interesting variant builds on LIF sampling (LIFS) neurons 12,13, as they represent a de-facto standard model across the vast majority of neuromorphic platforms. In LIFS, the required stochasticity is generated by adding (or exploiting pre-existing) noise on the neuronal membranes, which consequently follow the Ornstein-Uhlenbeck process

τ m du j = (b j − u j) dt + ∑ k W kj κ(t − t spike k) dt + σ W dW t (3)

with an autocorrelation determined by the neuronal membrane time constant τ m (dW represents a Wiener process scaled by a noise amplitude σ W). The interaction between neurons is mediated by an additive postsynaptic potential (PSP) kernel κ that is triggered by an incoming spike at time t spike k.

Spiking nanolasers
The PSNs discussed in this work are semiconductor lasers composed of two sections: a gain section and an SA. The two sections form a single resonator and are therefore spanned by a single laser mode (Fig. 1a). The gain section is pumped to reach local electron population inversion and stimulated emission of photons. The SA is not pumped, such that here absorption always dominates, yet saturates as the flux of absorbed photons increases. As the pump rate increases beyond a certain threshold, the absorption is overtaken by the stimulated emission; akin to negative differential resistance in electronic oscillators, this introduces a positive feedback which initiates the emission of a pulse. In turn, this depletes the population of excited electron-hole pairs, hence the gain, and the laser is shut off. Therefore, the emission of a new pulse immediately following a previous one is strongly suppressed. After some time, the pumping re-establishes the initial gain, and the nanolaser can spike again. This process (Fig. 1b) bears a strong resemblance to the generation of action potentials and the subsequent refractoriness found in biological neurons, as described by the Hodgkin-Huxley model 40.

The spiking dynamics in a semiconductor laser with an SA is described by the Yamada model 35,41, which ignores spontaneous emission and therefore assumes a perfectly deterministic response. Here, we consider a laser with a single optical mode; its field is distributed over a volume comparable to λ 3, where λ is the wavelength. Moreover, the active region, where electron and hole pairs are created, is even smaller. Under these conditions, spontaneous emission is not negligible, as the fraction of it going into the mode – the spontaneous emission factor – is between 0.1 and 1.0, whereas in macroscopic semiconductor lasers it is 10 −4 or less 42. Therefore, in nanolasers, the average number of photons is much smaller, and the relative noise due to the granularity is much larger, which disrupts the otherwise deterministic spiking at regular intervals controlled by pumping strength.

A rigorous description of noise requires a quantum mechanical formalism, which is exceedingly complicated for semiconductor systems. Therefore, the semiconductor laser is often approximated as a homogeneously broadened two-level system, as intraband scattering is fast enough such that electrons and holes are in thermal equilibrium 43,44. The light-matter interaction within a semiconductor can therefore be described as a collection of dipoles interacting with the same optical mode. With some approximations, such as fast dephasing of the polarization with respect to the damping rate of carriers and photons, the quantum mechanical descriptions lead to rate equations 45,46.

The rate equations describe the evolution in time of the populations of excited dipoles n e, and photons S in a cavity mode due to multiple processes. Stimulated emission and absorption are described by a term γ r(2n e − n 0), where n 0 is the total number of dipoles, and γ r is the radiative transition rate. Since n e ≤ n 0, this term is limited to γ r n 0; this condition is ensured by an optical pumping term γ p(n 0 − n e), where γ p is the pumping rate and n 0 − n e represents pumping saturation. Spontaneous emission and photon damping are described by γ r n e and γ t S. The stochastic equations are formed by including the Langevin forces F i(t), which are random variables with zero mean and auto-/cross-correlation strengths 〈F i(t) F j(t)〉 = 2D ij. The Langevin forces are calculated consistently with the rate equations 46 based on the McCumber noise model 47 (see Sec. “Nanolaser noise model” in the Methods). Here, we extend the model by including two separate sections, one providing amplification with an excited population n e, and another one representing the SA, denoted with the suffix a, with an excited population n a, both interacting with the same mode:

Ṡ = GS + γ r n e + γ r,a n a + F S(t),
ṅ e = −γ r(2n e − n 0)S − γ t n e + γ p(n 0 − n e) + F e(t),
ṅ a = −γ r,a(2n a − n 0,a)S − γ t,a n a + F a(t), (4)

where G = γ r(2n e − n 0) + γ r,a(2n a − n 0,a) − γ t is the gain including photon damping, and other parameters are defined in Supplementary Table 1. For clarity, we normalize γ p by the threshold pumping strength in the absence of noise γ p thr (see Sec. “Nanolaser with quantum wells” in the Methods) which we later refer to as the threshold. It is important to note that with noise, spike emission can also occur with γ p < γ p thr.

Figure 2a shows a time trace of the membrane potential of a single LIFS neuron with parameters from ref. 13 and Gaussian noise. This serves as a reference for the behavior of PSNs, which we discuss below. In Fig. 2b, we show time traces of PSN gain with pumping far from and close to the spiking threshold, respectively. In this system, the gain is a stochastic variable as it depends on stochastic electron densities. As a result, when a PSN is not refractory, the gain undergoes a random walk as shown in Fig. 2c, similarly to the LIFS membrane potential. The inter-spike interval distribution, gain autocorrelation, and probability density function of a PSN are shown in Fig. 2d, e, and f, respectively. These are close to the corresponding properties of LIFS neurons, with small deviations explained by the additional nonlinear terms of the PSN state equations.

It has been pointed out that the granularity of light challenges numerical integration based on the Langevin forces 48–50. Therefore, we further cross-checked the results by implementing a rigorous discrete PSN model, as proposed in ref. 49 and concluded that both methods lead to the same statistical properties of the PSN (see Sec. “Discrete nanolaser model” in the Methods).

Networks of spiking nanolasers
While many different implementations of our proposed network can be conceived, we discuss several ideas in the following. Interference-based interaction between PSNs is challenging: their bistability, combined with a potentially large number, will complicate the necessary resonance alignment. For this reason, in this work we assume that PSNs are incoherent, and their interaction is mediated by photodiodes. Spikes are extracted from a PSN by coupling to a waveguide. All-to-all connection weights can be implemented using a photonic matrix (see Fig. 3c), a photonic integrated circuit developed to implement matrix-vector multiplication 24,25; the use of this technology here assumes incoherent superposition, as nanolasers may emit at different wavelengths, yet this still leaves a margin for setting the weights to their desired values. Wavelength diversity also allows the combined use of optical crossbars (see Fig. 3d) based on wavelength division multiplexing 51. In this respect, potential limitations due to the intentionally reduced interference and to the fact that the detector measures a positive quantity may need to be managed. In fact, in BMs, weights can be negative as well. A possible solution is the use of balanced photodetectors 52, which requires an N × 2N coupling matrix for N PSNs.

As a result, optical spikes are converted into an electrical current that changes the pumping strength by

Δγ p,k(t) = ∑ j κ kj ∫ LPF(t − t*) S j(t*) dt*, (5)

where κ kj represents a coupling strength incorporating relevant factors such as the coupling of a PSN to a waveguide and the photodiode responsitivity, and LPF(Δt) is the impulse response function of the photodiodes that we assume to be a low-pass filter with a timescale τ U due to their limited bandwidth. The change of current induces a change of gain shown in Fig. 3b; the change can be positive (excitatory) or negative (inhibitory), as shown in Fig. 3e. In the next section, we propose that the gain plays the role of the membrane potential in a PSN. Therefore, Eq. (5) describes the PSP; its shape is largely defined by the bandwidth of the electrical scheme (Fig. 3a, b).

This implementation of connections induces an almost alpha-shaped interaction kernel t/τ exp(−t/τ), again similar to typical interaction kernels in both biological and abstract (LIF) neuronal networks. For optimal performance of a network, the timescale of Δγ p,k(t) can require optimization. Here, it depends on τ U, which we assume to be tunable, as discussed further below.

Membrane potential of spiking nanolasers
Consider a network of nanolasers and its k-th PSN is not currently spiking, i.e., its gain does a random walk (see Fig. 2c). Then, Eq. (4) is approximated by the following stochastic equation for the gain (see Sec. “Simplified gain equation” in the Methods):

dG k = −(G k − G p,k) / τ G dt + 2γ r n 0 Δγ p,k dt + σ G dW, (6)

where G p is a drift term depending on the pumping strength, and the second term corresponds to the PSN interaction (see Eq. (5)). This equation is identical to that of the membrane potential of an LIFS neuron (see Eq. (3)) with two assumptions. First, n e ≪ n 0, which holds during non-spiking with moderate pumping strength. Second, changes of n e due to incoming spikes need to be small, otherwise the interaction between PSNs deviates from linearity (see Sec. “Simplified gain equation” in the Methods). Based on the formal equivalence between the two equations, we propose to use the gain G as the membrane potential of a PSN:

u k(t) = (G k(t) − G 0) = ∂ u G, (7)

where G 0 is the gain for which p(z = 1) = 0.5 (which requires the yet unknown refractory period τ), and ∂ u G is a proportionality coefficient. Therefore, we must to find these three parameters to relate the membrane potentials of LIFS neurons and PSNs.

Consider a single LIFS neuron with a logistic activation function p(z = 1 | u) = σ(u). We expect a similar behavior from a single PSN, where G = G p (see Eq. (6)). In Fig. 4a we sweep the pumping strength G p; for each G p, we solve Eq. (4) for a sufficiently long time and estimate the activation function using Eq. (2) and the law of large numbers:

p(z = 1 | G p) ≈ z̄(G p), (8)

where z̄ is an average PSN activation. Here and throughout the article, the PSN activation is discretized with a timestep τ before processing, i.e., z̄ = 〈z(m · τ)〉. The goal is then to ensure that

z̄(G p) ≈ σ(G p − G 0) = ∂ u G, (9)

which is an optimization problem for the three variables τ, G 0 and ∂ u G. With it solved, we find that the PSN activation function matches the logistic activation function of LIFS neurons well. Here, τ = 11 ns, which corresponds to a sampling rate of approximately 0.1 GHz. Based on the similarities between the expressions for membrane potentials of LIFS neurons (Eq. (3)) and PSNs (Eq. (6)) we find the following translation rules between the bias and the pumping strength:

b k = (G p,k − G 0) / ∂ u G. (10)

Connections between LIFS neurons induce a change of their membrane potentials according to Eq. (3). Consider two LIFS neurons connected unidirectionally, i.e., only W 12 ≠ 0, with b 2 = 2. This way, the second neuron acts as the source of spikes, while the first is the recipient. Sweeping b 1 and W 12 we obtain a set of curves p(z 1 = 1 | b 1, W 12) shown with lines in Fig. 4b. We now aim for a similar behavior in PSNs. We convert b 1 and b 2 to PSN pumping strength according to Eq. (10). We then sweep κ 12 to obtain a set of curves p(z 1 = 1 | G p,1(b 1), κ 12). The goal is to find a proportionality coefficient ξ such that

p(z 1 = 1 | G p,1(b 1), ξκ 12) ≈ p(z 1 = 1 | b 1, W 12). (11)

The result is given in Fig. 4b. We find that for limited biases and weights, PSNs replicate the behavior of LIFS neurons well. Based on Eq. (11), we find the translation rule for weights:

W kj = ξκ kj. (12)

To conclude, we have found that dynamics of PSNs and their interaction imitate that of neurons of a BM. In the following sections, we investigate the accuracy of the imitation by considering a series of tasks of increasing complexity.

Optical sampling from Boltzmann distributions
To investigate the accuracy of sampling with PSNs, we first consider sampling from predefined Boltzmann distributions over a small set of binary random variables. Following 13, biases and weights were drawn from Beta distributions, b k ∼ 1.2(B(0.5, 0.5) − 0.5) and W kj ∼ β(B(0.5, 0.5) − 0.5), where β controls the range of weights and is either 0.6, 1.2, or 2.4. The generated biases and weights are translated to PSN parameters using Eqs. 10 and 12. The sampled distributions are compared to the target distributions p* by means of the Kullback-Leibler divergence D KL(p ∥ p*).

First, we optimize the PSP timescale controlled by τ U. On one hand, τ U must be small enough such that a PSP does not last longer than the refractory period τ. On the other hand, too short τ U will make the PSP short, but strong, which can break the operating regime assumptions (see Sec. “Simplified gain equation” in the Methods). In LIFS neurons, the PSP timescale is ideally slightly shorter than the refractory period 13; we use this as a starting point. In Fig. 5a, we sweep τ U and track the sampling accuracy for all considered β. We find that the optimal τ U is approximately 0.37τ. Figure 3b shows a PSP with such a timescale, and indeed, the PSP becomes negligible after t = τ.

We proceed to sample from sets of 10 random Boltzmann distributions for different β. Figure 5b shows the convergence of the sampling procedure towards the target distribution. We find that within each set, the convergence is almost identical. In Fig. 5d, we compare a distribution of samples to the exact distribution for β = 2.4, and we note a very close match. Figure 5c shows a raster of spikes during this sampling.

For the next test, we split the five neurons into two arbitrary groups, Y and X. The first group is clamped to an arbitrarily chosen state y, and the neurons in the second group are free; their probability distribution is the conditional probability distribution p(X | Y = y). In PSNs, we clamp the state by significantly reducing or increasing the pumping strength. During sampling, we collect samples of the free neurons and this way, estimate p(X | Y = y). In other words, the PSN network samples from the conditional probability distribution of the BM. Figure 5e shows the close match between the correct conditionals and those sampled with our PSNs. We therefore conclude that a network of PSNs can sample accurately from a wide range of Boltzmann distributions and their conditionals over small state spaces.

Optical sampling from arbitrary distributions
Boltzmann distributions are only a subset of all possible probability distributions over binary variables. However, the fully visible sampling networks described above can be extended by adding hidden neurons that are not observed during sampling. Consequently, the probability distribution of the visible layer becomes a marginal distribution over the full state space, which can, in principle, take any shape, given a large enough hidden space. To simplify training and improve convergence, a hierarchical network structure is preferable, with no horizontal connections within individual layers 15. Here, we emulate such a two-layer restricted Boltzmann machine (RBM) with an equivalently structured PSN network (see Fig. 6a). Given enough hidden neurons, an RBM can sample from any distribution with arbitrary precision 53. Consequently, by implementing such an RBM with PSNs, optical sampling from arbitrary distributions can be achieved.

Here, we choose a target probability distribution p*(z) with 4 binary variables, which can be fully defined by 2 4 = 16 probabilities, which we generate by sampling 16 numbers from the inverse continuous uniform distribution over [0, 1] (whose probability density is f(x) = x −2 for x ∈ [1, ∞] and zero elsewhere) and normalizing them such that their sum is unity. The probability distribution is shown in Fig. 6c. Using the contrastive divergence algorithm (see Sec. “Contrastive divergence” in the Methods), we train an RBM with 4 visible and 10 hidden neurons. Its sampled distribution over visible neurons p RBM(z) is shown in Fig. 6c.

The parameters of the trained RBM were then transferred to a network of 14 PSNs. Its sampled distribution over visible neurons p PSN(z) is also shown in Fig. 6c. We find that the accuracy of the PSN network is very close to that of the RBM, and sampling from the target distribution p*(z) and several conditionals shown in Fig. 6d is correspondingly accurate.

Optical probabilistic inference
In Bayesian statistics, a probability represents the degree of belief in an event. Bayesian inference is an update of this belief based on additional information about this event. For example, assume a scenario in which the receiver of an image expects any of three digits (here: 0, 3, and 4) to arrive with equal probability. If the received image is highly corrupted, such that only a few of the received pixels are known to be correct, what can the receiver infer about the transmitted image? In this section, we let PSN networks provide the answer.

More specifically, we chose three images of digits 0, 3, and 4 from the MNIST dataset 54, rescaled to 12 × 12, and with brightness rounded to zero or unity. These images were mapped to a fully connected network of 144 PSNs, with one pixel assigned to each neuron. We then proceed to train the PSN network to store the prepared images, meaning the network will have three most likely states that form the images of the digits, with each state being equally probable. For this, we first train an equivalent BM using wake-sleep (see Sec. “Probabilistic inference training” in the Methods) and map the resulting parameters to the PSN network.

First, we assess the mixing capability of the network by observing its dreaming phase (corresponding to the sleep phase during training). Without external input, the PSN network switches randomly between states, forming the three images with approximately equal probability, which represents our prior belief about the possible images. Figure 7a shows a two-dimensional projection of PSN samples onto vectors corresponding to the images (see Sec. “Probabilistic inference visualization” in the Methods). We found that the result is close to that of the BM, and in most cases, the network switches between the numbers every few samples (Fig. 7b), indicating good mixing between the states.

Next, we assess the inference capability of the network in a pattern completion scenario. By applying additional bias to a few neurons that are only active for two out of the three patterns, we provide informative, but limited input to the network. We chose five pixels that are black for 0 and 3, but not 4 (Fig. 7c), and applied an additional positive bias to the corresponding neurons. The biases of other neurons are unchanged, i.e., no information is given. Such an input is ambiguous with respect to 0 and 3, but incompatible with 4. As a result, the PSN network only generates complete images of 0 and 3, and randomly switches between them.

For guided dreaming, label neurons were clamped for 200τ. To enforce top-down control (from labels to pixels), we strengthened the weights between the hidden and label layers by a factor of 2, while marginally reducing the others by 10% (see Discussion). In Fig. 8e, we show how the PSN network can thereby be used for generating images from all learned classes.

Learning from data
Sampling from arbitrary distributions implies the ability to sample from distributions dictated by real-world data. Here, we learn a generative model of handwritten digits based on the MNIST dataset 54. Following ref. 20, we round the brightness of each pixel to the minimum or maximum.

To work with this dataset, we use a hierarchical sampling network: an RBM with three layers: visible, hidden, and label (Fig. 8a). The brightness of pixels is mapped to the activity of a neuron in the visible layer. The activity in the label layer shows which digit is represented in the visible layer at the current point in time. For example, if visible neurons form an image of 0, only the label neuron corresponding to a 0 will spike.

In this section, we consider three tasks: completion, guided dreaming and classification. For completion, visible neurons are clamped to the brightness of corresponding pixels except for a bottom-right quadrant, which is assumed obscured and remains free alongside other neurons. The task is for the free visible neurons to complete the obscured quadrant. For guided dreaming, one label neuron is clamped, and the visible neurons are expected to form an image of a corresponding digit. For classification, visible neurons are clamped to the brightness of corresponding pixels, and the most active label neuron on average is taken as the answer (winner-takes-all).

For these tasks, we use the BM parameters from ref. 55 (Fig. 8b, c), as they were already optimized for LIFS networks and are therefore expected to be favorable for an implementation with PSNs. This network is composed of 1194 neurons: 784 in the visible, 400 in the hidden, and 10 in the label layers, respectively. Its implementation with PSNs thus represents a scaling test for our general approach.

The simulation of PSNs during the tasks starts with a burn-in phase, where no input is provided. Then, the inputs are provided sequentially, with each subsequent input following immediately after the previous one.

For pattern completion with PSNs, 20 samples were drawn for each image. This is a difficult inference task, as the network should not simply produce an average image, but instead need to adapt to the style of each individual input sample. Figure 8d shows the results for a few images from the MNIST testing dataset. We found that in most cases, the obscured parts were completed to a large degree of accuracy. In Fig. 8e, we show how the PSN network can thereby be used for generating images from all learned classes.

While BMs are designed primarily as generative networks, especially when trained purely with contrastive Hebbian methods and without dedicated backpropagation-based fine-tuning, input classification can still be viewed as a form of Bayesian inference; thus, it is instructive to compare classification performance between the original BM and its PSN implementation.

Simulating the processing of 10,000 images of the testing MNIST dataset with a PSN network of this size is computationally intensive; we therefore drew samples until the saturation of the classification convergence curve (Fig. 8f). In this case, we drew 50 samples for each image, yielding an average classification accuracy of 87.0%. The RBM is much less demanding; we could thus draw 500 samples, obtaining 87.8% accuracy. Their confusion matrices are compared in Fig. 8g; we note their similarity, as well as the only marginal performance loss caused by the exchange of substrates.

Interestingly, during the completion task, the networks perform classification as well. However, the accuracy of the PSN network is reduced to 75.0%, compared to 85.6% for the RBM (Fig. 8f); we expect this to be mitigated by further fine-tuning or direct in-situ training of the PSN network (see Discussion). However, in both cases, we observe that the PSN network only needs a few samples to converge to a solution, which is faster than the RBM and can prove beneficial in time-constrained scenarios.

Methods
Nanolaser with quantum wells
The rate Eqn. (4) Describe the interaction of light with quantum dots 45,46. These semiconductor nanostructures, akin to artificial atoms, localize carriers within a few nanometers and therefore are well modeled with a two-level electronics system. However, this work builds on the experimental results reported in ref. 39, where the gain material consists of quantum wells. There, unlike in quantum dots, the gain depends nonlinearly on the density of carriers in the conduction band (i.e. the population of excited dipoles in our model divided by the volume of the section). This nonlinear dependence is approximated by a piecewise linear function. The slope is larger at lower carrier density, with a ratio χ g. To ensure that n e ≤ n 0, we corrected the pumping term in Eq. (4) by replacing n e → 2n e / (χ g + 1).

The pump rate at threshold γ p thr is computed from Eq. (4) by assuming a steady state without noise, G = 0 and S ≈ 0, which leads to γ p thr = γ t n e thr / (n 0 − 2n e thr) = (χ g + 1) γ t, where n e thr = (γ t + γ r n 0 + γ r,a n 0,a) / 2γ r. Substituting values from Supplementary Table 1, we find n e thr ≈ 1.47 × 10 6 and γ p thr ≈ 1.31 × 10 9 that corresponds to a current of γ p thr n 0 e ≈ 215 μA, where e is the elementary charge.

Nanolaser noise model
The stochastic Eq. (4) are composed of deterministic (drift) and stochastic (diffusion) parts that can be represented in matrix form:

du = μ(u, t) dt + σ(u, t) dW t, (13)

where u = [S, n e, n a] T, and W is the Wiener process. In this work, the diffusion term follows the approach in [ref. 47, A.13.1.2] based on the McCumber noise model 63 and is comprised of five Langevin forces corresponding to the following groups of processes:
1. electron-photon interaction in the gain section,
2. same, in the SA,
3. electronic processes inside the gain section,
4. intrinsic optical loss,
5. electronic processes inside the SA.

The Langevin forces are stochastic processes represented with Wiener processes with zero average and cross-correlation 2D ij; for i = j, the latter represents the autocorrelation or noise spectral density. For the groups of processes described above, we find

2D SS e = γ r n e S + γ r n e,
2D SS a = γ r,a n a S + γ r,a n a,
2D ee = γ p n 0 + (γ r − γ p) n e / (1 + χ g),
2D SS o = γ t S,
2D aa = (γ t,a − γ r,a) n a,

where n e, n a and S are the averaged (not stochastic) values. This way, the diffusion matrix becomes:

σ(u, t) = [ 2D SS e 2D SS a 2D SS o 0 0; 0 0 0 2D ee 0; 0 0 0 0 2D aa ] 1/2. (15)

The equations are integrated using the SKSROCK solver from the DifferentialEquations.jl library 64 in the Julia programming language 65.

Simplified gain equation
Here, we derive the Eq. (6). Consider a network of PSNs. A derivative of the gain of a PSN is given in Eq. (4):

Ġ = 2γ r ṅ e + 2γ r,a ṅ a. (16)

Assume the PSN is not currently emitting a spike, i.e. its gain does a random walk shown in Fig. 2c, in which case, S ≈ 0. Then, substitute the derivatives of electron populations from Eq. (4):

Ġ ≈ −2γ r γ p n e + 2γ r(γ p + Δγ p(t))(n 0 − n e) − 2γ r,a γ t,a n a + 2γ r F e(t) + 2γ r,a F a(t) = −2γ r n e(γ p + γ t) + 2γ r γ p n 0 − 2γ r,a γ t,a n a + 2γ r Δγ p(t)(n 0 − n e) + 2γ r F e(t) + 2γ r,a F a(t). (17)

where Δγ p(t) is given in Eq. (5). Then, replace 2γ r n e = G + γ t − γ r,a(2n a − n 0,a) + γ r n 0:

Ġ ≈ −(γ t + γ p)(G + γ t − γ r,a(2n a − n 0,a) + γ r n 0) + 2γ r γ p n 0 − 2γ r,a γ t,a n a + 2γ r Δγ p(t)(n 0 − n e) + 2γ r F e(t) + 2γ r,a F a(t). (18)

Rearranging the terms we find

Ġ ≈ −(γ t + γ p)G + 2γ r,a(γ t + γ p − γ t,a) n a − (γ t + γ p)(γ t + γ r n 0 + γ r,a n 0,a) + 2γ r γ p n 0 + 2γ r(n 0 − n e) Δγ p(t) + 2γ r F e(t) + 2γ r,a F a(t). (19)

The second term on the first line is negligible compared to the first, as during the random walk n e ≪ n 0. The terms on the second line can be considered a drift term:

G p = −γ t − γ r n 0 − γ r,a n 0,a + 2τ G γ r γ p n 0,

where τ G = 1 / (γ t + γ p). This way, the dynamical equation becomes

dG = −(G − G p) / τ G dt + 2γ r(n 0 − n e) Δγ p(t) dt + σ G dW. (21)

Here, the interaction term can be nonlinear as n e changes due to incoming spikes; we assume that such an interaction is weak. Moreover, during the walk, the gain section is close to transparency, i.e. n e ≈ n 0 / 2. Finally, we find

dG = −(G − G p) / τ G dt + 2γ r n 0 Δγ p(t) dt + σ G dW. (22)

Contrastive divergence
We use contrastive divergence to train an RBM to sample from an arbitrary distribution p*(z). In this section, we provide a brief summary of the procedure; for a more detailed introduction, we refer to ref. 67.

We start with arbitrary weights W, visible biases b v and hidden biases b h. Then, we iterate over all possible states of z and set the target probabilities to p*(z). For each z m, we set the initial state of visible neurons v(0) = z m. Following the RBM sampling procedure, we then sample the hidden neurons h(0) given v(0), followed by a resampling of the visible states v(1) given h(0), etc. This procedure is repeated K times, resulting in a sequence of states v(0, ..., K) and h(0, ..., K). The final states v(K) and h(K) are then used for updating the parameters. In our case, K was also gradually increased from 1 up to 20 as the sampling accuracy improved.

For each z m, the state-specific parameter updates are then calculated, but not yet applied. The weight updates are given by

ΔW m = p*(z m)(v(0) ⊗ h(0) − v(K) ⊗ h(K)), (23)

where ⊗ is the outer product; for the bias updates, we have

Δb v m = p*(z m)(v(0) − v(K)), (24)
Δb h m = p*(z m)(h(0) − h(K)).

After iterating through all m, the weights and biases are updated according to ΔW ∝ ∑ m ΔW m, Δb v ∝ ∑ m Δb v m, and Δb h ∝ ∑ m Δb h m.

Probabilistic inference training
The approach used here follows 13. A fully visible BM was trained to store three digits – 0, 3, and 4 – taken from the MNIST dataset and scaled down to 12 × 12. The intensity of each pixel, which ranged from zero to unity, was rounded to 0.05 or 0.95. This way, the intensity of the pixels defines the target statistics 〈z k〉* and 〈z k z j〉* for the BM. We start with arbitrary weights W kj and biases b j and collect a sufficient number of samples to estimate 〈z k〉 and 〈z k z j〉. Then, we refine the BM parameters using the following update rules: Δb k ∝ 〈z k〉* − 〈z k〉 and ΔW kj ∝ 〈z k z j〉* − 〈z k z j〉. Then, sampling and refinement is repeated until a satisfactory result is achieved.

Probabilistic inference visualization
Figure 7a, d visualize high-dimensional states of the network using a star plot of their projections. The projection procedure follows 13 with a minor modification.

The three axes show the basis vectors B representing pixel intensities of the target images:

〈z〉 = (B 0, B 3, B 4) T. (25)

The basis vectors are normalized: ∥B i∥ = √∑ j B ij 2 = 1. After a simulation, the network states are discretized once per τ and are projected onto this basis:

z 034(t m) = (B 0 z(t m), B 3 z(t m), B 4 z(t m)) T. (26)

This projection is, in turn, projected onto the two-dimensional plane z proj(t m) = M proj z 034(t m), where

M proj = [sin(φ 0) sin(φ 3) sin(φ 4); cos(φ 0) cos(φ 3) cos(φ 4)] (B T B) −1. (27)

where B = (B 0, B 3, B 4) and (φ 0, φ 3, φ 4) = (0, 2π/3, 4π/3). The (B T B) −1 term accounts for the non-orthogonality of the basis vectors and was not used in 13.

Due to a limited number of unique projection results, a standard scatterplot cannot represent clustering of a large number of projections in a single spot. Instead, we opt for a scatterplot with large near-transparent markers. This way, multiple overlaid projections result in a darker marker.

The proximity of the network states to the target states is represented by the distance between z proj(t m) and 〈z〉. Judging by the clustering of dark markers around the projections of the target images, the network spends the majority of time around the corresponding states.

MNIST power consumption
First, we estimate the power necessary to reach the operating point, which, in PSNs, is close to the threshold (see Fig. 4a). For the parameters considered here, it is γ p thr ≈ 215 μA (see Sec. “Nanolaser with quantum wells” in the Methods). The current due to BM bias will change the electrical current by a few percent (see Fig. 4a), which we consider negligible. The typical bias voltage for III-V semiconductor lasers is 1 V, which results in approximately 215 μW per PSN. For MNIST with 1194 nanolasers, 250 mW would be required, although roughly 20% of PSNs never spike and require no pumping.

Next, we estimate the power necessary for interaction between PSNs. Assume the j-th PSN has emitted a spike that is then routed to the k-th PSN. According to Eq. (4), Δn e,k = ∫ Δγ p,k(t)(n 0 − n e) dt. Assume the k-th PSN is undergoing a random walk, during which n e ≈ n 0 / 2, and with Eq. (5), we find

Δn e,k ≈ n 0 / 2 κ kj ∫ ∫ LPF(t − t*) S j(t*) dt* dt. (28)

Since the low-pass filter preserves energy, the two integrals become the total photon count of a spike; overall, |κ kj| n 0 / 2 of the photons are absorbed.

Now, assume the spike is subject to losses in the transmission chain, and the k-th PSN only receives a fraction ζ of the expected photons; thus, ζ −1 times more photons must be sent to compensate. Summing over all PSNs connected to the k-th, we find that ζ −1 ∑ j |κ kj| n 0 / 2 of the spike photons must be sent to the coupler. If this value exceeds unity, more photons are required than the spike contains, and amplification is therefore required. With the parameters considered in this work, the translation rule for weights is κ kj ≈ (0.2 / n 0) W kj. Then, the k-th neuron must send 0.1 · ζ −1 ∑ j |W kj| of the spike power, and the required amplification multiplier is

max(1, 0.1 · ζ −1 ∑ j |W kj|) − 1. (29)

The weighting scheme is expected to be the dominant source of optical losses. Here, we follow the conclusion in the previously mentioned work on a 240 × 240 switch 60. There, losses are dominated by waveguide crossings: 0.016 dB for each. For MNIST, a spike traveling between the visible and the hidden layers would pass over 784 + 400 − 1 crossings, which corresponds to 18.9 dB loss, and, similarly, 6.5 dB for the hidden and label layers. This loss is expected to be significantly reduced by the use of multilayer bus waveguides. As a result, propagation losses become dominant, but the use of shallow-etched silicon rib waveguides allows losses down to 0.1 dB cm −1. For the visible-hidden coupler, it corresponds to 0.9 dB (= 0.1 dB × 75 μm cell −1 × 1184 cells) and 0.3 dB for the hidden-label coupler. Still, with the technology in 60 and such waveguides, for the visible-hidden coupler we find ζ ≈ −19.8 dB and for the hidden-label coupler ζ ≈ −6.8 dB. Then, using Eq. (29). We find the average amplification
- 24.2 dB from the visible to the hidden layer,
- 27.1 dB same, in reverse,
- 0.002 dB from the hidden to the label layer,
- 14.2 dB same, in reverse.

To get amplification in power units, we must first find the power emitted by a PSN. In Fig. 1b, we find that S(t) can be approximated with a Gaussian with FWHM ≈ 30 ps. The peak photon count in the cavity is S max ≈ 8e4. Assuming that intrinsic optical losses are negligible compared to losses due to waveguide coupling, photons escape to the waveguide at a rate γ S(t). Considering that photons escape in both directions of the waveguide (see Fig. 1a), the number of useful photons is halved, γ S max FWHM π / 2 ≈ 4.2e5 which at C-band (1550 nm) corresponds to 54 fJ. The emission power depends on emission frequency, which can vary significantly depending on the dynamics of the entire network. For simplicity, assume that all PSNs operate with zero membrane potential, i.e., p(z = 1 | b = 0) = 0.5. Then, the average power emitted by a PSN is 54 fJ × 0.5 × 0.1 GHz ≈ 2.7 μW. Therefore, amplification of the signals from the visible layer to the hidden layer would require 2.7 μW × 784 × 24.2 dB ≈ 560 mW. Following this logic, we find:
- 560 mW from the visible to the hidden layer,
- 1100 mW same, in reverse,
- 2 mW from the hidden to the label layer,
- 55 mW same, in reverse.

Therefore, for MNIST classification, where the visible layer is clamped, amplification of spikes coming from the hidden to the visible layer is unnecessary, and 560 mW + 2 mW + 55 mW = 617 mW is required for amplification in total. In contrast, during guided dreaming, the label layer is clamped and 560 mW + 1100 mW + 55 mW = 1715 mW are required.

For photonic accelerators, floating point operations (FLOP) per second (FLOPS) are often provided as a figure of merit. For the network of PSNs, it is appropriate to find how many FLOPs is required by the BM it replicates. For simplicity, the hierarchical sampling network used here can be considered an RBM, where visible and label neurons are in one layer, and hidden neurons are in the other layer. The neuron activation function is a sigmoid σ(x) = exp(x) / (1 + exp(x)). Assuming that exp can be computed in 1 FLOP, computing the sigmoid costs 3 FLOPs. A neuron connection is a multiplication and requires 1 FLOP.

During MNIST classification, 65% of neurons are clamped, and the estimation of FLOPS will not represent the true capability of the PSN network. Instead, consider drawing one sample during the guided dreaming task. First, the activations of the visible and the label neurons are weighted by a 400 × (784 + 10) matrix, which requires 2 × 400 × 794 FLOPs. Then, the activations of the hidden neurons are computed by adding 400 biases and spending 3 × 400 FLOPs for the sigmoid. The same logic holds for the backward step. Therefore, one sample requires 2 × 400 × 794 + (1 + 3) × 400 + 2 × 794 × 400 + (1 + 3) × 794 = 1275176 FLOPs. Considering the network of PSNs draws one sample per τ = 11 ns, it performs an equivalent of 1275176 FLOP / 11 ns ≈ 116 TFLOPS during this task. With the previously found approximate power draw of 2 W, we also find the energy efficiency of 58 TFLOP J −1 or 17 fJ FLOP −1.

To conclude, the power consumption is composed of PSN pumping and spike amplification. For MNIST tasks, the former is 250 mW for the parameters used in this work. The latter, assuming crossbar technology used in ref. 60, varies between 600 mW and 2 W depending on the task. The performance is estimated to reach 116 TFLOPS with the energy efficiency of 58 TFLOP J −1 or 17 fJ FLOP −1.

t-Stochastic Neighbor Embedding
t-Stochastic Neighbor Embedding (t-SNE) is a technique for dimensionality reduction of highly-dimensional datasets 68. In particular, it performs a nonlinear projection of data onto a two-dimensional plane, where the proximity of the projections represents the similarity between the features of the data points.

During guided dreaming, the PSN network generates images of digits. Due to stochasticity, the generated images are not identical, but still have similar features, which are used by t-SNE to cluster the images for better visualization.

Simulation code details
In order to apply a BM bias b to a neuron in the PSN network, instead of Eq. (10) We translate the bias to PSN pumping directly using

γ p = γ p,0 + ∂ u γ p b, (30)

which is equivalent as G p and γ p are linearly dependent. This is a more practical approach that allows using Fig. 4a directly.

All tasks considered in this work follow the same pattern (see also Supplementary Fig. 2):
1. train a BM to solve a problem,
2. translate BM parameters to PSN parameters,
3. simulate the PSN network,
4. discretize the PSN output,
5. process the result as if the output is from a BM.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript proposes and analyzes the concept of spiking nanolasers, i.e. semiconductor lasers that emit ultrafast\noptical pulses, and demonstrates a model for computing based on sampling-based neural computation. The core idea is to\nuse small photonic crystal nanolasers, each behaving as a stochastic spiking neuron, to perform Bayesian inference through\nsampling from learned probability distributions (e.g., Boltzmann machines). The authors provide a theoretical treatment\ncomparing these nanolaser-based neurons to conventional leaky integrate-and-fire sampling (LIFS) neurons, show how to\nmap typical Boltzmann machine parameters onto spiking nanolasers, and then illustrate the approach on generative tasks\nand pattern completion (e.g., small-scale MNIST).\nKey contributions include:\ni) A derivation that links nanolaser gain dynamics and standard stochastic neuron models.\nii) Numerical experiments suggesting significant speed and potential power-efficiency gains over purely electronic spiking\nsystems.\niii) Demonstrations of sampling accuracy for small toy problems (MNIST-like tasks).\nOverall, the paper provides an interesting proposal for ultrafast, optically based neuromorphic systems. The results are\ntimely, given the growing interest in photonic computing for machine learning, and the manuscript’s focus on sampling-\nbased models (rather than purely feedforward neural nets) highlights unique advantages of spiking neural substrates.\n2. Significance and Relation to Established Literature\nOptical spiking neurons could become a promising direction in neuromorphic engineering due to their potential for high-\nspeed operation and parallel data transmission. Compared with other photonic neuromorphic proposals, the authors’\nemphasis on sampling, as opposed to deterministic classification, addresses an important niche. This aspect could be\nfurther empathized in the introduction.\nThe paper references prior demonstrations of spiking semiconductor lasers and places the present approach in context with\nexisting neuromorphic devices like BrainScaleS and other hardware for spike-based sampling.\nProvided these speed and energy claims hold up in experimental settings, this work may have significant impact on the\ndesign of specialized photonic or hybrid photonic-electronic chips for Bayesian inference, optimization, and other tasks that\nbenefit from fast random sampling. It builds upon established knowledge of Boltzmann machines and spiking neural\nnetworks but extends it into the realm of ultrafast photonic hardware.\n3. Support for Conclusions and Claims\nSpiking Mechanism & Neural Analogy: The theoretical equations linking the gain parameter of the laser to a membrane-\npotential-like stochastic variable appear consistent with known rate equations for semiconductor lasers. The authors’\nsimulations further support that the excitability and refractoriness can mirror spiking neuron behavior.\nSampling Accuracy: Simulation results show that mapped Boltzmann parameters lead to correct sampling distributions in\nsmall networks, as well as decent performance on pattern completion for MNIST-like inputs. These results support the claim\nthat photonic spiking neurons can implement sampling-based inference.\nEnergy and Performance Estimates: The authors provide speed estimates (on the order of gigahertz spiking rates), as well\nas approximate energy calculations that suggest orders-of-magnitude gains over purely electronic solutions. While these\nnumbers are plausible, an experimental prototype or a more explicit breakdown of power budgets (e.g., how pumping,\nwaveguide coupling, photodiode readouts, and any amplification circuitry scale with network size) would further validate the\nclaims.\nAdditional Evidence: The paper convincingly motivates the theoretical viability of photonic sampling. However, real-world\nconstraints, fabrication tolerances, optical losses, device mismatch, and large-scale routing complexities, might introduce\ndeviations from idealized modeling. Some more in-depth discussion of those potential challenges, or references to relevant\nintegration technologies, would strengthen the claims.\nPlease address in the discussions points b, c, d and how these could be overcome and or limit the current work.\nPotential Flaws and points that require revision\na) Energy Estimation Details: The manuscript’s power estimates for large-scale networks could benefit from more specificity\non amplifier or detector overheads, along with tolerance analysis of waveguide losses. As it stands, the rough estimates\nseem encouraging but may need tighter bounds or a more explicit breakdown.\nb) Definition of Abbreviations: There are multiple instances where abbreviations appear before being defined. A\ncomprehensive glossary or ensuring every acronym is defined upon first appearance would significantly improve clarity. (eg,\nSPS, PNS page 2 etc..)\nc) Methodological Completeness: While the theory and simulation approach are well described, readers from the photonic\nhardware community may want more details on device geometry or coupling schemes. Additional clarity about how well the\nauthors expect the assumptions (e.g., linear summation of spikes, negligible inter-laser coherence) to hold in practice would\nhelp reproducibility.\n5. Methodology and Reproducibility\nSoundness: The approach, using the known rate equations for semiconductor lasers augmented with noise terms, is\nstandard and appears methodologically sound. Mapping Boltzmann parameters onto laser parameters follows from typical\nsampling-based spiking neuron theory.\nReproducibility: The methods section largely includes the relevant equations and standard parameter sets, enabling future\nresearchers to replicate the simulations. Including a fully annotated list of parameters for each experiment (e.g., for the\nMNIST demonstration) and clarifying how the code is structured would enhance reproducibility.\nIn summary, the authors present a theoretically sound proposal. The authors make a strong case for how such ultrafast\noptical spiking neurons could perform sampling-based inference with notable improvements in speed and potentially power\nefficiency. Despite some points that would benefit from more details, particularly around practical power budgets and large-\nscale integration, the paper is overall a valuable contribution.\nRecommendation: I recommend revisions to polish the manuscript, clarify energy estimates, ensure consistent definitions of\nabbreviations, and possibly strengthen the discussion around potential integration and scaling challenges.", "reviewer_2": "This work reports on the mapping from a photonic spiking neural network (PSN) to a Boltzmann machine (BM), and\nnumerically demonstrates efficient computational tasks based on Bayesian inference. Among the myriad of today’s optical\ncomputing approaches, this paper addresses the neuromorphic spinking response of semiconductor lasers with saturable\nabsorbers, as a consequence of excitability. While excitability in semiconductor lasers is quite an old topic, intensively\ninvestigated in the past two decades, the advent of integrated photonic devices on the one hand, and the groundbreaking\nadvancements in AI on the other hand, have boosted the research in this domain. In this respect the present manuscript\nrepresents quite a significant advancement. In particular, the small footprint and low consumption of photonic crystal\nnanolasers are claimed to offer a clear advantage over in-silico substrates, both in terms of speed and energy budget.\nNevertheless, there are a few important issues that prevent me to positively recommend this work for publication in Nat.\nComm. (points 1-4 below).\n1) Broad audience and general interest. This manuscript is difficult to be understood by a non-specialist. It assumes a\nreader’s background on Bayesian statistics and Boltzmann machines, which is far to be the case even in the community of\noptical neuromorphic computing. It is then of primary importance to introduce the main concepts of Bayesian inference,\nBoltzmann machines and sampling (both non-spiking and spiking), with a basic example (prior to spiking models) which\nwould help the reader to a large extent. Otherwise, we might consider this work to be better suited to a more specialized\njournal (JCP or CPC). In my opinion, introducing these concepts is substantially more important than the details of the\nmapping (“Membrane potential of spiking nanolasers”): most of this subsection could be moved to a supplementary material,\nretaining only the fundamental result in the main text.\n2) Overall clarity. Many points remain obscure to me (many of them in connection to point 1). For instance: it is not clear how\nthe training is performed, i.e. what parameters (weights?) are tuned, and how. Another example is the 2D representations of\nnetwork states (Fig. 7.a and 8.e): the vast majority of readers ignore what T-SNE is, so a hint to interpret the representation\naxes is needed. Furthermore, sentences like “…the PSN network correctly samples from the prior” (l 410) or “…by sampling\nfrom conditional probability distributions” (l. 345), or “the probability of each state sampled from the inverse…” (l. 378) are\nvery technical and specific from Bayesian statistics, so they need further clarification here. How PSP is implemented is not\nclear neither.\n3) The optical spiking advantage. It is claimed that the spiking dynamics is advantageous both in terms of energy saving and\nrobustness to noise. While I acknowledge this as general statement, which may hold true in neuroscience, it is not clear to\nwhat extent it applies to the particular photonic system studied here. Specifically:\na. Power. CW laser-threshold powers are even smaller than those of pulsed devices.\nb. Bandwidth. Frequency multiplexing is not clear in the context of the proposed\nphotonic architecture, and should be commented.\nc. Losses. Interconnect matrix is based on waveguide crossbar arrays. Then the scaling of crossing losses is a potential\nissue in large networks, which is not discussed.\n4) The context of excitability in lasers. The references in the Introduction do not fairly reflect the historical developments. The\nfirst cited works on excitability (Refs. [23-25], Wünsch et al., Romeira et al., Prucnal et al.) do not accurately retrace the\nseminal studies on optical excitability. The first report on excitability was by F. Plaza et al. (EPL 38, 85, 1997) in a CO2 laser\nwith SA, followed by experiments in semiconductor lasers (M. Giudici et al., Phys. Rev. E 55, 6414, 1997) and Nd-YAG\nlasers with SA (Larotonda et al., PRA 65, 033812, 2002). In the context of photonic crystals, the first demonstration in a\nnanocavity was reported in Brunstein et al, PRA 85, 031803 (2012).\nMinor comments:\ni. Abstract. Saying that speed is an advantage of optics w.r.t electronics is not true in general, as an electronic signal also\npropagates at the speed of light in metallic wires. In addition, in the abstract (l. 5), what do authors mean by “conventional\nsemiconductor devices”?\nii. Possible confusion in the notation: $\\beta$ is used as spont. emission factor of nanolasers, and as a parameter for\nprobability distributions.\niii. Absolute/relative refractoriness are assumed to be well-known concepts, yet they’re most likely ignored by the\ngeneral reader.\niv. Line 264 should read: “Second, changes of $n_e$ due to incoming spikes need to be small…”", "reviewer_3": "The manuscript demonstrates the feasibility of general spike-based sampling using a network of spiking lasers and discuss\nseveral examples of this sampling operation for the realization of complex machine learning tasks. The work is very\ncomplete, very well written and should be interesting to a large community ranging from machine learning specialists to\nphotonics experts. The results include foundational points such as how to map the spiking probability to BM weights\nand very advanced tasks such as \"guided dreaming\". For the completeness of the work, its novelty, informative content and\ngeneral quality, I warmly recommend publication of the manuscript.\nHowever, a few points deserve clarification in my opinion:\nThe coupling between nodes is considered to take place via embedded photodetectors, aiming at incoherent transmission of\nspikes. Is coupling through saturation of the absorber section impossible or impractical?\nThe model uses variabe S to describe a photon number, which must therefore be positive. What prevents S from becoming\nnegative in presence of the stochastic term F? This is particularly relevant considering that most of the time the system is not\nemitting a spike thus in absence of noise S=0 (figure 1b).\nFigure 2b shows time traces of the normalized gain, which is, as expected mostly slightly below 0. However, at the onset of\nspikes, the gain depletion is anticipated by a very large increase of the gain (particularly visible for instance on spikes 3,4,5),\nwhich takes place at a time scale which is comparable to the following depletion (which is mediated by the photon number\ngrowth). What is the origin of this very sharp increase? Could it be related to the previous point?\nThe \"wake sleep\" algorithm *must* be briefly described, with a reference for more information, if the work is to be readable\nalso to non-experts in stochastic neural networks training.\nPoints to be considered:\nThe \"optical probabilistic inference\" and \"learning from data\" sections seemed to me as \"logical\" specific applications of the\nsection \"optical sampling from arbitrary distribution\". If it is the case, I fear these sections do not really help the reader.\nI understand they constitute \"bright\" points and their implementation is probably worth a description but if they are really direct\nconsequences of \"optical sampling from arbitrary distribution\" the reader would benefit from being informed of that very\nclearly.\nMinor points:\nfigure 1a: is a spike really shorter than the laser resonator as suggested by fig 1a? If yes, a PDE system would really be\nrequired to describe the physics. If not, then a single mode model as used here is fine but then figure 1a is a bit misleading.\n\"alpha shaped\" is used here and there in the manuscript but not defined\nThe connection between BMs and spiking networks cites 11-13, I believe Hinton, G. E., & Brown, A. (1999). Spiking\nboltzmann machines. Advances in neural information processing systems, 12. should be credited."}, "author_response": "Response to Reviewers\n\nReviewer #1\n\nAuthors’ response: Thank you for this very encouraging assessment!\n\nSignificance and Relation to Established Literature\n\nAuthors’ response: Thank you for the suggestion. We have amended the introduction with several additional comments and links to recent and ongoing research (lines 65 ff., 75 ff. and 89 ff.).\n\nSupport for Conclusions and Claims\n\nAuthors’ response: Thank you for your suggestion. We now provide a detailed breakdown of power budgets in the “MNIST power consumption” section in the Methods and have amended the Discussion accordingly (line 655 ff.).\n\nAdditional Evidence: The paper convincingly motivates the theoretical viability of photonic sampling. However, real-world constraints, fabrication tolerances, optical losses, device mismatch, and large-scale routing complexities, might introduce deviations from idealized modeling. Some more in-depth discussion of those potential challenges, or references to relevant integration technologies, would strengthen the claims.\n\nAuthors’ response: It is true that our simulations did not explicitly account for fabrication tolerances. However, we expect the emission and absorption of spikes to be resistant against usual amounts of fixed-pattern noise. The former is explored extensively in [1], and the regions of parameter that provide excitable behavior are wide; fabrication procedures can be tuned accordingly. Spike reception is carried out by photodiodes, which are broadband and therefore resistant against spike wavelength variation.\nWhile it is true that certain implementations of the coupling matrix can have issues, this does not pose a significant problem to our general approach. Therefore, we view this particular discussion as lying outside of the scope of our article.\nFinally, in this work we assumed that the parameters of the system do not change over time. Therefore, the approach of calibration followed by training is sufficient. In practice, such parameter drifts are usually both weak and very slow, especially in comparison to the sampling time scales of PSN networks. Nevertheless, should this become necessary, it is straightforward to train or fine-tune the hardware system with the same wake-sleep algorithm as done in our simulations here. We expect such training to even improve the performance values reported here, as mentioned in the Discussion (line 621 ff.).\nYour question about scaling is particularly interesting and relevant. While extremely difficult to answer comprehensively, we can extrapolate from the current state of the art in optical circuitry. Our PSN networks assume the existence of an optical circuit providing an all-to-all connection among the nanolasers. Ideally, nanolaser could be connected with free space optics and programmable arrays of micro-mirrors and spatial phase modulators [2, 3]; in this setting, a large number of channels can be controlled independently with a manageable level of losses. Spiking nanolasers such as those reported in [4] could be modified for free space emission using compact grating couplers or free-form optics with very low losses (∼ 1 dB). On the other hand, MEMS technology offers a photonic integrated platform with ultra-low-loss switches connecting up to 240 ports [5]. The detailed breakdown of losses and required gain discussed in section “MNIST power consumption” is based on specifications given there. We note that the MEMS technology implemented there offers digital switching, which is unsuitable as such for our purposes. However, a continuous control of the weights is possible by using MEMS actuated directional couplers instead, with similar energy and loss specifications, as is the case for photonic matrices [6, 7], which feature 64 channels and are fully reconfigurable and low power. Finally, spectral multiplexing is possible owing to the incoherent nature of the interconnections, which could drastically increase the number of connections, as existing multiplexing/demultiplexing technology based on weighting filter banks handles several spectral channels [8, 9].\nA part of this answer has also been added to the end of the Discussion section (line 680 ff.).\n\nPlease address in the discussions points b, c, d and how these could be overcome and or limit the current work.\n\nPotential Flaws and points that require revision\n\na) Energy Estimation Details: The manuscript’s power estimates for large-scale networks could benefit from more specificity on amplifier or detector overheads, along with tolerance analysis of waveguide losses. As it stands, the rough estimates seem encouraging but may need tighter bounds or a more explicit breakdown.\n\nAuthors’ response: Regarding one specific detail of your comment, we note that nanoscale detectors have been shown not to need transimpedance amplifiers. In fact, their much reduced capacitance allows the use of large load resistance [10, 11]. Otherwise, we thank you for your suggestions and point to the more detailed power estimates provided now in the Discussion (line 655 ff.) and especially in the “MNIST power consumption” section of the Methods.\n\nb) Definition of Abbreviations: There are multiple instances where abbreviations appear before being defined. A comprehensive glossary or ensuring every acronym is defined upon first appearance would significantly improve clarity. (eg, SPS, PNS page 2 etc..)\n\nAuthors’ response: Already in the previous version, we tried to ensure that all acronyms are properly defined upon their first appearance, and every instance of every acronym also has a hyperlink pointing to its first appearance and definition. Upon rechecking, we found and amended a few missing hyperlinks and hope we haven’t missed anything this time. To improve readability, we have further eliminated some acronyms that only appeared infrequently in the manuscript (ANN, SBS and ISI).\n\nc) Methodological Completeness: While the theory and simulation approach are well described, readers from the photonic hardware community may want more details on device geometry or coupling schemes. Additional clarity about how well the authors expect the assumptions (e.g., linear summation of spikes, negligible inter-laser coherence) to hold in practice would help reproducibility.\n\nAuthors’ response: We have amended the Figure 1a and its caption to better represent the actual approach to coupling a PSN and a waveguide.\nThe nanolaser technology is described in detail in [12]. We provide this reference in the introduction (line 83 ff.). The spiking dynamics has been reported recently in a conference [4], and a manuscript is in preparation by our collaborators. We address this work in the introduction as well (line 89 ff.). Our model is based on parameters extracted from measurements of these devices.\nThe point you raise about coherence is indeed very important. Coherent coupling between lasers could happen due to residual (unwanted) reflection in the switching matrix, although in principle this should be very low. Yet, we believe this effect could be made negligible as nanolasers will be designed to have a large spread of the emitted wavelengths. It is perfectly possible to engineer laser individually on the same chip, with the nanolaser technology described in [12]. The emitted frequencies could ideally cover the full C-band (∼10 THz), while the spectral width of the emitted spikes is well below 100 GHz. This will also ensure linear summation of the spikes.\nA part of this answer has also been added to the end of the Discussion section (line 710).\n\nMethodology and Reproducibility\n\nAuthors’ response: We have added a “Simulation code details” section to the Methods clarifying the code structure and giving values of essential parameters.\n\nIn summary, the authors present a theoretically sound proposal. The authors make a strong case for how such ultrafast optical spiking neurons could perform sampling-based inference with notable improvements in speed and potentially power efficiency. Despite some points that would benefit from more details, particularly around practical power budgets and large-scale integration, the paper is overall a valuable contribution.\n\nAuthors’ response: Thank you again for your encouraging comments, helpful suggestions and constructive criticism.\n\nRecommendation: I recommend revisions to polish the manuscript, clarify energy estimates, ensure consistent definitions of abbreviations, and possibly strengthen the discussion around potential integration and scaling challenges.\n\nReviewer #2\n\nAuthors’ response: Thank you for your appreciation of our manuscipt. Please find detailed replies to your reservations below, which we hope will alleviate your remaining concerns.\n\n1) Broad audience and general interest. This manuscript is difficult to be understood by a non-specialist. It assumes a reader’s background on Bayesian statistics and Boltzmann machines, which is far to be the case even in the community of optical neuromorphic computing. It is then of primary importance to introduce the main concepts of Bayesian inference, Boltzmann machines and sampling (both non-spiking and spiking), with a basic example (prior to spiking models) which would help the reader to a large extent. Otherwise, we might consider this work to be better suited to a more specialized journal (JCP or CPC). In my opinion, introducing these concepts is substantially more important than the details of the mapping (“Membrane potential of spiking nanolasers”): most of this subsection could be moved to a supplementary material, retaining only the fundamental result in the main text.\n\nAuthors’ response: Accessibility is an important matter to us as well and we have thus taken your suggestions to heart. In particular, we have substantially extended our introduction to Boltzmann machines and sampling in the “Preliminaries” section (line 105 ff.) and to Bayesian statistics in the “Optical probabilistic inference” section (line 467 ff.).\n\n2) Overall clarity. Many points remain obscure to me (many of them in connection to point 1). For instance: it is not clear how the training is performed, i.e. what parameters (weights?) are tuned, and how. Another example is the 2D representations of network states (Fig. 7.a and 8.e): the vast majority of readers ignore what t-SNE is, so a hint to interpret the representation axes is needed. Furthermore, sentences like “...the PSN network correctly samples from the prior” (l 410) or “...by sampling from conditional probability distributions” (l. 345), or “the probability of each state sampled from the inverse...” (l. 378) are very technical and specific from Bayesian statistics, so they need further clarification here. How PSP is implemented is not clear neither.\n\nAuthors’ response: As for the point above, we have taken your concerns to heart. We have revised and expanded multiple explanations of our methods, both in the main text and in the Methods. Beyond the already existing “Probabilistic inference training” section, we have now added detailed descriptions of the following procedures to the Methods:\n• the contrastive divergence method used in training to sample from arbitrary distributions (section ”Contrastive divergence”, line 828 ff.);\n• the creation of the star plots used for visualization of probabilistic inference (section ”Probabilistic inference visualization”, line 867 ff.);\n• t-Stochastic Neighbor Embedding (t-SNE) used for visualization of the guided dreaming results (section ”t-Stochastic Neighbor Embedding”, line 1010 ff.).\nAs for the other points, we have revised the surrounding text; together with the additional information in the Methods (see above), we hope that these explanations are clearer now:\n• on “sampling from prior”, see line 491 ff.;\n• on “by sampling from conditional probability distributions”, see line 411 ff.;\n• on “sampled from the inverse...”, see line 445 ff.;\n• on the PSP: the explanation was given in the caption of Figure 3, which, we concur, might be missed by readers; we have thus added a short explanation and a reference in the main text (line 310 ff.);\n• for more clarifications and explanations, see also lines 344 ff. and 481 ff.\n\n3) The optical spiking advantage. It is claimed that the spiking dynamics is advantageous both in terms of energy saving and robustness to noise. While I acknowledge this as general statement, which may hold true in neuroscience, it is not clear to what extent it applies to the particular photonic system studied here. Specifically:\n\na. Power. CW laser-threshold powers are even smaller than those of pulsed devices.\n\nAuthors’ response: We agree with the argument that all else being equal, a saturable absorber increases losses and threshold pumping as a result. That said, we do not believe the argument is applicable in the context of this work: CW lasers would not be suitable for stochastic neural networks due to a lack of prominent stochasticity.\n\nb. Bandwidth. Frequency multiplexing is not clear in the context of the proposed photonic architecture, and should be commented.\n\nAuthors’ response: We have added an explanation after the multiplexing is mentioned (line 286 ff.).\n\nc. Losses. Interconnect matrix is based on waveguide crossbar arrays. Then the scaling of crossing losses is a potential issue in large networks, which is not discussed.\n\nAuthors’ response: We agree with this point, which was also brought up by the Reviewer #1. In general, losses can be compensated for by associating each port (each of them grouping lasers emitting in separated spectral channels) with a nanoscale semiconductor optical amplifier (i.e., based on the same technology) providing gain yet with very low power requirements (e.g., raising the average power emitted by the laser from 1 µW to 100 µW). For many further details, we kindly refer to our responses to Reviewer #1, sections ”Additional Evidence”, ”Energy Estimation Details” and ”Methodological Completeness”.\n\n4) The context of excitability in lasers. The references in the Introduction do not fairly reflect the historical developments. The first cited works on excitability (Refs. [23-25], Wünsch et al., Romeira et al., Prucnal et al.) do not accurately retrace the seminal studies on optical excitability. The first report on excitability was by F. Plaza et al. (EPL 38, 85, 1997) in a CO2 laser with SA, followed by experiments in semiconductor lasers (M. Giudici et al., Phys. Rev. E 55, 6414, 1997) and Nd-YAG lasers with SA (Larotonda et al., PRA 65, 033812, 2002). In the context of photonic crystals, the first demonstration in a nanocavity was reported in Brunstein et al, PRA 85, 031803 (2012).\n\nAuthors’ response: We thank you for sharing these references, which we have now added to the Introduction (line 75 ff.).\n\nMinor comments:\n\ni. Abstract. Saying that speed is an advantage of optics w.r.t electronics is not true in general, as an electronic signal also propagates at the speed of light in metallic wires. In addition, in the abstract (l. 5), what do authors mean by “conventional semiconductor devices”?\n\nAuthors’ response: We agree with this point, as “speed” is a rather unspecific term that may lead to confusion. Rather, we meant “bandwidth” and/or “latency”, that now replace “speed” where appropriate. As for “conventional semiconductor devices”: we meant “digital electronics”. See our amendments in lines 2 ff. and 5 ff. in the Abstract and 62 ff. and 65 ff. in the Introduction.\n\nii. Possible confusion in the notation: $\\beta$ is used as spont. emission factor of nanolasers, and as a parameter for probability distributions.\n\nAuthors’ response: Good catch! We have removed $\\beta$ as the spontaneous emission factor (line 204 ff.) as it is not used elsewhere.\n\niii. Absolute/relative refractoriness are assumed to be well-known concepts, yet they’re most likely ignored by the general reader.\n\nAuthors’ response: Agreed, we have added a brief explanation in the Discussion (line 605 ff.).\n\niv. Line 264 should read: “Second, changes of $n_e$ due to incoming spikes need to be small...”\n\nAuthors’ response: Good catch, thank you (line 330 ff.).\n\nReviewer #3\n\nAuthors’ response: Thank you for your very encouraging assessment of our work!\n\nHowever, a few points deserve clarification in my opinion:\n\n• The coupling between nodes is considered to take place via embedded photodetectors, aiming at incoherent transmission of spikes. Is coupling through saturation of the absorber section impossible or impractical?\n\nAuthors’ response: We believe it is possible, but impractical as it would require the coupling of a spike inside the receiving nanolaser. A sufficient coupling requires the alignment of resonances, which is challenging for a large number of nanolasers.\n\n• The model uses variabe S to describe a photon number, which must therefore be positive. What prevents S from becoming negative in presence of the stochastic term F? This is particularly relevant considering that most of the time the system is not emitting a spike thus in absence of noise S=0 (figure 1b).\n\nAuthors’ response: Indeed, noise can move the variables outside the “physical” range of values. During our simulations, S was forced to always be above a small positive number, and $n_e$ and $n_a$ to always be at least 0. However, this way we may have significantly altered the dynamics of the system. For this reason, we compared the continuous model to its discrete counterpart which does not suffer from such issues in the “Discrete nanolaser model” section of Methods, which now also includes this discussion (line 802 ff.). We found that both models give very close results.\n\n• Figure 2b shows time traces of the normalized gain, which is, as expected mostly slightly below 0. However, at the onset of spikes, the gain depletion is anticipated by a very large increase of the gain (particularly visible for instance on spikes 3,4,5), which takes place at a time scale which is comparable to the following depletion (which is mediated by the photon number growth). What is the origin of this very sharp increase? Could it be related to the previous point?\n\nAuthors’ response: The sharp increase is due to a positive feedback loop created by the saturable absorber (SA). The increase of the photon count due to the spontaneous emission reduces the electron density in the pumped section, which, in turn, reduces its gain. However, in our structure, the SA is a nanolaser itself, which shares the pool of photons with the pumped section. The SA will absorb the newborn photons to create free electrons. As the SA is not pumped, its differential gain is higher than in the pumped section, so the change of gain in the SA is positive and its magnitude is higher than the loss of gain in the pumped section. As a result, the overall gain of the structure has increased, which will create more photons, which will further increase the gain – a positive feedback loop. We observe this behaviour with the discrete model as well, therefore it is unrelated to noise-induced S < 0.\nAs for the “particularly visible ... on spikes 3,4,5” it is an unfortunate outcome of insufficient density of points in the figure: some gain peaks were lost between the subsequent points. We have increased the density of points around the spikes, and all spikes look similar now. We have also added a comparison of gain timetraces in the “Discrete nanolaser model” section (Figure 9) in the Methods.\n\n• The ”wake sleep” algorithm must be briefly described, with a reference for more information, if the work is to be readable also to non-experts in stochastic neural networks training.\n\nAuthors’ response: We absolutely agree. The explanation of the wake-sleep algorithm can be found in the “Probabilistic inference training” section in the Methods. However, we had unfortunately omitted a corresponding reference in the main text, which is now given at the first occurrence of the term (line 489 ff.).\n\nPoints to be considered: The “optical probabilistic inference” and “learning from data” sections seemed to me as “logical” specific applications of the section “optical sampling from arbitrary distribution”. If it is the case, I fear these sections do not really help the reader. I understand they constitute “bright” points and their implementation is probably worth a description but if they are really direct consequences of “optical sampling from arbitrary distribution” the reader would benefit from being informed of that very clearly.\n\nAuthors’ response: We have tried to better motivate our progression of functional applications the end of the “Membrane potential of spiking nanolasers” section (line 377 ff.).\n\nMinor points:\n\n• figure 1a: is a spike really shorter than the laser resonator as suggested by fig 1a? If yes, a PDE system would really be required to describe the physics. If not, then a single mode model as used here is fine but then figure 1a is a bit misleading.\n\nAuthors’ response: This is quite a subtle point, you are absolutely right and thank you for pointing it out. The spikes are indeed actually “long”, as their FWHM of ∼30 ps corresponds to ∼3 mm, whereas a typical photonic crystal is roughly 10 µm long. We have therefore removed the spikes from the Figure 1a, as they were also not very instructive at that point in the first place.\n\n• ”alpha shaped” is used here and there in the manuscript but not defined\n\nAuthors’ response: We have added a definition at the first mention (line 313 ff.).\n\n• The connection between BMs and spiking networks cites 11-13, I believe Hinton, G. E., & Brown, A. (1999). Spiking boltzmann machines. Advances in neural information processing systems, 12. should be credited.\n\nAuthors’ response: Thank you for bringing our attention to this article. While the proposed implementation is quite different from ours, we were very happy to include its citation in the manuscript (line 46 ff.)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for thoroughly addressing all of the substantive points raised: they have strengthened the introduction by emphasizing the sampling niche with new citations; provided a detailed “MNIST power consumption” section to clarify energy budgets; discussed fabrication tolerances and routing challenges by covering noise tolerance, photodiode resilience, and options such as free-space/MEMS/spectral multiplexing; cleaned up infrequent acronyms, fixed hyperlink definitions, and added a “Simulation code details” section with annotated parameters to improve methodological clarity and reproducibility.\nA single follow-up question remain: do the authors intend to publicly release their simulation code, and how should readers interpret the statement “Simulation code is available upon a reasonable request”? Could the authors clarify under what conditions and through what process one might request and obtain the code?", "reviewer_2": "Thank you for having modified the manuscript according to my comments and questions. Most of the “obscure points” have now been clarified in the text. Still there are a couple of issues that need to be addressed:\n1) The “optical spiking advantage”. The authors’ response to my point 3a) is: “we do not believe the argument is applicable in the context of this work: CW lasers would not be suitable for stochastic neural networks due to a lack of prominent stochasticity.”. This is not true for a nanolaser: close enough to laser threshold, the amplified spontaneous emission noise dominates the dynamics. Therefore, the question remains open: are spikes in optical NNs advantageous enough –as they are in the brain–, justifying its implementation over other (stochastic) approaches in a photonic context? I think a short discussion on this point in needed.\n2) The authors have revised their use of the word “speed”, and now focus on \"bandwidth\" and/or \"latency\". Unfortunately, their meaning in the present context, together with their advantages, remain unclear:\n- The “bandwidth” advantage remains non-justified in the present context. In particular, the paragraph from lines 279 to 298 (of the redlined document) is confusing. Does wavelength multiplexing present an advantage in terms of parallelism? Or is it just a way of selectively tuning the weights (Fig. 3d)?. By the way, Fig. 3d has been added in the new version, but no explanation is provided in the main text. Is this implementation instrumental for the device operation? I observe that, adding microrings to the phtotonic chip makes the architecture even more complex: hundreds of nanolasers, photonic cross-bars, amplifiers, balanced-detectors, and now even microrings… A realistic implementation of this can be extremely challenging, and even criticized. In my opinion, the authors should keep their designs realistic, at least they should clarify what are the fundamental ingredients needed for a proof-of-principle experimental demonstration.\n- What do the authors mean by “latency”?. Do they mean that extra delays in signal processing which could be alleviated in the photonic implementation? I note that the use of “latency” in the context of excitability might be confusing, because refractory times is a form of latency. The authors should clearly state what they mean by “latency” in this context and in what sense the nanolaser PSN are beneficial.\nAdditional remarks:\n- Equation (3) is not written correctly: since the Wiener process dW_t is a differential quantity, the other terms should also be expressed as differentials (as in Equation (6)).", "reviewer_3": "The authors have adequately addressed all my comments.\nOne clarifying remark: In their reply about the \"gain\" becoming positive in Figure 2b, the authors discuss a \"net gain\", which includes the role of the amplifying and the absorbing section. That is now very clear to me but was initially not, perhaps for good reason: lines 225 and following state \"Stimulated emission and\nabsorption are described by a term GS, where G = γr (2ne − n0 ) represents the gain, n0 is the total number of dipoles, and γr\nis the radiative transition rate\". This is the \"gain\" I had in mind while parsing figure 2b. Instead, in their discussion, the\nauthors refer to another \"gain\", defined line 243 as \"G = γr (2na − n0,a ) + γr,a (2na − n0,a ) − γ is the net gain including photon\ndamping\". Just adding \"net gain\" instead of \"gain\" where applicable (in particular caption of fig 2) and perhaps choosing\nanother form of \"G\" for the one defined line 225 and the one defined line 243 would remove the possibility of\nmisunderstanding."}, "author_response": "Reviewer #1\n\nAuthors’ response: While the manuscript offers a full and detailed description of the performed simulations, we are open to receiving individual requests for direct access to the code. We will meet any such request with our full support, but our institutional regulations require us to obtain case-by-case approval by the management of Thales.\n\nReviewer #2\n\nAuthors’ response: We agree that our previous response about CW lasers not being suitable for this task was not nuanced enough. We definitely do not mean to claim that our specific implementation of Boltzmann sampling is the only way to perform stochastic Bayesian computation. The important point here is that the solution presented in our manuscript represents a near-perfect match for neural sampling as a demonstrated means of carrying out sampling from Boltzmann distributions. This framework relies crucially on stochastic, pulsed neuronal responses followed by pronounced refractory periods. We have amended the introduction (see line 92 ff.) to hopefully clarify this point better.\n\n2. The authors have revised their use of the word “speed”, and now focus on “bandwidth” and/or “latency”. Unfortunately, their meaning in the present context, together with their advantages, remain unclear:\n\n– The “bandwidth” advantage remains non-justified in the present context. In particular, the paragraph from lines 279 to 298 (of the redlined document) is confusing. Does wavelength multiplexing present an advantage in terms of parallelism? Or is it just a way of selectively tuning the weights (Fig. 3d)?. By the way, Fig. 3d has been added in the new version, but no explanation is provided in the main text. Is this implementation instrumental for the device operation? I observe that, adding microrings to the phtotonic chip makes the architecture even more complex: hundreds of nanolasers, photonic cross-bars, amplifiers, balanced-detectors, and now even microrings… A realistic implementation of this can be extremely challenging, and even criticized. In my opinion, the authors should keep their designs realistic, at least they should clarify what are the fundamental ingredients needed for a proof-of-principle experimental demonstration.\n\nAuthors’ response: We believe the bandwidth advantage is justified, but not necessarily in the sense of handling faster inputs, but facilitating data transfer between internal components. The broad optical bandwidth allows for wavelength multiplexing which makes the coupling between nanolasers more compact without losing data transfer parallelism. For example, in digital spiking networks, several neurons share the same bus; each neuron must wait until the bus is free to perform I/O. In optics, one waveguide can transfer spikes from multiple nanolasers, provided their resonances are separated.\nWavelength multiplexing is not required for tuning the weights of nanolaser coupling. If spikes from different nanolasers interfere (which is a challenge in itself), one might consider a matrix-vector multiplication circuit such as in [Shen et al., Nat. Photon. (2017)].\nWe agree that the technology is extremely challenging, although integrated photonics is improving steadily and in the examples shown in the cited references, photonic chips with several hundreds of functional components exist already. We do not claim that our approach is “better” than others and the experimental implementation we suggest is by no means expected to be optimal, rather, we regard it as an educated guess of a feasible architecture. We have now addressed this point explicitly on line 266 ff. and added references to Figures 3c,d in the discussions that follow.\n\n– What do the authors mean by “latency”?. Do they mean that extra delays in signal processing which could be alleviated in the photonic implementation? I note that the use of “latency” in the context of excitability might be confusing, because refractory times is a form of latency. The authors should clearly state what they mean by “latency” in this context and in what sense the nanolaser PSN are beneficial.\n\nAuthors’ response: Indeed, latency can be interpreted differently depending on the context. In our text, latency only appears in the introduction, in the context of computing speed, as also addressed by the provided references [23-26]. We modified the sentence immediately following the first appearance of the term (line 61 ff.) to make this even clearer. In other words, “low latency” can also be translated to “short time-to-solution”. As we show later on, the sampling speeds achieved by our proposed methods are orders of magnitude higher than those of state-of-the-art electronic neuromorphic systems, translating to a corresponding decrease in the time required to obtain an approximation of the sought probability distribution with some desired degree of precision.\n\nAdditional remarks:\n\n– Equation (3) is not written correctly: since the Wiener process dW_t is a differential quantity, the other terms should also be expressed as differentials (as in Equation (6)).\n\nAuthors’ response: Good catch, thank you! We have amended the equation.\n\nReviewer #3\n\nAuthors’ response: That is true, thank you! We now refer to G = γr (2ne − n0 ) + γr,a (2na − n0,a ) − γ as “gain”. Please see lines 217 and 236 ff."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The manuscript proposes and analyzes the concept of spiking nanolasers, i.e. semicon- ductorlasersthatemitultrafastopticalpulses, anddemonstratesamodelforcomputing based on sampling-based neural computation. The core idea is to use small photonic crystal nanolasers, each behaving as a stochastic spiking neuron, to perform Bayesian inference through sampling from learned probability distributions (e.g., Boltzmann ma- chines). The authors provide a theoretical treatment comparing these nanolaser-based neurons to conventional leaky integrate-and-fire sampling (LIFS) neurons, show how to map typical Boltzmann machine parameters onto spiking nanolasers, and then illustrate the approach on generative tasks and pattern completion (e.g., small-scale MNIST).\nKey contributions include:\n• A derivation that links nanolaser gain dynamics and standard stochastic neuron models.\n• Numerical experiments suggesting significant speed and potential power-efficiency gains over purely electronic spiking systems.\n• Demonstrations of sampling accuracy for small toy problems (MNIST-like tasks).\nOverall, the paper provides an interesting proposal for ultrafast, optically based neuro- morphic systems. The results are timely, given the growing interest in photonic comput- ing for machine learning, and the manuscript’s focus on sampling-based models (rather than purely feedforward neural nets) highlights unique advantages of spiking neural sub- strates.", "reviewer_2": "Thankyouforhavingmodifiedthemanuscriptaccordingtomycommentsandquestions. Most of the “obscure points” have now been clarified in the text. Still there are a couple of issues that need to be addressed:\n1. The “optical spiking advantage”. The authors’ response to my point 3a) is: “we do not believe the argument is applicable in the context of this work: CW lasers would not be suitable for stochastic neural networks due to a lack of prominent stochasticity.”. This is not true for a nanolaser: close enough to laser threshold, the amplified spontaneous emission noise dominates the dynamics. Therefore, the question remains open: are spikes in optical NNs advantageous enough –as they are in the brain–, justifying its implementation over other (stochastic) approaches in a photonic context? I think a short discussion on this point in needed.", "reviewer_3": "The authors have adequately addressed all my comments.\nOne clarifying remark: In their reply about the “gain” becoming positive in Figure 2b, the authors discuss a “net gain”, which includes the role of the amplifying and the absorbing section. That is now very clear to me but was initially not, perhaps for good reason: lines 225 and following state “Stimulated emission and absorption are described by a term GS, where G = γr(2n −n ) represents the gain, n is the total number of e 0 0 dipoles, and γ is the radiative transition rate”. This is the “gain” I had in mind while r parsingfigure2b. Instead,intheirdiscussion,theauthorsrefertoanother“gain”,defined line 243 as “G = γ (2n −n )+γ (2n −n )−γ is the net gain including photon r e 0 r,a a 0,a damping”. Just adding “net gain” instead of “gain” where applicable (in particular caption of fig 2) and perhaps choosing another form of “G” for the one defined line 225 and the one defined line 243 would remove the possibility of misunderstanding."}, "author_response": ""}]

Decision letter