Designing Polaritonic Integrated Circuits for Quantum Processing

Venue: NeurIPS_cc_2024_Workshop_MLNCP
Published: 2024-09-12 12:10:02

Abstract: We propose photonic integrated circuits augmented with a semiconductor exciton-polariton nonlinearity to accomplish two fundamental tasks in quantum processing: quantum state tomography and quantum state generation. We demonstrate in simulations that the design of the circuit can be optimized to great effect and showcase the efficacy of the optimized polaritonic circuits for the cases of i) fully identifying a family of emblematic quantum states and ii) stabilizing an accurate train of single photons.

Main paper: 
 introduction: 
 the heart of quantum optics [3] , thus paving the road towards large-scale optical quantum computing [4, 5] . While in the classical approach to optical computing [6, 7] information is captured solely in the relatively easily measurable and modulatable amplitude and/or phase of different components of light, a quantum state of light is fully determined by (i) full-counting photon statistics and (ii) the relative phase of each photon number outcome -features of fundamental importance for quantifying quantum entanglement and state superposition [8] .Strongly coupled matter-light systems, such as exciton-polaritons [9, 10] , have attracted growing interest in the quantum-computing community. Thanks to the flexible state manipulation provided by the incident laser field and direct access to the emitted light, experiments with exciton-polaritons have accomplished a plethora of impressive results -e.g., the creation of spontaneous long-range coherence (Bose-Einstein condensation) [11, 12] , realization of a superfluid flow [13, 14] , or the stabilization of an optical topologically protected insulator [15, 16] . It is, therefore, not surprising that recently these platforms have also been proposed for optical quantum computing [17] [18] [19] .In this paper, we introduce a photonic integrated circuit (PIC), built upon a platform of nonlinear photonics, that creates a novel pathway for achieving universal quantum state generation and characterization. See Fig. 1 for a schematic depiction.Recently, a series of papers showcased the potential of utilizing quantum reservoirs for the tasks of state tomography -an all-encompassing characterization of a photonic quantum state [20] [21] [22] [23] -and state generation [24] . Although the presented results were innovative and notable from a conceptional point of view, limited efforts were dedicated to the actual implementation, the design, or the potential deployment in state-of-the-art experimental settings. The primary goal of this paper is to do precisely that.The realization of polaritonics in integrated photonics [25] served as a motivational cornerstone for this work. We demonstrate a feasible implementation of both quantum tomography and quantum state generation using polaritonic PICs. Moreover, while reservoir computing is, by definition, a non-tunable technique to enlarge the feature space, we demonstrate the use of machine-learning methods for optical circuit design to shape an optimal feature space for the task at hand.Our contributions can be summarized as follows:1. we showcase the potential of integrated nonlinear photonics to perform quantum state tomography and quantum state generation 2. we provide evidence that one can efficiently and to great effect optimize the design of the polaritonic PIC for tomographical tasks 3. we show the efficacy of our polaritonic device for the task of photonic state generation. Specifically, we focus on the technologically relevant case of generating highly accurate trains of single photons, at maximal emission rate. 
polaritonic photonic integrated circuit optimization: 
 On-chip quantum dynamics. The full quantum dynamics, describing photons traversing the chip of interconnected nonlinear waveguides illustrated in Fig. 1 , is governed by the interplay of energy-conserving Hamiltonian dynamics and dissipation in the form of, primarily, polariton losses.Under the Born-Markov approximation, this is fully captured by the Lindblad master equation for the photonic density matrix -see, e.g., Ref. [26] (ℏ = v g = 1, with ℏ Planck's constant and v g the polariton group velocity, see Appendix A for information on units),EQUATIONThe first term describes the coherent dynamics, defined by the HamiltonianEQUATIONand the second the photonic losses. Here, a i and a † i are the photon creation and annihilation operators, respectively, which create or remove a photonic quantum from a quantum state |ψ⟩, and L is the number of waveguides. Let us clarify the parameters that were used in Eqs. ( 1) and (2) . Some of them are semiconductor material properties and are, therefore, not easily accessible for tuning. These include, the detuning ∆ of the incoming light with respect to the exact polariton resonance -in what follows, we assume exact resonance, so ∆ = 0; the polariton loss rate γ; and the strength of the nonlinear photon interaction U that is mediated by excitons, which are dipolar semiconductor excitations.Of fundamental importance is the dimensionless ratio U/γ, quantifying the effective nonlinearity experienced by the photon field. While in traditional microcavity platforms a modest optical nonlinearity has been reported [27] , it is believed that much larger effects can be engineered in integrated waveguide designs [25] .Besides the fixed properties of the material stack ∆, U, and γ, the set of parameters J ij (t) in (2) remains directly accessible for the purpose of optimizing the circuit design. These parameters describe the graph matrix of waveguide couplings, as fabricated using directional couplers (DCs) or multi-mode interferometers (MMIs), in combination with phase shifters (PSs). We choose to keep the matrices J ij real and will always choose a brick-like pattern of couplings, see Fig. 1 .In this work, it is precisely the flexibility in the design of the PIC, parametrized by J ij (t) in ( 2), that we exploit to optimize the performance of the chip on two exemplary tasks: (i) the quantum tomography of incoming Gaussian squeezed states of light and (ii) the emission of a single-photon train at maximal emission rate.Quantum simulation. In principle, full exact quantum simulation of the dynamics governed by Eqs.(1) and ( 2) can be achieved using the toolbox of linear algebra, by representing quantum states |ψ⟩ as vectors and the operators as matrices. Noteworthy, even though the method is, in principle, exact, the vector space dimension grows exponentially as D = N L , for L waveguides carrying photon states up to occupation number N . Needless to say, this relates back to the fundamental hardness of classically simulating quantum mechanics. Also worth noting is that for photonic systems, the occupation number can be any positive integer. That's why we rely on computationally more efficient techniques with manageable memory overhead, to wit, a Gaussian variational approximation and the method of Monte Carlo unravelings of the master equation. For details on these simulation techniques we refer the reader to the supplementary material, see sections B and C.Quantum state tomography and generation. For quantum state tomography, the photonic integrated circuit establishes a quantum extreme learning machine [28] . Recall that a classical extreme learning machine pre-processes data via a memoryless (i.e., non-recurrent) non-linear map, while the final output is obtained via a trainable linear layer. In its quantum version, the preprocessing of the input data -in our case a quantum state of light -is governed by a quantum-mechanical time-evolution, followed by the quantum measurement of a collection of observable quantities. In our realization, the quantum dynamics is governed by the above-described effects, while we measure the expected photon count at the end of each waveguide.• Polariton-mediated interaction As these coefficients govern the quantum dynamics as in (2), they appear as parameters in the ordinary differential equations describing the time evolution of the quantum state. Choosing a gradient descent approach to perform the optimization, we can rely on the adjoint state (or sensitivity) method to find the updates to the ODE parameters. Within the machine learning community, this approach was reintroduced in [29] , and we leverage the associated PyTorch package. Note that the linear layer is always fitted perfectly optimally given the quantum measurements of the occupation numbers.For state generation, we optimize the parameters of the unitary mixing, which we parameterize asV ( ⃗ θ) := exp iG( ⃗ θ) .Here, G is an L × L Hermitian matrix, representing the inter-waveguide couplings that we want to optimize to obtain stable high-probability single-photon outcomes -this is perfectly similar to the optimization of the coefficients J ij (t) for tomography. For this usecase, the nonlinear PIC is kept random and fixed. The rationale is that we ultimately aim for a device that is universal in its quantum state generation capabilities. There is, however, no obstruction to designing a device tailored to generating a specific state by leveraging the techniques of the previous paragraph. 
results: 
 In this section, we present our results, based on numerical simulations, on optimizing the tasks of quantum tomography and state generation. See Appendix E for the simulation details.Gaussian quantum states are multi-mode optical states characterized by the property that number statistics -as expressed by expectation values of strings of creation and annihilation operators -are governed by Gaussian statistics. Amongst this class of states, we focus on single-and two-mode thermal squeezed states, often abbreviated as SMTSS and TMTSS respectively. As squeezed versions of thermal states, they are specified by three real parameters: the equal thermal average number of photons n in each mode and the complex squeezing parameter ξ = re iθ . These parameters fully determine the covarience matrix of the above-mentioned Gaussian statistics. We do not displace the state, i.e., we keep the distribution's first moments zero. In formulae, the tomographical information of an SMTSS can be summarized by the following expectation values [3] :⟨a⟩ = 0 , n = ⟨a † a⟩ = n + 1 2 cosh(2r) - 1 2 , c = ⟨aa⟩ = -n + 1 2 sinh(2r)e iθ .(3) For two-mode thermal squeezed states similar expressions determine the equal diagonal elements of N ij = ⟨a † i a j ⟩, while the off-diagonal entries are zero, and the equal off-diagonal elements ofC ij = ⟨a i a j ⟩, whose diagonal is zero.Naturally, the indices i, j run over the values 1 and 2. We implement the quantum state characterization of thermal squeezed states as the task of determining n and the real and imaginary parts of c.Concretely, we inject the single/two-mode thermal squeezed state into one/two input waveguides of the chip of Fig. 1 . We choose the chip to have five waveguides, and send fixed coherent laser light, described quantum mechanically by coherent states, into the remaining four/three waveguides.We consider five layers of interference regions in a brick-pattern. We created a training and test dataset, each consisting of 250 random samples. We measure the total mean squared error (MSE) loss after the linear layer. To efficiently optimize the MSE during training, we implement a two-stage optimization. First, we optimize the linear layer given the features provided by the quantum circuit.As this is standard linear regression, it is an easily solved convex optimization problem that, in fact, has an analytical solution. Next, we freeze the weights of the linear layer and optimize the coupling parameters of the quantum circuit, simulated using the Gaussian variational approximation, by taking a gradient descent step. In practice, we use SciPy's L-BFGS-B algorithm to perform the gradient descent.In Fig. 3 , we summarize our results. Panel (a) contains a boxplot showing statistical information on the improvement in the test mean squared error that can be achieved by optimizing a randomly initialized quantum circuit; we do so as a function of the nonlinear interaction strength. We start by noting that a random circuit is the waveguide equivalent of the currently state-of-the-art quantum reservoirs of [20] [21] [22] [23] . In particular we observe a more than ten thousand-fold improvement in median MSE for the currently experimentally feasible interaction strength of U/γ = 0.1 [27] . The necessity of the photonic nonlinearity is also made manifest by the very poor results in the bin log 10 U/γ = -∞, i.e., U = 0. In panel (b), we showcase an example of a circuit for SMTSS characterization at U/γ = 0.1. Finally, in panel (c), we show several circuit optimization curves, initialized with random circuits, as measured by the test MSE on the two-mode characterization task. We easily achieve an order of magnitude improvement. As we didn't change the overall architecture of the PIC, it is natural that the MSE of the two-mode characterization cannot fully match the single-mode MSE. Indeed, two-mode characterization requires the system to disentangle the statistical cross-correlation between two strongly entangled modes, instead of merely identifying single-mode, uncorrelated, statistics. Nevertheless, we obtain perfectly satisfactory values for the MSE.Our setup for quantum state generation is illustrated in Fig. Density-density correlations. The correlation function g (2) (τ ) defines the density-density correlations between a pair of quantum modes i and j. It is defined as, gEQUATIONWe focus on single-photon generation in one output arm l (i.e., i = j = l), at zero time delay (τ = 0). It can be shown that the quantity gll (0) is always positive. Moreover, it is only zero for single-photon states. In fact, minimizing gll (0; ⃗ θ), where θ i are the PIC configuration variables, amounts to maximizing the single-photon generation efficiency [30] [31] [32] .Density matrix. We also formulate a FOM in terms of the entries of the reduced density matrix ρl of output arm l,EQUATIONNow, the goal is to maximize the matrix entry ρl,11 , the single-photon probability, or otherwise filling the vacuum mode ρl,00 (0 ≤ w 0 ≤ 1). All other matrix entries ρl,00 are suppressed by virtue of the last term. In (5), we explicitly aim at maximizing the single-photon probability, which was absent in g (2) (0). Indeed, in (4), very accurate single-photons can be generated, at the cost of obtaining a very low intensity n l .Nevertheless, measuring gll (τ ) is a standard task, accomplished in a Hanbury Brown -Twiss interferometric setup (e.g., Ref. [27, 33] ), opening the possibility of real-time PIC optimization based on experimental outcomes. In stark contrast, experimentally obtaining matrix entries of ρ l for (5) is much more complicated [34] . for g (2) (0; V opt ), this comes at the expense of low single-photon output intensities and thus low probabilities -only a few percent chance of single-photon emission, see yellow dots.On the right-hand side of Fig. 4 , we show the results for using the objective F(V ) of Expr. 5. These results were obtained for the same nonlinear PIC configuration, but allowing for larger interaction coefficients U . The input intensity is set to α = 1. The curve in panel (d) shows some wiggles, but we observe that generally better values of the FOM are obtained at higher U/γ. In panel (e), the obtained optimal density matrix for U = 35γ is shown and it is clear that, up to a good approximation, only the matrix entries ρ 11 and ρ 00 are occupied. Finally, in panel (c), we compare the probability of producing vacuum and the probability of single-photon generation. We find a single-photon probability of about p 1 = 0.4 for large values of U/γ. Clearly, higher photon numbers (inducing spurious signals ) are strongly suppressed, since p 0 + p 1 ≈ 1 for all U -values (black dotted line). 
conclusions and outlook: 
 We have proposed a novel polaritonic PIC to perform universal quantum state characterization and generation. Moreover, we introduced and successfully brought to bear a methodology to optimize the circuit so as to maximize the efficiency of the respective tasks. We provided evidence, based on numerical simulations, of the efficacy of the optimized chips in the tasks of characterizing single-and two-mode squeezed thermal states and the highly accurate generation of single photons.We are currently actively corresponding with experimental research groups about realizing our proposed PICs, using GaAs-based polariton quantum wells. Once implemented, the PIC will be tested and benchmarked against our numerical simulations. Furthermore, for both state tomography and generation, we aim at scaling up our simulations to larger system sizes, using matrix-product based methods [35] , to the end of characterizing and generating states of higher complexity -e.g., for performing full density-matrix tomography or generating Bell-pair states. These results will be benchmarked and compared against the original (conceptual) results of Refs. [20, 24] . 
 
 Appendix: In Eqs. ( 1) and ( 2) we briefly indicated that we choose to set ℏ = v g = 1, with ℏ the reduced Planck constant and v g the polariton propagation group velocity. Here we give an overview of the motivation and the conversion to actual experimental units.A photon with frequency ω has an energy equal to E = ℏω. Thus, by choosing ℏ = 1, energy and frequency scales have, by construction, the same units. One can similarly argue for the effect of setting the polariton group velocity to one by considering d = v g t, the distance a polariton propagates, in a time t, through the waveguide after injection. The value for v g can be determined from the waveguide polariton dispersion ϵ(k), v g = ∂ϵ(k) ∂k . Note that it is strongly dependent on the angle of laser incidence and the Rabi coupling, the coupling between the photon and exciton field -see Ref.[9].Unifying the dimensions of time and length, using v g as the conversion factor, facilitates the simulations but has, on the other hand, a few important assumptions; all photon backscattering, due to roughness in the waveguides, is neglected, as well as potential inter-waveguide dephasing effects, caused by small length mismatches -see the paraxial approximation [36] .Our simulations are generally performed with polariton decay rate set to unity, γ = 1. Consequently, experimentally obtaining a value for polariton lifetime of, e.g., τ = γ -1 = 10 ps (ps: picoseconds) and for group velocity v g = 20 µm • ps -1 allows for the conversion to actual energy, frequency and length units, using ℏ = 0.658 meV • ps (meV: mili-electronvolt). As an example, a value U γ = 0.1 from the main text, would, in this case, correspond to U = 0.1 • ℏγ = 66µeV and a circuit length ofT = 0.5 γ -1 to T = 0.5 • v g τ = 200µm.Simulating the dynamics of a quantum system is achieved by integrating the Schrödinger equation,EQUATIONHere, |ψ(t)⟩ is conveniently represented as the quantum-state vector, |ψ⟩ = ⃗ ψ, where each vector element ψ ⃗ n is a complex numbers (amplitude and phase). Evaluating p ⃗ n = |ψ ⃗ n | 2 gives the probability of detecting the system in a number state ⃗ n, upon performing a photon-resolved measurement. The number states are labeled as ⃗ n = n 1 , n 2 , . . . , n L T , with n i the number of photons in waveguide i and we consider L waveguides in total. The core obstacle of quantum simulation on classical hardware is the exponential increase of the vector space (Hilbert space). Indeed, suppose we want to track up to N photons in one channel, and we have L channels in total, the full vector-space dimension is given by D = N L .The Hamiltonian H(t) appearing in Eq. 6 is the energy operator applied to the quantum state |ψ⟩. It is represented as a D × D matrix, acting on the state vector ψ ⃗ n . For the photons traversing a set of L coupled nonlinear waveguides, the space axis of photon propagation can be mapped to an effective time dimension, for which the Hamiltonian is given by,EQUATIONHere, a i and a † i are, respectively, the bosonic annihilation and creation operators. The former annihilates one photon from waveguide i with n i photons,a i |n i ⟩ = √ n i |n i -1⟩, the latter injects a photon, a † i |n i ⟩ = √ n i + 1|n i + 1⟩, and their combination gives the photon number,a † i a i |n i ⟩ = n i |n i ⟩.The parameter ∆ = ω L -ω LP is the detuning between frequency of the incoming light ω L and the lower-polariton resonance ω LP . The matrix J(t) is the graph matrix of waveguide couplings, so each matrix entry J ij (t) describes the tunneling amplitude of one photon from waveguide i to waveguide j. The parameters J ij (t) are kept explicitly time-dependent, since, by design, couplings between waveguides change as photons cross different layers imprinted on the chip. Finally, the parameter U quantifies the photonic nonlinearity.Integrating the Schrödinger equation ( 6) assumes coherent energy-conserving dynamics in a perfectly isolated environment, without any effects of dissipation or noise. However, in polaritonic systems, polariton losses (due to exciton non-radiative decay, photon scattering...) are unavoidable and must be included in the dynamics. This is described by the the Lindblad Master Equation [26] ,EQUATIONin which the first term describes the coherent unitary dynamics, analogous to Eq. ( 6), and the second the photonic losses in the waveguides (operators a i ) with a (uniform) rate γ. Often this is described by the Liouvillian L, a superoperator of which the spectrum shows interesting details about the dissipative dynamics and convergence to the steady state [37] .Crucially, now the quantum simulation must be evaluated on the full density matrix ρ, which is a D × D matrix that can be composed in a set of pure states,EQUATIONwith p k the probability of encountering the pure quantum state k. Importantly, the representation in terms of a quantum state ensemble |ψ k ⟩ is not unique and this relates back to a fundamental invariance of the Lindblad master equation ( 8) [38, 39] .Fortunately, the full simulation of the master equation can be circumvented by using a Monte Carlo (MC) sampling technique to collect a set of k pure quantum states|ψ k ⟩ (so D-dimensional, not D 2 ),which approximately describe the density matrix in some unraveling,EQUATIONImportantly, the density must never be explicitly evaluated from the sampled states |ψ k ⟩ if only expectation values of operators are important (that is, measurement outcomes),EQUATIONThe method was originally introduced, independently, in Refs. [40, 41] , see Ref. [26] for a good overview. We resorted to the MC quantum trajectory technique for the quantum simulation needed for state generation in Sec. 3.2.Numerically, the quantum simulations were performed using the python package qutip [42] .Rather than performing a full quantum simulation of the Lindblad master equation ( 8), the dynamics can be captured using the ansatz that the quantum states attain Gaussian statistics.In quantum mechanics, a Gaussian quantum state obeys Gaussian number statistics and is, as a direct consequence, fully defined by the following quantities:• Operator expectation value: α i := ⟨a i ⟩,• Normal operator correlations:N ij := ⟨δa † i δa j ⟩,• Anomalous operator correlations:C ij := ⟨δa i δa j ⟩.Above, we defined the fluctuation operator δa i := a i -⟨a i ⟩.A coupled set of closed-form nonlinear ordinary differential equations (ODEs) can be derived to describe the lossy circuit dynamics, within the approximation of the photons obeying Gaussian number statistics -the so-called Hartree-Fock-Bogoliubov method (see, e.g., Ref. [43] for a clean derivation),EQUATION)EQUATIONEQUATIONHere, we defined n i := N ii and c i := C ii and only nearest-neighbor couplings J i,i±1 are considered.Using standard methods (e.g., Runge-Kutta-45), the system of equations ( 12)-( 14) can be integrated to simulate the Gaussian quantum state of photons emitted by the coupled nonlinear waveguides.Evaluating the waveguide intensities at time t is easily done by reading outI i (t) := |α i (t)| 2 + n i (t).Notably, assuming Gaussian statistics of the photonic quantum state greatly reduces the complexity of quantum simulation. Indeed, instead of requiring the quantum simulation of the Schrödinger equation (coherent) (6) or Lindblad master equation ( 8) in the full exponential Hilbert space with dimension D = N L , we can now resort to a set of coupled nonlinear ODEs ( 12)-( 14) -polynomial scaling, not exponential. Of course, this gain in efficiency comes at the cost of omitting part of the quantum correlations that builds up during quantum time evolution. For short enough times and sufficiently low nonlinearity U , however, we see the results can be reliable [43] .The task of designing a PIC to obtain an optimal train of single photons as output consists of two parts, as was depicted in the main text in Fig. 2 . First, coherent laser light is coupled into a first non-tunable PIC, with a polaritonic nonlinearity, to generate a multimode entangled output state.This generated state is next injected in a fully linear PIC for the single-photon optimization task. We describe in detail the simulations below.The input laser light is defined as a coherent state |β⟩, which have as defining property that they are eigenstates of the bosonic annihiliation operator, a|β⟩ = β|β⟩. As a direct consequence, coherent states show Poissonian number statistics, ⟨a † a⟩ = Var(a † a) = |β| 2 , with n = ⟨a † a⟩ the number expectation value. A single-frequency continuous wave (CW) or, equivalently, a laser beam is described by a coherent state |β⟩, where the complex number β = |β|e iθ defines the amplitude and phase of the laser.The output light of the first circuit is represented by the density matrix, expressed using the sampled MC trajectory states,EQUATION)For a given optical quantum circuit, characterized by the interaction strength U , the dissipation rate γ and a set of waveguide couplings J ij (see Eqs. ( 2) and (1)), a set of sampled quantum states |ψ s ⟩ was stored on local hardware for the ensuing task of optimizing the linear circuit.The generated set of MC quantum states |ψ k ⟩, describing the density matrix ρ out , are now injected in a second, fully linear PIC. We aim at maximizing the probability of generating a single-photon as output in one of the output channels of this second PIC. Given that (i) the PIC is now fully linear and(ii) we neglect losses during optimization, we can circumvent a cumbersome full quantum simulation of its dynamics and resort to an equivalent operator formalism -see [44] . The output modes of the PIC are now a linear transform of the original bosonic operators,EQUATIONwith V the unitary matrix to be optimized and representing the PIC. V is of size L × L (so quadratic in the number of waveguides, not exponential). We define it as the matrix exponent of a Hermitian matrix G that represents the graph of waveguide connections,EQUATIONSince the only restriction to G is it being Hermitian, this guarantees a straightforward and efficient expression of V in terms of optimization variables θ i .Expressing the optimization objective in terms of the transformed operators b i (V ) from Eq. ( 17) is essential. For single-photon generation, we present two possible objectives,The density-density correlations. The coincidence rate of photons is quantified by the g (2) (0) function, defined as, gEQUATIONwhere ⟨•⟩ is the expectation value of a quantum state ρ, here evaluated over the MC sampled trajectory states (see Eq. ( 11)). We optimize for one output mode l. Essentially, (19) represents the ratio of fourth-order correlations (∼ density variation) over density squared. For a coherent laser beam, thanks to the Poissonian number statistics, gl (V ) = 1, regardless of the choice of V . In contrast, gl (V ) = 0 is the correlation value of a single-photon state. Thus, we want to minimize Expr. (19) , to approach a single-photon state as closely as possible. Any value 0 ≤ gl (V ) < 1 is considered a genuine signature of quantum statistics, referred to as antibunched statistics (∼ the particle nature of light)Density matrix entries. A considerable problem with the optimization of gl (V ) = 0 is that optimal values are typically obtained at very low photon density (intensity) n l -there is no a priori constraint for ensuring high n l .Therefore, besides optimizing the accuracy of a single-photon outcome, we also aim at maximizing the probability for this to occur. We propose an objective expressed in terms of density-matrix entries of the density matrix ρ of output mode l,EQUATIONThat is, we want to maximize the single-photon matrix entry ρ11 and, with a lower weight (0 < w 0 < 1), the vacuum entry ρ00 . The occupation of any other matrix entry is penalized (last term). Note that a non-zero value vacuum entry ρ 00 is not penalized in (20) . The relatively low input intensities, α i ≤ 1, necessitate a non-zero occupation of the vacuum state. Moreover, if no single photon is emitted, we prefer to have a clear-cut zero signal, instead of generating spurious signals with higher photon numbers.The matrix entries are expressed using the transformed operators from (17) ,EQUATIONwith ρ out the density matrix representing the output state of the first, nonlinear optical network. C is some constant of normalization -diagonal entries of ρ are probabilities, so trρ = 1 by definition.Both objectives are conveniently computed as tensor contractions. For this, the stored MC quantum states |ψ s ⟩ are loaded and appropriately contracted with the PIC unitary V , each iteration in the optimization. As an example, obtaining the fourth-order correlation from Expr. V * l,i1 V * l,i2 V l,i3 V l,i4 a † i1 a † i2 a i3 a i4 s .The fourth-order correlator can be precomputed and stored as a rank-4 tensor, to be reused for every optimization step,EQUATIONStarting from the generated set of quantum Monte Carlo trajectory states |ψ k ⟩, the PIC unitary V was optimized using PyTorch [45] with the LBFGS solver.Quantum Tomography. The single/two-mode squeezed thermal state is injected into the third (and fourth) waveguide out of the five waveguides of the PIC. Into the remaining four/three waveguides, we inject coherent states with displacement amplitudes and phases, in order, 0.5, 0.1, (0.1,) 0.5 and -π/2, -π/4, (π/4,) π/2 respectively. We consider five layers of interference regions in a brickpattern, each of duration ∆t = 0.1/γ, and initialize the hopping coefficients randomly in the range For the optimization task for single-photon generation in the second linear PIC, we used the LBFGS solver of the PyTorch library [45] . The graph matrix G was represented in the basis of trace zero Hermitian matrices, so that V is a unitary from the special unitary group SU(L), having the condition detV = 1. The coefficients of the basis matrices, N 2 -1 real numbers, were chosen randomly from the normal distribution (so µ = 0 and σ 2 = 1) for the initial value. Each datapoint (α, U ) started from the same presampled random condition to guarantee consistency for comparison. 



Reviews: review: title: title: Timely study on the role of polaritons in quantum architectures optimized using ML
main_review: main_review: I have read the paper "Designing Polaritonic Integrated Circuits for Quantum Processing" with interest. This is a well written study with a  clear message on the potential of polaritonic circuitry for quantum information processing enhanced by ML techniques (extreme learning techniques). The timing is quite appropriate since a recent report from the Rapaport group has shown impressive progress in using dipolaritons in waveguides for enhanced nonlinearities, https://journals.aps.org/prx/abstract/10.1103/PhysRevX.14.031022 , also there are several EU funded projects aimed towards bringing polaritons closer towards neuromorphic platforms [e.g. EIC Pathfinder Open – 2023 „QUONDENSATE” ; EIC Pathfinder Open – 2023 „PolArt” ; EIC Pathfinder Challenges – 2022 „Q-One”]. Also, https://doi.org/10.1021/acs.nanolett.0c04696 is a relevant reference (albeit classical polaritonics) which the authors might benefit from.

From a modeling perspective, the equations are quite well known to a quantum optics system with a local chi-3 nonlinear term and losses (Lindblad master equation approach). None-the-less, I think the study is nice because the authors are not pumping the system aggressively and creating a lot of dephasing exciton background, which permits them to work with standard quantum optical Hamiltonian with interactions. But I am curious why the authors don't mention the electric field (gates) which is important to get dipolaritons (see PRX). Would electrical contacts be a problem from a design perspective in the context of directional couplers and phase shifters?

The authors study has the word "polariton" in the title, but I struggled to understand ... why. It seems that the only "polariton" aspect of the study is the choice of parameter U/gamma = 0.1 based on recent experiments. I wonder if this justifies calling the study "polaritonic" without explaining in more details what are the materials needed to achieve strong light-matter coupling, is it a cryogenic setup to stabilize the excitons or is it room-temperature, etc. Otherwise it feels like the study is about arbitrary lossy quantum optical system with waveguides. This is not necessarily a bad thing, but it gives of the feeling that polaritons are being used as a "hook".

Another comment is regarding the training of the structure described in lines 123-129. I got a bit lost there. Do I understand correctly that the nonlinear layer is trained for tomographic purposes, but for quantum state generation it is kept fixed and random? And of course, the linear layers are trained in each case? This could maybe be made a bit more explicit/clear in the text since the schematic figure says that the nonlinear layer is fixed but then the text seems to say otherwise.

The results of the study look excellent. I have no criticism here. It seems that the system performs well and the underlying quantum-optical equations and ML methodology applied is well founded.

Lines 186-187 are rather strange with a statement that the g_2 is always positive and then a following statement syaing that it can be zero (contradictory). It would be better to say that g_2 is strictly "non-negative".