Scalable unitary computing using time-parallelized photonic lattices

2.1 Synchronized PPTCs

Implementing high-dimensional and universal unitary operations using light signals [13], [14], [23], [24], [25], [26], [27], [28] has become an urgent issue in achieving photonic deep learning acceleration [29], [30] and quantum gate operations [7], [31], [32]. Although recent studies have enabled the use of higher-dimensional units for more compact design [23], [24], [25], [26], [27], [28], the most conventional approach for analytically deterministic design procedures relies on the factorization of a matrix U _N ∈ U(N) into a set of SU(2) matrices and a diagonal matrix [13], [14]. Each step of factorization is designed to achieve the nulling of an off-diagonal component, by applying a universal SU(2) operation to the unitary matrix U _N . Such a nulling process continues sequentially until the matrix becomes diagonal, thereby allowing for the deterministic configuration of U _N through the cascaded multiplications of physically achievable SU(2) and diagonal matrices.

The sequence of the nulling processes can be implemented in different forms, for example, by nulling the components of each row (or column) sequentially – the Reck design [13] – or by nulling from higher- to lower-order off-diagonals – the Clements design [14]. These sequences lead to distinct circuit depth, as shown in the Clements design with a half of the circuit depth compared to the Reck design. However, both approaches eventually possess O(N) scaling in the circuit depth, leading to O(N ²) scaling in the overall device footprint.

To address the scalability issue, recent studies have proposed the utilization of wave evolutions along alternative physical axes, such as frequency synthetic dimension [33], [34], [35], [36] or temporal perturbations [15], [16], [37], [38], [39], [40]. In these approaches, one popular candidate is to utilize coupled-resonator optical waveguides (CROWs) [41], [42], [43], which enable dynamical modulation of the system for temporal- [15], [16] or frequency-domain [33], [35] tailoring of optical states. In the time-domain computing approach based on the concept of PPTCs [15], [16], the unit gate consists of a pair of resonators coupled through the zero-field loop couplers [42], [43] (Figure 1a), which can be modulated with refractive index changes. While each resonator supports the clockwise and counter-clockwise resonance modes, we focus on the clockwise mode, which leads to the Hamiltonian H governing the gate dynamics [15]:

where ω ₀ is the reference resonance frequency, Δω(t) is the resonance perturbation, ξ ^L (or ξ ^U) characterizes the phase shifts of the lower (or upper) loop coupler, τ is the lifetime determined by the coupling between a resonator and a loop coupler, and σ ₀ and σ _x,y,z are the identity matrix and Pauli matrices, respectively (Figure 1a).

Figure 1:

Synchronized operation of PPTC. (a) SU(2) gate composed of two resonators coupled via lower and upper loop couplers. Each loop coupler includes two phase shifters (colored boxes): ±ξ ^L and ±ξ ^U for lower and upper loop couplers, where +ξ ^L and +ξ ^U (or –ξ ^L and –ξ ^U) denote the lower (or upper) phase shifters of each coupler, respectively. Each resonator includes the phase shifter for resonance perturbation Δω(t) (gray boxes). (b,c) Universal unitary operations U _N using the PPTC for N = 4: serial execution protocol (b) and parallel execution protocol using buffering (c). D ∈ U(4) denotes the designed 4-dimensional diagonal unitary matrix. (d) Buffering function achieved with the SU(2)-gate structure. In (a, d), dashed boxes and solid arrows denote the coupling regions and the direction of wave propagation, respectively. Dashed arrows represent the propagation length of an integer multiple of the effective wavelength of light.

When we set Δω(t) = 0 and ξ ^L = ξ ^U = ξ (Figure 1a), the system Hamiltonian becomes H = −ω ₀ σ ₀ − (1/τ)[σ _x cosξ + σ _y sinξ] according to Eq. (1). Therefore, the gate yields an arbitrary and reconfigurable rotation operation on the Bloch sphere about the axis n = (cosξ, sinξ, 0), where the amount of the rotation is determined by the time Δt maintaining the system parameters and the corresponding H [15]. For the modulation range of 0 ≤ ξ < 2π, the SU(2) gate operations required for decomposing arbitrary U(N) matrices are guaranteed by controlling Δt.

These gate operations serve as the building blocks for realizing universal U(N) operations. To implement a given U(N) operation, we can construct a one-dimensional (1D) resonator lattice composed of the resonators coupled via the coupler in Figure 1a, while mapping each resonator mode to a basis. For an initial resonance state [ψ ₁ ^I, ψ ₂ ^I, … ψ _N ^I]^T, the cascaded, designed SU(2) gate operations for pairs of resonators (blue boxes in Figure 1b) implement the target matrix U _N ∈ U(N) with the final diagonal matrix (the red box in Figure 1b) that can be achieved with elementwise resonance detuning [15], [16]. We note that the SU(2) operations are executed sequentially in time [15], [16], in contrast to the spatial arrangement in the spatial PPCs [3], [4], [5], [6], [7], [8], [9], [10], [11], [12].

However, the critical hurdle of the PPTC described in [15], [16] is the poor scaling along the temporal axis, which exhibits O(N ²) scaling. This limitation originates from the operation principle of an SU(2) gate: the state rotation on the Bloch sphere of which the speed is determined by the coupling strength between resonators and loop couplers. Due to the static coupling determined by the stationary location of optical elements (dashed boxes in Figure 1a), the necessary times to achieve the target SU(2) operations are differentiated for each nulling process and the form of U _N . Consequently, the temporal depths of SU(2) gates are not to be the same, leading to asynchronization and inevitably requiring sequential rather than parallel computing process (Figure 1b). The resulting O(N ²) temporal circuit depth not only degrades the operation speed but also enforces the necessity of ultrahigh Q-factor resonators to store light during the computation [15].

To resolve this poor time scaling issue, we develop the synchronization scheme for the PPTC. Figure 1c shows the proposed design protocol to resolve the time-scaling issue, by synchronizing multiple independent gate operations to the longest-latency one using the buffering functions (white boxes with 2 × 2 identity operations). This can be implemented by the odd-parity configuration with ξ ^U(t) = −ξ ^L(t) = π/2, which leads to decoupling between resonators through the system Hamiltonian H = −ω ₀ σ ₀ [15].

2.2 U(N) operation with buffering

To demonstrate the proposed scheme, we design the phase shifters of the PPTC under the protocol described in Figure 1c. As an example, we investigate the state evolution for a U ₄ matrix. The decomposition based on the Clements design for U _N [14] requires the circuit composed of (N + 1) layers: N layers for SU(2) operations and 1 diagonal-operation layer (Figure 2a). Each layer includes n (≈N/2) SU(2) gates in parallel, executing SU(2) matrices {T ₁, T ₂, … T _n } during distinct operation times {Δt ₁, Δt ₂ … Δt _n }. Because Δt _i ≠ Δt _j in general, we introduce an additional buffer duration of max(Δt ₁, Δt ₂ … Δt _n ) − Δt _i for each operation T _i after completing the designed gate operation. This ensures that all gates in the set complete the operation synchronously with a common cycle period of max(Δt ₁, Δt ₂ … Δt _n ).

Figure 2:

Synchronized U ₄ PPTC example. (a) Schematic of the PPTC and its time coded modulations: the resonance shift for the n-th resonator, Δω _n , and the phase shifts ξ _n,(n+1) ^L,U for the loop coupler between the n-th and (n+1)-th resonators [15]. Each modulation is scheduled to satisfy synchronization. Considering the device bandwidth in real implementation, an ideal set of modulation signals are then low-pass filtered with the cut-off frequency ω _c = 0.01ω ₀. (b, c) Bloch-sphere representations of the pseudospinor states ψ _1,2 = [ψ ₁, ψ ₂]^T (b) and ψ _3,4 = [ψ ₃, ψ ₄]^T (c), which are defined by the resonance modes of the neighboring resonators 1 and 2, and 3 and 4, respectively. The Stokes parameters S _i∈x,y,z = ψ _n,(n+1) ^T σ _i ψ _n,(n+1) are plotted on the unit sphere. The evolutions of the states are plotted during time t ₀ ≤ t ≤ t ₁, where t _b denotes the completion time of the SU(2) operation between the resonators 3 and 4. Black stars denote the target states. The buffering function is observed in-between t _b (marked in purple) and t ₁ (marked in red).

Figure 2a illustrates the designed temporal modulation of the phase shifters under this protocol. In cycle 1, the gate operations T ₁ and T ₂ are executed in parallel, leading to the required SU(2) operations for the resonator pairs (1, 2) and (3, 4), respectively. Upon completion of T ₂ at time t _b, the (3, 4) resonator pair is on the buffering mode, freezing the pseudospinor state without any energy exchange or the emergence of relative phase shift.

From the designed modulation, we calculate the PPTC operation following the same numerical method in [15], [16]: the application of the sixth-order Runge–Kutta (RK6) method [44] to solve the time-dependent equation idΨ/dt = HΨ with Eq. (1). The operations based on buffering are visualized on the Bloch sphere in Figure 2b and c. Because T ₁ has the longest latency execution time in this cycle, the designed SU(2) operation on [ψ ₁(t ₀), ψ ₂(t ₀)]^T completes at time t ₁ (Figure 2b). Meanwhile, the operation T ₂ completes earlier at time t _b, and the state of [ψ ₃(t _b), ψ ₄(t _b)]^T remains stationary on the buffer mode until t ₁ (Figure 2c). As a result, the pseudospinor states of both pairs [ψ ₁(t ₀), ψ ₂(t ₀)]^T and [ψ ₃(t ₀), ψ ₄(t ₀)]^T reach their target output state simultaneously at t ₁, demonstrating the state synchronization via buffering.

3.1 Time scaling

To assess the impact of the proposed synchronization scheme, we evaluate the operation time T required to fully implement an arbitrary unitary matrix U _N . Figure 3a presents the time scaling evaluated using 50 random Haar unitary matrices for each N [45], where the power-law fitting provides the exponent α = 2.10 for the serial scheme and α = 1.01 for the synchronization scheme. We observe an apparent improvement of time scaling from O(N ²) to O(N) under the proposed synchronization scheme, reflecting the reduction in required cycles from M = N(N – 1)/2 + 1 to M = N + 1 for U _N implementation.

Figure 3:

Time and Q-factor scaling. (a) Time scaling for N. The power-law exponents are α = 2.10 and α = 1.01 for the serial and synchronization schemes, respectively. (b) Q-factor scaling for N. In (a, b), blue and pink dots represent the 50 realizations for the serial and synchronization schemes, respectively. The colored dashed lines denote the power-law fittings of the data using the mean square error. The black dashed line in (b) corresponds to Q = 4 × 10⁷ experimentally demonstrated with the TFLN ring resonator [46]. The achievable N estimated with the fitted Q is marked as a black dot each for the serial and synchronization schemes.

Because our time-domain computation is conducted while light is stored in the resonators, reducing the operation time directly alleviates the requirement for high-Q optical elements. In Figure 3b, we estimate a lower bound of the intrinsic Q factor of each resonator via Q = ω ₀ T, where the full computation time T is obtained using the lifetime parameter τ = 10T _M to ensure the conservation of sufficient optical signals during the computation. We calculate the bound for 10 GHz modulation bandwidth of T _M ∼ 0.1 ns timescale, which is an accessible number with thin-film lithium niobate (TFLN) modulators [47]. The calculated Q in Figure 3b demonstrates the significant alleviation for implementing high-dimensional unitary operations under the synchronization scheme, enabling U ₂₁ operations with the Q = 4 × 10⁷ resonators, in contrast to the U ₈ limit in the prior serial scheme [15].

3.2 Fidelity analysis

Given the substantial reduction in operation time achieved with the synchronization scheme, we can reasonably expect improved operational fidelity at the same level of system defects, because shorter operation times reduce the exposure of light to defects. To examine this prediction, we evaluate various types of defects relevant to our architecture: fluctuations in the resonance frequency ω ₀ and the coupling lifetime τ, the altered zero-field condition of loop couplers, and resonator radiation loss. We characterize all these defects using a uniformly random distribution unif(a, b): ω ₀′ = [1 + unif(−σ, σ)]ω ₀, ω ₀″ = [1 + i × unif(0, 2σ)]ω ₀, 1/τ′ = [1 + unif(−σ, σ)]/τ, and Δφ = unif(−σ, σ)π, where ω ₀′ and ω ₀″ denote the real- and complex-valued resonance frequencies accounting for resonance fluctuations and radiation loss, respectively, τ′ denotes the perturbed lifetime defined for the uniformly random perturbation of the coupling between a resonator and a loop coupler, and Δφ is the unwanted phase shift exerted on each arm of the loop couplers.

For the target unitary operation U _N and its PPTC realization V _N , we examine the operational fidelity with the following Frobenius-norm error:

where Tr[A] denotes the trace of matrix A. Although the error ε satisfies 0 ≤ ε ≤ 2 for the unitary output V _N , this bound no longer applies to nonunitary V _N . The defects characterized by the altered zero-field condition of loop couplers and resonator radiation loss lead to nonunitary V _N because optical energy is transferred into the loop coupler and environment, respectively, and therefore, the resonator-system Hamiltonian is no longer Hermitian.

To estimate ε statistically, we follow the procedure in the previous work [15]. In obtaining the V _N for one target U _N , we calculate the PPTC operations for the N input states ψ _n ^I (n = 1, 2, …, N), which are obtained by normalizing N × 1 complex vectors that have the amplitudes unif(0, 1) and the phases unif(0, 2π). For initial field distributions of {ψ _n ^I} inside the PPTC, we obtain N different outputs ψ _n ^O using the RK6 method, which allow for estimating V _N directly through the following inversion, as

In evaluating ε, we examine 10 random Haar unitary matrices as the target unitary operations.

Figure 4a–d shows the σ-dependent errors for the serial and synchronized implementations at N = 5 and N = 20. Although the detailed dependencies vary with the type of defects – for example, more sensitive responses to resonance perturbations than to the lifetime variations – the synchronization scheme consistently exhibits superior robustness against defects compared to the serial implementation. Notably, such an improvement in robustness by the synchronization is evident especially in larger systems (N = 20 in Figure 4(b, d, f, h) versus N = 5 in Figure 4(a, c, e, g)). This trend also holds for other values of N as shown in Supplementary Note S1.

Figure 4:

Robustness against defects. (a–d) Variations of errors with respect to the degree of defects σ: (a, b) resonance fluctuations, (c, d) lifetime fluctuations, (e, f) resonator radiation loss, and (g, h) phase-shift defects in loop couplers. Blue and pink dots denote ε obtained with 10 realizations of random Haar matrices at each value of σ. Solid lines denote their averages. We compare the serial (blue) and synchronization (pink) scheme for small- (N = 5) and large-scale (N = 20) systems for each type of defects.

3.3 Analysis for real implementation

We evaluate the computing performance of our PPTCs in terms of real implementation, including the comparisons with the prior serial configuration [15], [16] and with other electronic [48] and photonic [29] platforms. Our target calculation – U _N matrix multiplication to an N × 1 vector – includes N ² multiplications and N(N − 1) additions of complex scalars. Because a multiplication and an addition of complex scalars require six and two operations regarding real-valued scalar calculations, respectively, the target calculation requires 8N ² – 2N real-valued operations.

The operation times of the PPTCs are governed by the SU(2) gate operation time, which is comparable to τ. Because the time scalabilities of the serial [15], [16] and our synchronized implementations are O(N ²) and O(N), respectively (Figure 3), the operations (OPs) per second (OPS) of PPTC become ≈ (8N ² – 2N)/(N ² τ) ≈ 8/τ OPS in the serial design [15], [16] and ≈ (8N ² – 2N)/(Nτ) ≈ 8N/τ OPS for our design. This estimation demonstrates that the parallelization enables scalable photonic computing with the operation speed proportional to the system size N. To reach the cutting-edge electronic GPU performance close to 10³ tera-OPS (TOPS) [48], the synchronized PPTC requires N/τ ≈ 10¹⁴, which can be achieved with 10⁴ resonators under ultrafast electro-optic modulation (τ ≈ 0.1 ns, [47], [49], [50]). Considering the difficulty of achieving the stable coupling of a huge number of resonators, the application of systolic matrix multiplications [51], [52], [53] can be a candidate using a lattice of about 10² coupled resonators.

We also estimate energy efficiency, by assuming the modulator in [49], which supports 0.63 mW power consumption and 5 mm² footprints. According to its O(N) scaling of the gate number, the PPTCs support ≈ 8,000/(Nτ) OP/J with the serial process [15], [16] and ≈ 8,000/τ OP/J with the synchronization design. Notably, our scheme using ultrafast electro-optic modulation surpasses the energy efficiency of the cutting-edge electronic GPU that supports a few TOP/J.

In terms of substituting the spatial axis with temporal axis, we further examine a comparison with the spatial PPCs [29]. We note that PPCs provide superior operation-speed scalability with ≈ (8N ² – 2N)/τ _D OPS and comparable power efficiency ≈ 8,000/τ _D OP/J, where τ _D is the characteristic time of photodetection and the power efficiency is governed by the O(N ²) phase-shift bias power for large N. When assuming τ _D = 0.1 ns detection time, the PPC achieves 10 TOPS and 100 TOP/J with N ≈ 10, which leads to the footprint of 500 mm² – comparable to the cutting-edge GPU die size of 10³ mm². However, the PPC suffers from poor scalabilities in the footprint and gate number, leading to 5 × 10⁴ mm² footprint and 10⁴ phase shifters for N ≈ 100, hindering the increase of computing dimension despite superior operation speed of 10³ TOPS. Therefore, the major advantage of the PPTC lies in the scalability of the gate number and footprint: 500 mm² footprint and ∼500 phase shifters for N ≈ 100.

We note that the spatio-temporal computing circuit composed of waveguides and resonators, which has been very recently studied for the different purpose [54], [55] – nonlinear expressivity in photonic deep learning – would be a candidate platform resolving both the dimensionality limit of PPCs and the speed limit of PPTCs. This configuration also maintains the potential advantage of PPTCs as all-optical programmable devices – the use of resonance modes for energy localization and the consequent enhanced nonlinearity – which has been partly examined in [54], [55].