Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J.,
et al. 2010, , 404, 1639, 10.1111/j.1365-2966.2010.16407.x

[Viaene et al., 2017]
Viaene, S., Sarzi, M., Baes, M., Fritz, J., & Puerari, I. 2017,
, 472, 1286, 10.1093/mnras/stx1781

[Walsh et al., 2008]
Walsh, J. L., Barth, A. J., Ho, L. C., et al. 2008, , 136, 1677,

[Wright, 2006]
Wright, E. L. 2006, , 118, 1711, 10.1086/510102

[Yıldırım et al., 2017]
Yıldırım, A., van den Bosch, R. C. E., van de Ven, G., et al.
2017, , 468, 4216, 10.1093/mnras/stx732

[Yoon, 2017]
Yoon, I. 2017, , 466, 1987, 10.1093/mnras/stw3171

[Young et al., 2011]
Young, L. M., Bureau, M., Davis, T. A., et al. 2011, , 414, 940,

[Zabel et al., 2019]
Zabel, N., Davis, T. A., Smith, M. W. L., et al. 2019, , 483,
2251, 10.1093/mnras/sty3234

[Zabludoff & Mulchaey, 2000]
Zabludoff, A. I., & Mulchaey, J. S. 2000, , 539, 136,

[Zhao et al., 2021]
Zhao, Y., Ho, L. C., Shangguan, J., et al. 2021, , 911, 94,

§ SURFACE BRIGHTNESS PROFILES AND ISOPHOTAL ANALYSIS RESULTS

Here, we provide more detailed and individual analysis of the $H$-band surface brightnesses. In Figures <ref>-<ref>, we provide isophotal analysis results for each ETG as a function of radius nearly out to the edge of the WFC3/IR mosaic FOV. These figures also include 1D surface brightnesses extracted along the major axis from both the $H$-band data and the best-fitting 2D MGE decomposition using .

These surface brightness measurements were made after masking out the most dust-obscured regions with $\Delta (J-H) \gtrsim 0.08$ mag, which leads to (sometimes large) gaps in major-axis plots but only a reduced number of data points for the 2D fits. The majority of our MGE fits show good agreement to the $H$-band data, with fractional differences generally at the 10% level or less across all (relevant) radii. The major axes of the dust disks tend to be fairly well aligned with the respective stellar photometric axes, although some targets do show significant ($>$20) divergence between PA$_\star$ and PA$_{\mathrm{dust},H}$. The isophotal analyses for these targets are unreliable within the CND radii but reveal smoother variations in the fitted ellipse PA and ellipticity beyond the respective dust extents. The $a_4/a$ deviations are likewise small, with some tending toward boxy or disky isophotes at large radii. In nearly every case, the best-fit MGE reproduces the central stellar light distribution quite well, and should therefore be useful in future/ongoing ALMA CO dynamical modeling efforts. For just under a third of the sample, however, constraints on the MGE components worsen the global goodness-of-fit. Below, we briefly discuss particulars for all fits individually.

* Hydra A (Figure <ref>): Despite showing close agreement between the stellar and disk orientation at $R\sim R_{\mathrm{dust},H}$, the isophotal analysis shows a rapid shift in PA$_\star$ and a steady increase in $\varepsilon_\star$ starting just afterwards. The $a_4/a$ deviations are generally small but show some preference for boxy isophotes. Such gradients and boxiness are consistent with the dense galactic environment and merger history for BCGs. The best-fit MGE has PAs that align with $\overline{\mathrm{PA}}_\star$, which here hinges on the greater number of measurements at larger radii. Nevertheless, the generally good agreement (with data$-$MGE fractional differences of $\lesssim$5%) is due to fairly circular isophotes at small $R$.

* NGC 612 (Figure <ref>): The MGE shows passable agreement to the $H$-band data out to $R\sim 70\arcsec$ despite the challenges inherent in fitting such a large dust disk. Of course, the high $\Delta \mathrm{PA}_\star > 45\degr$ cannot be reproduced by the MGE, but the low $\varepsilon_\star$ over most of this $\Delta \mathrm{PA}_\star$ range allows for MGE components with a uniform PA to still fit the data without large discrepancies. The $a_4/a$ values show clear preference for disky isophotes between 40$-$50 with no clear trend beyond, although higher $\varepsilon_\star$ towards the mosaic edge suggests overall greater flattening. Interior to these radii, the outer dust disk features are generally aligned to within $\sim$20 of PA$_\star$.

* NGC 997 (Figure <ref>): This E galaxy shows close agreement between dust and stellar orientation until $R \sim 20\arcsec$. This radius coincides with the angular distance of a possible companion, PGC 200205. NGC 997 is tidally interacting with NGC 998, and these two other galaxies may both contribute to the apparent shift in PA$_\star$. The MGE shows excellent agreement in the central couple $\times$10 pc but the fit remains good out to $R\sim70\arcsec$. While $\varepsilon_\star$ remains flat, the limiting $q^{\prime}$ still limits most of the MGE components, although the global fit only appears slightly affected. In addition, the $a_4/a$ values shows a clear preference for boxy isophotes with increasing radius. In either case, the global fit appears to be only slightly affected. We note that Dominiak et al., 2024 ran an independent MGE model in the H-band for this galaxy, finding an additional compact component (with $\sigma^{\prime} < 0\farcs06$) at lower luminosity that does not appear necessary based on the MGE residuals in Figure <ref>. However, these differences can likely be explained by the choice of PSF and lack of masking.

* NGC 1332 (Figure <ref>): The dust disk orientation remains closely aligned to PA$_\star$ at all radii, even while the isophotes become increasingly flattened. Despite having disk-dominated outskirts, this lenticular galaxy shows slightly boxy isophotes. The MGE reproduces the observed $H$-band stellar light distribution well, although a slight PA$_\star$ shift and flattened isophotes beyond $R \sim 60\arcsec$ accentuate this discrepancy along the major axis.

* NGC 1387 (Figure <ref>): The measured PA$_\star$ reveals one of the most misaligned ($>$50) dust disks in our sample, with the stellar isophotes lining up well with the stellar bar instead of PA$_{\mathrm{dust},H}$. Given the nearly face-on CND orientation, this 2D MGE cannot match the flattened isophotes around the bar-dominated regions. Because of the CND shape and uniformly low $\Delta (J-H)$ observed for this ETG, we do not create any dust mask. The nearly face-on CND should not affect either the isophotal results or the MGE $\sigma^{\prime}$ or $q^{\prime}$ values; however, MGE central surface brightness $I_H$ may be slightly suppressed, especially for the innermost components that lie almost entirely within the large-scale $R_{\mathrm{dust},H}$.

* NGC 3245 (Figure <ref>): The dust disk is generally well-aligned with PA$_\star$ except in the first couple arcseconds after $R_{\mathrm{dust},H}$. However, the more circular isophotes at these radii limit the data$-$MGE discrepancies. At larger radii, the more flattened isophotes do result in a couple of MGE component $q^{\prime}$ values settling to $q^{\prime}_\mathrm{min}$, but the overall fit does not seem to be affected.

* NGC 3258 (Figure <ref>): The best-fit MGE is good at all radii, as is demonstrated along the major axis, although the data do show a moderate PA$_\star$ shift from PA$_{\mathrm{dust},H}$ past $R \sim 20\arcsec$. The data–MGE discrepancy that would otherwise arise from this $\Delta \mathrm{PA}_\star$ is lessened for surface brightness measurements between $\sim$20$-$50 because of more circular isophotes at these radii.

* NGC 3268 (Figure <ref>): This galaxy shows nearly perfect power-law surface brightness behavior beyond $R_{\mathrm{dust},H}$ with only minor changes in PA$_\star$ and $\varepsilon_\star$ and no significant $a_4/a$ deviations from pure elliptical isophotes. As a result, the MGE provides a very good global fit.

* NGC 3271 (Figure <ref>): The consistent and moderate-amplitude $\mathrm{PA}_\star - \mathrm{PA}_{\mathrm{dust},H}$ offset arises in part due to the MGE fitting to a prominent, misaligned stellar bar, although the PA$_\star$ of the stellar disk remains persistently high out to $R \sim 50\arcsec$. The generally poorer MGE for this cuspy galaxy also reflects the $a_4/a$ values changing from very disky behavior around the bar to boxy features with increasing radius.

* NGC 3557 (Figure <ref>): Good alignment between PA$_{\mathrm{dust},H}$ and PA$_\star$ at all radii helps this cored radio galaxy be fit well with an MGE. While we include a central PSF component in the fit, at these wavelengths the AGN contributions are essentially negligible. Beyond $R \sim 50\arcsec$, the $\varepsilon_\star$ decrease is coincident with a gradual increase in $a_4/a$ towards disky behavior.

* NGC 3862 (Figure <ref>): This galaxy shows a nearly uniform PA$_\star$ gradient while becoming increasingly disk-dominated beyond . Despite the large change in PA$_\star$, the MGE does not result in a poorer fit until $R \gtrsim 20\arcsec$, when the isophotes become less circular. The measured $a_4/a$ values show some preference for boxy shapes only near the edge of the WFC3/IR FOV. The central surface brightnesses for this radio galaxy are consistent with a point source, and including a PSF component recovers the expected core behavior.

* NGC 4061 (Figure <ref>): Out to $\sim$2, the stellar isophotal behavior remains consistent and the measured PA$_\star$ aligns with the observed dust disk orientation. At larger radii, an abrupt shift in PA$_\star$ is accompanied by boxy isophotes. The more disturbed features appear to be due to interactions with the nearby NGC 4065. However, the tidal influence is limited to the outer regions of the WFC3/IR FOV, and the  MGE reproduces the inner 2D stellar surface brightnesses well.

* NGC 4261 (Figure <ref>): The relatively close agreement between stellar isophotal and outer dust orientations, together with overall low $\varepsilon_\star$, allows for a good global fit without any preference for disky or boxy isophotes. Interestingly, the CND and stellar PAs remain offset by $\sim$5 out to . Boizelle et al., 2021 show that the outer CO kinematics appear to match PA$_{\mathrm{dust},H}$ reported in Table <ref>. However, the CO kinematic line-of-nodes PA approaches $\sim$0 towards the disk center, which only increases the $\mathrm{PA}_\star - \mathrm{PA}_{\mathrm{dust},H}$ offset. We note that the MGE includes a central PSF component fitted simultaneously in , but this point source is negligible at these wavelengths.

* NGC 4373a (Figure <ref>): The MGE fit shows generally good agreement out to $R\sim13\arcsec$. At larger radii, high $\varepsilon_\star > 1-q_\mathrm{min}^{\prime}$ for $R\gtrsim 15\arcsec$ and isophotes with $a_4/a > 0.01$ for $15\arcsec \lesssim R \lesssim 40\arcsec$ lead to a poorer overall fit. Indeed, most of the MGE components have $q^{\prime}$ that are constrained by $q_\mathrm{min}^{\prime}$. Inside the disky isophotes, the stellar isophotal analysis returns $a_4/a \approx -0.01$ and the inner MGE fit around $R\sim 10\arcsec$ struggles to match the more boxy features. The dust disk is aligned to within $\sim$10 of the stellar photometric axis at most radii; at large $R$, PA$_\star$ may begin to diverge significantly from PA$_{\mathrm{dust},H}$, although the certainty of these final measurements is not high.

* NGC 4429 (Figure <ref>): Overall, the MGE shows generally good agreement out to $R\sim45\arcsec$. $\mathrm{PA}_\star \sim \mathrm{PA}_{\mathrm{dust},H}$ just beyond $R_{\mathrm{dust},H} \sim 13\arcsec$, although PA$_\star$ increases thereafter to a maximum $10\degr$ offset. The isophotal analysis shows an $I(R)$ bump for $50\arcsec \lesssim R \lesssim 90\arcsec$ that corresponds to increased PA$_\star$, $\varepsilon_\star$, and $a_4/a$ values. Figure <ref> reveals this feature arises from the (inner) stellar ring. All of the outer MGE components have $q^{\prime} \rightarrow q_\mathrm{min}^{\prime}$, and the global MGE does not reproduce the outer stellar light distribution as well as it does the inner regions. This  fit includes a PSF component to account for moderate AGN contamination, although the central MGE component(s) exceed(s) the central surface brightness  of this unresolved source.

* NGC 4435 (Figure <ref>): This barred lenticular galaxy shows only modest PA$_\star$ changes and moderate $\varepsilon_\star$ fluctuations that are below $1-q_\mathrm{min}^{\prime}$ for $R < 70\arcsec$, allowing the MGE to reproduce the inner stellar light distribution well. In this region, fluctuating $a_4/a$ indicates both disky and boxy isophotal behavior at different radii. In the outer regions (at $R \gtrsim 80\arcsec$), the PA$_\star$ abruptly shifts while the isophotes become very disky. Some portion of this shift may be due to the stellar light distribution of (or perhaps interactions with) neighboring NGC 4438.

* NGC 4697 (Figure <ref>): The MGE provides very good fits to both inner and global stellar light distributions, with $\mathrm{PA}_\star \sim \mathrm{PA}_{\mathrm{dust},H}$ at all $R$. The measured $a_4/a$ values approach boxy behavior towards the outer parts of this galaxy without a noticeable impact on the MGE goodness-of-fit. Positive data$-$MGE residuals along the major axis for $R\gtrsim 15\arcsec$ are evidence of a disk-like stellar component ($\varepsilon_{\mathrm{disk},\star} > 0.7$) for this fast-rotating ETG [Spiniello et al., 2015] that is not fully fit with this MGE approach.

* NGC 4751 (Figure <ref>): This cuspy lenticular galaxy is moderately well-fit by an MGE at all radii, with $\mathrm{PA}_\star \approx \mathrm{PA}_{\mathrm{dust},H}$ across most of the WFC3/IR FOV. At larger radii, the $a_4/a$ values hint at slightly more boxy features, which result in negative data$-$MGE fractional residuals along the major axis.

* NGC 4786 (Figure <ref>): For this more cored galaxy, the $\mathrm{PA}_\star \sim \mathrm{PA}_{\mathrm{dust},H}$ and $\varepsilon_\star$ values are consistent across nearly the entire FOV, only showing discrepancies for $R\gtrsim 90\arcsec$ coincident with the start of disky isophotes. The MGE reproduces the stellar light distributions well on both inner and global scales.

* NGC 4797 (Figure <ref>): This fairly compact lenticular galaxy shows consistent PA$_\star$ out to $>$4 even as the stellar ellipticity fluctuates. Beyond that radius, PA$_\star$ and $a_4/a$ both increase, with the isophotes transitioning from boxy to disky behavior. The data–MGE residuals are around the 10% level for many of the major-axis measurements, although the 2D fits show generally good agreement throughout the WFC3/IR FOV.

* NGC 5084 (Figure <ref>): The MGE shows good agreement to the data only out to $\sim$7, beyond which the highly elliptical stellar light distribution and constrained $q^{\prime}$ lead to a poor fit, especially along the major-axis direction. Beyond the dust extent, the isophotal analysis shows fairly constant PA$_\star$. However, PA$_{\mathrm{dust},H}$ is nearly orthogonal to $\overline{\mathrm{PA}}_\star$, indicating a polar-ring CND. In such a case, the MGE may still be useful but one cannot simply assume the CND lies in the galaxy's midplane when calculating $v_{\mathrm{c},\star}$. At the transition between disk-dominated regions (out to $R \sim 70\arcsec$) and the stellar halo, the stellar isophotes show a clear trend towards boxy behavior.

* NGC 5193 (Figure <ref>): This peculiar E galaxy is well-fit by the standard MGE. The stellar and dust-disk orientation angles agree just beyond $R_{\mathrm{dust},H}$, although an increasing PA$_\star$ leads to a $\sim$15 offset between $R\sim 20-30\arcsec$. This discrepancy may be due to a companion galaxy to the southwest whose light may not be fully masked. This PA$_\star$ shift does not result in a noticeably worse MGE goodness-of-fit, due to a coincident drop in $\varepsilon_\star$ to nearly circular isophotes.

* NGC 5208 (Figure <ref>): For this cuspy galaxy with a large (projected) disk size, the remaining unmasked $H$-band data show only mild gradients in PA$_\star$ and $\varepsilon_\star$. Despite slightly larger data–MGE residuals along the major axis, the overall MGE is not a poor representation of the $H$-band mosaic.

* NGC 5838 (Figure <ref>): The MGE shows generally good agreement with the data out to $R \sim 35\arcsec$, beyond which highly elliptical isophotes and constrained $q^{\prime} \rightarrow q_\mathrm{min}^{\prime}$ lead to a poorer fit, especially along the major axis. Modest changes in PA$_\star$ for $1< R/R_{\mathrm{dust},H} \lesssim 3$ exacerbate the discrepancies. However, $\mathrm{PA} \sim \mathrm{PA}_{\mathrm{dust},H}$ both near the CND edge and for $R\gtrsim 5 R_{\mathrm{dust},H}$. At more intermediate radii, the isophotes show boxy behavior.

* NGC 6861 (Figure <ref>): The MGE for this BGG is able to reproduce the observed surface brightness distribution fairly well. Despite close agreement between PA$_\star$ and PA$_{\mathrm{dust},H}$ in the inner regions, PA$_\star$ begins to noticeably decrease for $R\gtrsim 50\arcsec$, leading to large-scale data$-$MGE discrepancies. The disky isophotes at large $R$ complicate the global fit, but not more than the observed $\Delta\mathrm{PA}_\star$.

* NGC 6958 (Figure <ref>): In this cD galaxy, the very face-on [$(b/a)_{\mathrm{dust},H}\sim 1$] disk orientation severely limits the $q^{\prime}$ values that the MGE components can assume at all radii, leading to high major-axis data$-$MGE residuals. As a result of a large PA$_\star$ decrease and a moderate $\varepsilon_\star$ increase for $R\gtrsim 40\arcsec$, the MGE solution beyond this radius becomes less reliable.

Isophotal analysis (left panels) and MGE (shown along the major axis; right panels) for Hydra A after masking the most dust-obscured regions of the CND. Each panel gives the dust disk extent (vertical dotted line), beyond which the fitted PA$_\star$ and ellipticity $\varepsilon_\star$ from the  task are reliable and show generally smooth variations. The CND PA$_{\mathrm{dust},H}$ is also shown for reference (horizontal dotted line). The $a_4 / a$ deviations from a perfect ellipse are likewise generally small. The full MGE shows generally good agreement (lower right panel) to the major-axis surface brightness measurements.

Same as for Figure <ref> but for NGC 612.

Same as for Figure <ref> but for NGC 997.

Same as for Figure <ref> but for NGC 1332.

Same as for Figure <ref> but for NGC 1387.

Same as for Figure <ref> but for NGC 3245.

Same as for Figure <ref> but for NGC 3258.

Same as for Figure <ref> but for NGC 3268.

Same as for Figure <ref> but for NGC 3271.

Same as for Figure <ref> but for NGC 3557.

Same as for Figure <ref> but for NGC 3862.

Same as for Figure <ref> but for NGC 4061.

Same as for Figure <ref> but for NGC 4261.

Same as for Figure <ref> but for NGC 4373a.

Same as for Figure <ref> but for NGC 4429.

Same as for Figure <ref> but for NGC 4435.

Same as for Figure <ref> but for NGC 4697.

Same as for Figure <ref> but for NGC 4751.

Same as for Figure <ref> but for NGC 4786.

Same as for Figure <ref> but for NGC 4797.

Same as for Figure <ref> but for NGC 5084.

Same as for Figure <ref> but for NGC 5193.

Same as for Figure <ref> but for NGC 5208.

Same as for Figure <ref> but for NGC 5838.

Same as for Figure <ref> but for NGC 6861.

Same as for Figure <ref> but for NGC 6958.

§ MULTI-GAUSSIAN EXPANSIONS WITH A FREE PA

We report  MGE solutions that allow for PA to vary between components in Table <ref> and show the resulting contour plots in Figure <ref>. Photometric PA twists are common in this sample, and accounting for these PA twists results in better overall fits. For uniformity with our previous MGE solutions, we restrict these models to an appropriate range of inclination angles, which can still result in some discrepancies. Accounting for PA twists can also be useful for stellar-dynamical efforts that explore triaxality <cit.>. However, allowing for a variable PA between Gaussian components prevents a simple deprojection. Isophotal centroids are not always consistent for different radii <cit.>, likely due to ongoing settling of a recent merger or tidal interaction. However, non-concentric series expansions are not currently viable for standard gas-dynamical modeling efforts, so we opt to force a concentric MGE solution. Compared to the MGE solutions with a fixed PA, changes in  are minimal in most cases, reaching up to $\sim4\%$ improvement at best.

MGE Parameters (with a Free PA)

j     $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA

$(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)

(1)     (2)     (3)     (4)     (5)      (2)     (3)     (4)     (5)     (2)     (3)     (4)     (5)

$q^{\prime}_{min} = 0.23$                 $q^{\prime}_{min} = 0.24$                 $q^{\prime}_{min} = 0.91$

3cHydra A         3cNGC 612         3cNGC 997

2-5 6-9 10-13
1     3.4568     1.2335     0.8913     $-$79.83     4.3281     0.1903     0.7587     $-$39.44     4.5662     0.2499     0.9100     36.57

2     2.9176     2.9311     0.9040     $-$30.96     3.8095     0.9655     0.2893     $-$23.97     4.2771     0.7129     0.9127     34.68

3     2.1985     6.6571     0.9179     $-$6.58     3.6600     1.0974     0.2400     $-$40.32     3.9215     1.4906     0.9309     35.48

4     2.0283     7.7789     0.6524     $-$44.11     3.7754     0.8971     0.5199     $-$87.74     3.4213     3.2687     0.9100     30.44

5     2.1393     13.591     0.7972     $-$34.43     3.1582     2.3074     0.8134     $-$84.61     2.9231     7.1770     0.9100     29.83

6     1.3498     27.250     0.7011     $-$29.73     3.8082     3.8546     0.2658     $-$4.70     2.3616     15.141     0.9100     21.84

7     0.9367     59.514     0.6606     $-$33.63     3.0907     4.7088     0.6825     $-$9.43     1.7550     36.268     0.9100     10.67

8                     2.5885     5.9420     0.9211     $-$86.23

9                     2.7923     7.1097     0.2400     $-$8.91

10                     2.3449     12.666     0.4365     $-$9.20

11                     2.2925     16.064     0.8224     $-$22.72

12                     1.1335     52.983     0.6825     $-$31.82

13                     1.1010     53.895     0.5797     $-$84.96

$q^{\prime}_{min} = 0.17$                 $q^{\prime}_{min} = 0.95$                 $q^{\prime}_{min} = 0.52$

3cNGC 1332         3cNGC 1387         3cNGC 3245

2-5 6-9 10-13
1     5.8054     0.1279     0.3348     $-$63.74     5.3558     0.1549     0.9500     $-$7.45     5.6183     0.1000     0.7331     $-$39.10

2     4.9808     0.4945     0.9725     $-$49.77     4.8527     0.5285     0.9744     $-$12.38     5.1971     0.1626     0.5200     $-$13.55

3     4.6118     1.5343     0.7284     $-$62.34     4.5401     1.6483     0.9922     $-$42.47     4.9520     0.4882     0.8925     $-$3.04

4     4.2028     3.4563     0.7279     $-$64.51     4.0631     4.0942     0.9500     $-$73.88     4.3391     1.0742     0.8129     $-$36.30

5     3.8253     7.4695     0.7682     $-$60.44     3.5116     6.8412     0.9500     $-$71.07     4.2465     1.1928     0.9226     $-$19.57

6     3.1793     18.073     0.3133     $-$64.62     2.9294     14.219     0.9500     $-$70.62     4.0351     2.5475     0.6392     $-$3.27

7     3.0388     34.972     0.2839     $-$66.38     2.4982     39.993     0.9599     $-$73.20     3.8187     3.2562     0.7724     $-$7.94

8     2.1287     55.919     0.2547     $-$52.75                     3.5250     8.6738     0.5200     $-$3.53

9     2.3712     60.489     0.3236     $-$71.54                     2.8619     26.303     0.5200     $-$2.97

10     1.5591     86.725     0.6815     $-$55.92                     2.1463     42.513     0.5269     $-$4.76

11                                     2.0451     18.731     0.5200     $-$3.92

12                                     1.3133     109.00     0.9668     $-$15.42

$q^{\prime}_{min} = 0.72$                 $q^{\prime}_{min} = 0.67$                 $q^{\prime}_{min} = 0.73$

3cNGC 3258         3cNGC 3268         3cNGC 3271

2-5 6-9 10-13
1     4.1544     0.7817     0.9534     81.66     3.8018     0.2295     0.6700     $-$64.47     4.9283     0.1000     0.7300     $-$56.05

2     3.9772     1.1312     0.7200     75.57     3.9040     1.0116     0.9786     $-$44.81     5.1004     0.1420     0.7300     $-$71.75

3     3.9493     1.9658     0.7573     77.43     3.9745     1.0911     0.7336     $-$70.80     4.5431     0.3628     0.8602     $-$85.05

4     3.6280     3.0048     0.8116     75.10     3.8481     1.9805     0.7185     $-$66.30     4.2114     0.9708     0.7300     $-$64.29

5     3.4937     4.7532     0.8475     74.75     3.7059     2.2801     0.8904     $-$73.63     4.0583     2.4108     0.7300     $-$60.23

6     2.7343     8.3235     0.8346     84.57     3.4643     3.8046     0.7788     $-$68.70     3.7306     3.7406     0.7300     $-$66.85

7     2.9480     11.477     0.9222     67.16     3.2627     6.3286     0.8091     $-$67.99     3.3797     7.0013     0.7300     $-$69.50

8     2.2619     20.600     0.9724     61.57     2.8096     12.068     0.7836     $-$65.37     2.9196     17.385     0.7300     $-$63.39

9     2.0426     47.745     0.7563     56.47     2.4636     21.539     0.8038     $-$69.14     1.7699     33.410     0.7300     $-$73.50

10                     2.0176     51.615     0.7267     $-$64.49     1.9335     50.589     0.7300     $-$81.49

11                     1.5354     89.674     0.8499     $-$45.60     0.72101     119.88     0.7390     $-$86.53

$q^{\prime}_{min} = 0.75$                 $q^{\prime}_{min} = 0.99$                 $q^{\prime}_{min} = 0.66$

3cNGC 3557         3cNGC 3862         3cNGC 4061

2-5 6-9 10-13
1     4.3101     0.8964     0.8939     40.58     4.1623     0.7526     0.9900     $-$24.23     4.2038     0.1643     0.6600     $-$7.04

2     4.2852     1.3284     0.7500     32.82     3.8002     1.4249     0.9913     $-$16.53     4.2872     0.5382     0.6600     $-$4.16

3     3.9338     2.0829     0.7500     33.04     3.4945     2.9143     0.9900     $-$1.82     3.7492     1.1051     0.9019     $-$5.77

4     4.0034     3.0179     0.7500     34.30     2.7925     7.0073     0.9900     $-$16.24     3.6006     1.2674     0.6600     $-$3.32

5     3.7192     4.5527     0.7606     32.91     2.2356     17.893     0.9900     $-$21.13     3.5601     2.1423     0.8069     $-$6.01

6     3.4630     7.0489     0.7509     34.41     1.6486     63.145     0.9900     $-$47.07     3.1868     3.7987     0.8128     $-$5.24

7     2.6865     17.322     0.7500     32.86                     2.4864     6.3725     0.9544     $-$11.23

8     3.1593     12.247     0.7500     31.70                     2.3310     8.5164     0.6600     $-$9.11

9     2.6519     31.261     0.7500     32.73                     2.1361     13.776     0.8540     $-$6.35

10     2.0201     47.821     0.8375     33.35                     1.6802     24.945     0.6938     $-$2.47

11     1.9106     91.520     0.7986     38.25                     1.4371     50.149     0.6600     $-$26.66

12                                     0.55051     144.60     0.7209     $-$50.73

Individual Gaussian components from the best-fitting MGE for each galaxy in this H-band sample, after masking out neighboring galaxies, foreground stars, and the most dust-obscured regions of the CND. The PA for each component was allowed to vary freely. Projected terms are indicated by a $^\prime{}$. During these fits, the individual $q^\prime$ values were constrained to be equal to or greater than the listed limit $q^{\prime{}}_\mathrm{min}$.

MGE Parameters (with a Free PA)

j     $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA

$(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)

(1)     (2)     (3)     (4)     (5)      (2)     (3)     (4)     (5)     (2)     (3)     (4)     (5)

$q^{\prime}_{min} = 0.71$                 $q^{\prime}_{min} = 0.52$                 $q^{\prime}_{min} = 0.57$

3cNGC 4261         3cNGC 4373a         3cNGC 4429

2-5 6-9 10-13
1     4.3300     1.1171     0.8200     $-$22.92     5.1906     0.1133     0.5200     $-$18.92     5.2744     0.1133     0.8697     $-$89.55

2     4.0998     2.1901     0.7100     $-$20.82     4.8712     0.3296     0.5200     $-$27.01     4.8345     0.3249     0.5700     $-$89.49

3     3.9659     3.7306     0.7257     $-$24.13     4.2773     0.9670     0.5505     $-$21.93     4.3016     0.7510     0.7012     $-$85.52

4     2.9954     6.5190     0.7100     $-$2.520     3.9327     2.1215     0.5788     $-$27.10     4.2403     2.4643     0.5700     $-$89.96

5     3.3170     8.3750     0.7100     $-$29.93     3.4518     4.6950     0.7496     $-$35.00     3.6811     3.1301     0.8392     $-$85.91

6     3.1818     12.620     0.8165     $-$17.42     3.1189     11.913     0.5200     $-$29.96     3.8262     6.6209     0.6155     $-$87.26

7     2.7878     19.976     0.8322     $-$28.76     2.2540     28.705     0.5200     $-$28.46     2.9329     19.703     0.5700     $-$75.08

8     1.7203     41.033     0.7560     $-$7.40     1.7247     54.713     0.5200     $-$29.52     3.2075     15.868     0.5700     $-$88.27

9     2.3330     45.969     0.8291     $-$20.43     0.6228     137.36     0.9964     $-$38.48     2.4565     52.720     0.5700     $-$81.74

10     1.7945     95.693     0.9336     $-$14.80                     2.2329     42.422     0.5700     $-$81.01

11                                     2.4152     71.394     0.5700     $-$83.12

$q^{\prime}_{min} = 0.41$                 $q^{\prime}_{min} = 0.40$                 $q^{\prime}_{min} = 0.35$

3cNGC 4435         3cNGC 4697         3cNGC 4751

2-5 6-9 10-13
1     5.1424     0.1773     0.4820     $-$9.47     5.7171     0.1000     0.4000     $-$32.04     6.1303     0.1136     0.3500     $-$1.66

2     4.6677     0.5383     0.6624     $-$32.01     5.4697     0.1805     0.4117     $-$63.62     5.2241     0.3007     0.3500     $-$4.49

3     4.3822     1.1710     0.8042     $-$4.91     4.8204     0.5119     0.7567     $-$65.81     4.4520     0.5333     0.9973     $-$27.53

4     4.2347     2.4201     0.6675     $-$16.38     4.5250     1.0984     0.7124     $-$66.88     4.4722     0.9159     0.6330     $-$3.20

5     3.9720     4.5892     0.7299     $-$6.98     3.7593     2.2197     0.9620     $-$63.31     4.2449     1.7674     0.6551     $-$5.38

6     3.2195     12.116     0.4100     $-$10.66     4.2598     2.5815     0.4380     $-$65.90     3.9386     4.1249     0.4735     $-$5.52

7     3.0735     16.582     0.4100     $-$11.94     3.8533     5.3273     0.4277     $-$65.67     3.3793     8.3846     0.4188     $-$4.49

8     2.5324     16.508     0.9569     $-$28.45     3.5599     5.5314     0.6935     $-$66.59     2.9959     15.765     0.4111     $-$5.25

9     2.4556     31.778     0.5339     $-$11.20     3.4864     9.5203     0.7144     $-$67.14     2.4976     34.652     0.3970     $-$6.04

10     1.7724     64.592     0.7168     $-$70.50     3.0765     12.411     0.4000     $-$65.41     1.7031     75.996     0.4970     $-$9.94

11     1.5648     88.297     0.4140     $-$9.91     2.8323     22.877     0.4000     $-$65.10     0.60477     181.39     0.8516     $-$19.76

12                     2.9443     23.326     0.5089     $-$66.05

13                     2.9126     34.276     0.6047     $-$66.95

14                     2.5195     58.878     0.6590     $-$66.11

15                     1.9977     118.41     0.8048     $-$69.42

$q^{\prime}_{min} = 0.69$                 $q^{\prime}_{min} = 0.59$                 $q^{\prime}_{min} = 0.53$

3cNGC 4786         3cNGC 4797         3cNGC 5084

2-5 6-9 10-13
1     4.3086     0.3298     0.9381     $-$75.16     4.9673     0.1404     0.5900     $-$43.03     5.2766     0.1000     0.5300     21.32

2     4.4393     0.5740     0.7206     $-$15.72     4.2167     0.4608     0.8116     $-$36.34     4.7053     0.4536     0.8297     87.62

3     4.1755     1.2669     0.8197     $-$18.02     3.6904     0.9806     0.9678     $-$43.52     4.5468     0.8996     0.8877     83.20

4     3.5620     2.7412     0.7348     $-$17.26     3.2996     2.0556     0.9337     $-$79.20     4.4677     1.8366     0.7501     82.48

5     3.4411     4.7184     0.8148     $-$17.94     2.7301     5.4252     0.8406     $-$50.53     4.0750     4.5063     0.5300     82.62

6     2.5653     5.9174     0.8593     $-$22.41     2.7371     7.4784     0.6212     $-$31.61     3.6607     8.4099     0.5300     82.35

7     2.5992     7.7617     0.6900     $-$14.24     2.3562     12.931     0.6168     $-$32.26     3.2294     18.160     0.5300     82.10

8     2.6956     12.806     0.6900     $-$16.27     1.9419     18.686     0.8588     $-$29.80     2.6726     47.305     0.5300     80.74

9     2.3145     14.855     0.9089     $-$16.79     1.3603     33.385     0.9183     $-$58.23     1.5243     105.20     0.5300     80.05

10     2.1400     23.637     0.6900     $-$14.79

11     1.7393     28.847     0.9223     $-$20.55

12     1.3589     55.304     0.6900     $-$17.07

13     1.1908     109.29     0.8174     $-$39.18

$q^{\prime}_{min} = 0.75$                 $q^{\prime}_{min} = 0.31$                 $q^{\prime}_{min} = 0.56$

3cNGC 5193         3cNGC 5208         3cNGC 5838

2-5 6-9 10-13
1     5.2531     0.1000     0.7500     61.39     4.7132     0.1000     0.6486     $-$3.62     5.3228     0.1000     0.7369     $-$30.22

2     4.3480     0.4486     0.7500     62.53     4.7639     0.2374     0.6454     $-$13.42     5.5671     0.1567     0.5600     $-$38.57

3     4.3466     0.9357     0.7500     67.88     4.2336     0.5364     0.6789     $-$13.20     4.9070     0.3916     0.9063     $-$73.40

4     3.9765     1.9775     0.8096     69.80     3.7478     1.3023     0.6911     $-$14.50     4.5814     0.7877     0.8646     $-$37.26

5     3.4464     4.6026     0.7500     73.33     3.9250     1.6661     0.3100     $-$16.77     4.4185     1.9833     0.7133     $-$36.51

6     3.0327     9.7373     0.8422     74.22     3.4542     2.3510     0.5494     $-$15.21     3.9404     4.5365     0.8312     $-$58.05

7     2.5116     18.578     0.9809     83.30     3.5274     5.4322     0.3100     $-$16.49     3.4456     10.621     0.6015     $-$41.87

8     1.6530     47.070     0.9401     64.60     3.0443     6.0382     0.5683     $-$14.00     2.7230     32.183     0.5600     $-$41.58

9                     2.5614     9.2453     0.3100     $-$17.57     2.0282     57.223     0.5600     $-$41.57

10                     2.5448     13.568     0.3100     $-$17.33

11                     2.2904     17.293     0.4537     $-$16.36

12                     1.8148     22.700     0.3100     $-$21.79

13                     1.5476     44.782     0.3775     $-$21.63

14                     0.8870     129.08     0.8838     $-$29.79

MGE Parameters (with a Free PA)

j     $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA      $\log_{10}~I_{H,j}$     $\sigma^{\prime}_{j}$     $q^{\prime}_{j}$     PA

$(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)     $(L_{\sun}~\rm{pc}^{-2})$     (arcsec)          (deg.)

(1)     (2)     (3)     (4)     (5)      (2)     (3)     (4)     (5)     (2)     (3)     (4)     (5)

$q^{\prime}_{min} = 0.38$                 $q^{\prime}_{min} = 0.95$

3cNGC 6861         3cNGC 6958

2-5 6-9
1     5.0115     0.1240     0.5642     $-$19.16     5.3939     0.1000     0.9500     $-$61.22

2     4.8017     0.2137     0.8592     $-$24.84     5.3091     0.1893     0.9500     $-$67.25

3     4.8535     0.5023     0.3800     $-$35.82     4.6775     0.6404     0.9500     $-$70.95

4     4.5926     0.6724     0.9826     $-$25.22     4.2927     1.4673     0.9500     $-$69.61

5     4.2607     1.5157     0.8176     $-$36.10     3.8961     3.0894     0.9500     $-$70.56

6     4.2174     3.4836     0.4681     $-$37.45     3.1265     7.1427     0.9500     $-$70.24

7     3.7291     4.4471     0.7497     $-$38.68     2.7974     14.133     0.9500     $-$77.93

8     3.5860     7.5404     0.4153     $-$37.57     1.9274     33.343     0.9500     $-$89.82

9     3.2756     11.558     0.5075     $-$38.29     1.5789     42.066     0.9500     $-$85.45

10     2.6651     13.092     0.6773     $-$37.38     0.87033     90.065     0.9942     $-$73.57

11     2.6056     23.592     0.4402     $-$39.64

12     2.3448     26.382     0.7077     $-$42.19

13     1.8380     43.342     0.5320     $-$39.05

14     1.6966     46.380     0.8709     $-$84.57

15     1.4120     153.99     0.5012     $-$70.61

Nearly full-frame HST WFC3/F160W mosaics, together with the (dust-masked)  MGE solutions (overplotted in red) that fits a variable PA for all components. At larger radii, some galaxies exhibit highly flattened stellar isophotes, resulting in unavoidable discrepancies. Contours are shown at logarithmic intensity intervals.


# Photonic next-generation reservoir computer based on distributed feedback in
optical fiber

Nicholas Cox<EMAIL_ADDRESS>Joseph Murray  Joseph Hart
Brandon Redding  U.S. Naval Research Laboratory, 4555 Overlook Ave, SW,
Washington, DC 20375, USA

###### Abstract

Reservoir computing (RC) is a machine learning paradigm that excels at
dynamical systems analysis. Photonic RCs, which perform implicit computation
through optical interactions, have attracted increasing attention due to their
potential for low latency predictions. However, most existing photonic RCs
rely on a nonlinear physical cavity to implement system memory, limiting
control over the memory structure and requiring long warm-up times to
eliminate transients. In this work, we resolve these issues by demonstrating a
photonic next-generation reservoir computer (NG-RC) using a fiber optic
platform. Our photonic NG-RC eliminates the need for a cavity by generating
feature vectors directly from nonlinear combinations of the input data with
varying delays. Our approach uses Rayleigh backscattering to produce output
feature vectors by an unconventional nonlinearity resulting from coherent,
interferometric mixing followed by a quadratic readout. Performing linear
optimization on these feature vectors, our photonic NG-RC demonstrates state-
of-the-art performance for the observer (cross-prediction) task applied to the
Rössler, Lorenz, and Kuramoto-Sivashinsky systems. In contrast to digital NG-
RC implementations, this scheme is easily scalable to high-dimensional systems
while maintaining low latency and low power consumption.

> Reservoir computing (RC) is a machine learning technique that has shown
> great success in chaotic systems analysis. For RCs implemented in digital
> electronics, the power consumption and computational latency are prohibitive
> for time and energy-sensitive applications. These limitations can be
> overcome by performing computations in the analog domain using a physical
> system. In particular, photonic RCs are especially promising as they exploit
> the high data encoding rates and parallelism of optical signals. This work
> demonstrates a new photonic reservoir computing architecture that uses
> Rayleigh backscattering in optical fiber as the main computational engine.
> In contrast to previous photonic RCs, which form a reservoir through a
> complicated feedback system, we achieve comparable computational accuracy
> using a simple feedforward approach that implements a virtual reservoir by
> applying a random mixing transformation to a fixed number of consecutive
> time samples at a given time step. We employ this architecture to create a
> time series observer which cross-predicts unknown state variables of a
> multi-dimensional chaotic system from a subset of measurable states. We
> optimize our experimental system using the additional degree of memory
> control provided by the virtual reservoir, achieving state-of-the-art
> performance on a series of benchmark observer tasks.

##  I Introduction

Recurrent neural networks (RNN) are a class of artificial neural networks
especially suited for dynamical systems analysis. RNNs derive their
computational power from the memory generated by an optimized feedback loop
between neuron layers. While this architecture has shown considerable success,
the learning process often suffers from vanishing or exploding gradients [1].
Reservoir computing (RC) is an RNN subclass developed around a recurrent
feedback network connected with fixed random weights. RCs avoid the training
pitfalls of standard RNNs because they require optimization of only a single
output layer that maps the reservoir state to a desired output. The idea of
employing a physical system as a reservoir dates back to the conception of RC
in the form of echo state networks [2] and liquid state machines [3]. Shortly
after, researchers realized RCs using water buckets [4], VLSI chips [5], and
electronic feedback loops [6].

Optics serves as a particularly promising platform for physical reservoirs
because the high modulation bandwidth and multiplexing capability of light
enable high-speed low-power computational inference. Paquot et al. [7] and
Larger et al. [8] developed the first photonic RCs using a delay-line
architecture originally proposed for electronic RCs [6]. This design, which
remains the most commonly investigated approach [9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], uses delayed feedback
following a single nonlinear transformation to define time-multiplexed virtual
nodes. Delay-line RCs have consistently demonstrated state-of-the-art
performance for tasks such as speech classification [29], time series
prediction [30], and distortion compensation for optical communications [18],
but the virtual node structure proves challenging to scale for the processing
of high-dimensional data. Multi-channel methods have been recently explored to
improve scalability [31, 32, 33, 34, 35, 36, 37, 38], but they come at the
cost of increased complexity. An alternative structure involving the coherent
mixing of spatial electromagnetic field profiles in free-space [39, 40, 41,
42, 43, 44] aims to overcome this limit on the maximum number of nodes, but
the computation speed is severely limited by the update rate of spatial light
modulators or digital micromirror devices.

A common feature of both the standard delay-line and free-space designs is the
existence of a nonlinear feedback loop that dictates the evolution of an
internal reservoir state. Performing a read-out operation on this internal
state yields an output feature vector. This feature vector is finally mapped
to the desired output by linear optimization using a training data set. As an
alternative to nonlinear feedback, Refs. [45, 46, 47, 48, 49] demonstrate the
feasibility of high-performance RC using a fully linear optical reservoir
followed by a nonlinear readout (e.g. from the photodetector response).
Recently, Gauthier et al. [50] introduced next generation reservoir computers
(NG-RC) as a further simplification of the linear reservoir concept. An NG-RC
removes the reservoir entirely by explicitly constructing output feature
vectors from various nonlinear combinations of data from current and previous
time steps. Although NG-RCs have shown impressive performance on a variety of
tasks with reduced training and warm-up time [51, 52, 53], the generation of
nonlinear feature vectors is computationally intensive for high-dimensional
systems or systems with high-order nonlinearities. While photonics might be
used to accelerate this process, creating feature vectors of the form proposed
in Ref. [50] on a photonic platform would require a complex optical routing
scheme and the ability to implement a variety of nonlinear transforms.

In this work, we show that state-of-the-art performance on benchmark tasks is
achievable using a different form of feature vector created by applying random
nonlinear projections to the input data. As a proof-of-concept, we present a
photonic NG-RC using a simple fiber optic platform that can efficiently
perform these projections on high-dimensional data with low latency and low
power consumption [54]. The photonic NG-RC operates by encoding an input
feature vector containing current and past data points onto the phase of an
optical pulse, injecting the pulse into standard single-mode fiber, and
recording the Rayleigh backscattering pattern measured by a photodetector. The
net nonlinear transformation, provided by coherent complex-domain mixing
followed by a quadratic readout, is a temporal analog to the one performed in
the spatial domain in Refs. [39, 40, 41, 42]. Comparing our fiber-optic design
to this spatial approach, we find advantages in that (1) our system is not
bottle-necked by the speed limitations of a spatial light modulator and (2)
the NG-RC framework removes the need for nonlinear optical-electronic-optical
feedback. Additionally, when compared to delay line RCs, we find better
scalability to high-dimensional data analysis while maintaining the benefits
of digital NG-RCs: short warm up time, reduced training duration, and
straightforward optimization and control of hyperparameters such as memory
design. Exploiting this control over memory structure, we introduce a novel
mask procedure that transforms elements of the input feature vectors at
different time delays. We find that tailoring this mask can improve
performance of the photonic NG-RC, especially for the analysis of high-
dimensional data.

We illustrate the versatility of the photonic NG-RC by applying it to a range
of time series observer tasks [55]. The role of an observer is to examine a
limited subset of state variables in a multi-dimensional dynamical system and
cross-predict the remaining states [56, 57, 58, 59, 60, 61]. This cross-
prediction task is important for both real-world systems (where some
quantities are unmeasurable but crucial for control or decision-making
purposes) and as a benchmark RC task (since successful cross-predictions for
chaotic systems relies on a delay-coordinate embedding [62] of the attractor
that is highly sensitive to the memory properties). In this work, we test the
photonic NG-RC observer performance on the chaotic Rössler [63], Lorenz [64],
and Kuramoto-Sivashinsky (KS) [65, 66] systems. For each system, we optimize
the memory by varying the number of delay samples present in the input feature
vector as well as the mask function applied to them. We use the Rössler
equations, characterized by the simplest chaotic attractor, to present initial
results and demonstrate the ability to achieve high performance with
relatively few training points. For the Lorenz system, we study the effect of
varying the system sampling rate after completing the initial memory
optimization. Finally, we use the spatiotemporal KS system to highlight the
ability of the fiber NG-RC to analyze high-dimensional data. Finding that the
photonic NG-RC achieves state-of-the-art performance for each task, we then
examine the potential of our design to perform these predictions with low
latency and low power consuption. We compare the power consumption to that of
digital electronic NG-RCs to identify tasks for which the photonic system is
best suited.

##  II Experimental design and operation

The photonic NG-RC takes an input feature vector comprising a fixed number of
consecutive time-delay samples and performs a nonlinear random projection to
generate an output feature vector in higher-dimensional space. This process is
depicted in Fig. 1 for an observer task in which the $x$ component of a
chaotic time series is known and the objective is to cross-predict $y$ and
$z$. An input feature vector is first constructed in software by concatenating
the current value of $x$ with $K-1$ previous time samples. Unlike a
traditional RC, for which the cavity lifetime dictates the system memory, the
NG-RC memory is easily controlled by adjusting the number of time steps
included in the input feature vector. In this work, we introduce a memory mask
function, $m_{k}$, to scale the “Observed Data” (Fig. 1(b)) in the input
feature vector. The mask provides an additional degree of control in the
photonic NG-RC, which we later find to be essential for optimizing the system
performance. The input feature vector, $\mathbf{v}_{\text{in},n}$, at the
$n^{\text{th}}$ time step is defined to be

$\mathbf{v}_{\text{in},n}=\mathbf{x}_{n}\oplus
m_{1}\mathbf{x}_{n-1}\oplus\cdots\oplus m_{K-1}x_{n-(K-1)}.$ (1)

The symbol $\oplus$ is the direct sum (concatenation) operator,
$\mathbf{x}_{n}$ is a vector composed of each observed state variable at
$n^{\text{th}}$ time step (in the case shown in Fig. 1 with a single observed
variable, $\mathbf{x}_{n}=x_{n}$), and $m_{k}$ is the multiplicative mask
applied to the NG-RC memory. Here, we considered three masks: (1) the unit
mask $m_{k}=1$, (2) the decaying mask $m_{k}=\beta^{k}$ (with $\beta=0.8$),
and (3) the random mask comprising pseudorandom numbers between 0 and 1. In
each case, the mask function remains fixed for each sample step $n$.

Figure 1: An overview of the photonic NG-RC system, detailing (a) the problem
definition, (b) the construction of input feature vectors, (c) the
experimental fiber-optic setup, and (d) the final prediction step in digital
electronics.

As shown in Fig. 1(c), the photonic NG-RC encodes input vectors into the phase
of an optical carrier using an electro-optic phase modulator ($\phi$-EOM). We
used phase modulation because it allows for a more natural encoding of
positive and negative values compared to amplitude modulation which requires a
voltage offset to accommodate negative values. An amplitude modulator (an
acousto-optic modulator, AOM, was used in this work) then carves a pulse
surrounding the phase modulated region. Note that the pulse duration was
chosen to be slightly longer than the phase-modulated region to provide a
fixed phase reference in each pulse. This encoded pulse was injected into a
stretch of standard single mode fiber and the Rayleigh backscattered (RBS)
light was measured by a photodetector. The resulting voltage, proportional to
the RBS pattern irradiance, was finally digitized using an analog-to-digital
converter (ADC) to yield the output feature vector. The length of this output
feature vector is proportional to the scattering fiber length, providing our
design with a simple method of tailoring the output dimensionality to a
specific problem.

We can model the above procedure to understand the nonlinear random projection
performed by the photonic NG-RC. To begin, input feature vectors are encoded
into the phase modulator driving voltage

$V_{\text{in},n}(t)\approx\sum_{\mu}v_{\text{in},n}^{\mu}\text{rect}\left(\frac{t-\mu\tau}{\tau}\right),$
(2)

where $v^{\mu}_{\text{in},n}$ is the $\mu$ component of
$\mathbf{v}_{\text{in},n}$ and $\tau$ is the duration of a single rectangular
pulse encoding an input element. After the AOM defines the envelope $|E(t)|$,
the modulated optical signal passes into the scattering fiber. Measuring the
subsequent RBS irradiance using a photodector yields the voltage

$V_{\text{out},n}(t)\approx\left||E(t)|e^{i(\omega t+\pi
V_{\text{in},n}(t)/V_{\pi})}\ast r_{\text{RBS}}(t)\right|^{2},$ (3)

where $\omega$ is the optical angular frequency, $V_{\pi}$ is the modulator
voltage to impart a phase shift of $\pi$, $(*)$ is the convolution operator,
and $r_{\text{RBS}}(t)$ is the Rayleigh backscattering response of the optical
fiber (See SI for a full mathematical description). The output vector is
finally formed by sampling Eq. 3 at a fixed rate $t_{s}$:

$v^{\mu}_{\text{out},n}\approx V_{\text{out},n}(\mu t_{s}).$ (4)

We finally see that the overall transformation
$\mathbf{v}_{\text{in},n}\rightarrow\mathbf{v}_{\text{out},n}$ is nonlinear
due to the voltage relationship in Eq. 3 including a convolution followed by
quadratic readout.

As in a traditional reservoir observer, the NG-RC operation assumes there is
an initial training period during which all state variables of a dynamical
system are available. This training period is used to optimize a output weight
matrix ($W$) containing one row for each predicted variable. After the
training period, the output weights are held fixed and used to predict the
unknown state variables ($y$,$z$ in Fig. 1) during the testing phase (e.g.
$y_{n}=W_{y}\mathbf{v}_{\text{out},n}+c_{y}$, where $c_{y}$ is a DC offset).
In this work, the prediction stage was performed via digital electronics after
digitizing the output feature vector. In the future, the weight matrix could
be applied in the analog photonic domain using the technique developed in
Refs. [14, 67] to minimize the system latency.

Experimentally, we constructed the photonic NG-RC using the basic architecture
shown in Fig. 1(c). A narrow-linewidth ($\sim$kHz) seed laser operating at a
wavelength of $\sim$1550 nm is directed to an electro-optic phase modulator to
encode the input feature vector. Using an arbitrary waveform generator (AWG),
each element of the input feature vector $\text{v}_{\text{in},n}$ is encoded
into the peak voltage of a $\tau=10$ ns rectangular pulse. These $\phi$-EOM
input voltages are scaled so that the maximum and minimum values over all
input vectors are $\pm V_{\pi}$. The length of the phase encoded input vector
is $KQ\tau$, where $K$ is the memory length (i.e. the number of time steps
included in the input feature vector) and $Q$ is the number of observed input
dimensions ($Q=1$ in the example shown in Fig. 1). The phase-encoded light is
then directed to an AOM which forms an amplitude envelope of width
$(KQ+6)\tau$ centered on the phase-encoded signal. The buffer time of 6$\tau$
allows for a portion of the unmodulated source to act as a zero phase-shift
reference (experimentally, we found that the NG-RC performance is relatively
insensitive to the length of this buffer within the range of 2$\tau$ to
20$\tau$). The pulse was then amplified using an erbium doped fiber amplifier
(EDFA) and coupled through a circulator to a section of single mode fiber of
length $L=500$ m. The RBS light was directed by the circulator to a second
EDFA before reaching a 120 MHz bandwidth photodetector. The detected signal
was digitized at 1 GS/s. In this proof-of-concept demonstration, the first
EDFA was used to compensate for insertion loss through the EOM and AOM, while
the second EDFA was used to mitigate the effect of detector and digitizer
noise. In the future, using a higher power seed laser, enhanced Rayleigh
backscattering fiber, lower insertion loss modulators, and a lower noise
detector could obviate the need for the EDFAs.

The total duration of the RBS pattern, and hence the output feature vector,
was $\tau_{\text{RBS}}=2L/(c/n_{g})+KQ\tau\approx 5$ $\mu$s, where $c$ is the
vacuum speed of light and $n_{g}$ is the group index. This entire process was
repeated for each input time step in a data set at a repetition period of
$T=6$ $\mu$s (set to be slightly longer than $\tau_{\text{RBS}}$). In each
experiment presented below, we sent a total of $N$ input feature vectors
including information from $N$ consecutive samples of a dynamical system.
After retrieving all output feature vectors, we partitioned into a set of
$N_{\text{train}}$ training points followed by $N_{\text{test}}$ testing
points. The output weight matrix ($W$) was then computed offline by ridge
regression and used to cross-predict output variables in the test phase.

To mitigate experimental noise, we repeated the measurements 25 times
(Rössler, Lorenz) or 5 times (KS) and averaged the resulting RBS patterns to
obtain a lower noise output feature vector (see SI for details on the effect
of averaging on the prediction accuracy and measurement SNR). Although the RBS
pattern was sampled at 1 GHz , the temporal correlation width of the RBS light
is set by the $\tau=10$ ns encoding bandwidth [68]. In this work, we performed
a 5-point moving mean filter and down-sampled the measured pattern at half the
correlation width, or a rate of 200 MS/s, retaining some partially correlated
information. This resulted in an output feature vector of length
$D\approx\tau_{\text{RBS}}/(\tau/2)=1000$.

##  III Rössler System

We first applied the photonic NG-RC to analyze the Rössler system, which is a
classic chaotic system often used to benchmark RC performance. The Rössler
system is governed by the coupled differential equations [63]

$\displaystyle x_{t}$ $\displaystyle=-y-z$ $\displaystyle y_{t}$
$\displaystyle=x+ay$ $\displaystyle z_{t}$ $\displaystyle=b+z(x-c),$ (5)

where the subscript $t$ denotes a derivative with respect to time. We use the
standard parameters $a=0.5$, $b=2.0$ and $c=4.0$, which are known to exhibit
chaos [55]. In this section, we train the fiber NG-RC to use $x$ as an input
to cross-predict the variables $y$ and $z$. Ground truth data was generated by
randomly selecting a set of initial conditions and integrating Section III
with a step size $\Delta t_{\text{R\"{o}ssler}}=0.125$ using the fourth order
Runge-Kutta method. A typical input feature vector with memory length $K=30$
and a random input mask is shown in Fig. 2a, and the resulting output feature
vector produced by the fiber optic platform appears in Fig. 2b. This process
was repeated for $N=3000$ time steps, of which the first 1000 are shown in
Fig. 2c.

Figure 2: Experimental data for the Rössler observer. a) the phase modulator
input voltage that encodes the input vector $\mathbf{v}_{\text{in},n}$ at
$n=800$ ($t_{\text{R\"{o}ssler}}=100$). b) the voltage output proportional to
the RBS irradiance that is downsampled to form the output vector
$\mathbf{v}_{\text{out},n}$ at $n=800$. c) a plot of the first 1000 retrieved
speckle patterns, with the dotted blue line marking the location of the
pattern shown in part b.

We allocated the first $N_{\text{train}}=1000$ points (corresponding to
$t_{\text{R\"{o}ssler}}=N_{\text{train}}\Delta t_{\text{R\"{o}ssler}}=125$)
for training and used ridge regression to calculate the weight matrix that
minimized errors in the cross-predicted values of $y$ and $z$. We then tested
the system performance on the remaining $N_{\text{test}}=2000$ data points. As
shown in Fig. 3, the NG-RC accurately predicted the $y$ and $z$ coordinates of
the Rössler system after the training period ($t_{\text{R\"{o}ssler}}>125$)
using an input feature vector with memory of length $K=30$ and a random memory
mask.

Figure 3: Output of the Photonic NG-RC observer (red dashed) with $K=30$ and a
random mask compared to ground truth (blue solid) for $N=3000$ samples of the
Rössler system. The grey shaded region marks the training set of
$N_{\text{train}}=1000$ points following a $K-1=29$-step warm-up time. The
errors for the data shown are $\text{NRMSE}_{y}=1.23\times 10^{-2}$,
$\text{NRMSE}_{z}=1.89\times 10^{-2}$.

To quantitatively evaluate the photonic NG-RC performance, we computed the
normalized root mean square error (NRMSE) for the variable $y$ (and similarly
for $z$) during the testing phase ($n>N_{\text{train}}$), as

$\text{NRMSE}_{y}=\frac{1}{\sigma}\sqrt{\sum_{n}\left(y_{n}-\hat{y}_{n}\right)^{2}},$
(6)

where $\hat{y}_{n}$ is the RC prediction and $y_{n}$ is the ground-truth value
at time step $n$. The error is normalized by the standard deviation, $\sigma$,
of the ground truth curve. For ease of comparison, we averaged the NRMSE for
outputs y and z, defining $\text{NRMSE}=(\text{NRMSE}_{y}+NRMSE_{z})/2$. For
$K=30$ case shown in Fig. 3, we found a NRMSE of $1.56\times 10^{-2}$
($\text{NRMSE}_{y}=1.23\times 10^{-2}$, $\text{NRMSE}_{z}=1.89\times
10^{-2}$).

We then studied the effect of the memory length and mask type on the photonic
NG-RC performance. The Rössler observer task was repeated with 5 different
initial conditions, formed by randomly choosing $x_{1}$, $y_{1}$, and $z_{1}$
between 0 to 1 and integrating over 5000 time steps to ensure the system has
settled into the chaotic attractor. The average NRMSE is shown in Fig. 4(a) as
a function of memory length, $K$, for either a unit, decaying, or random mask.
For further insight, the top axis displays $t_{\text{R\"{o}ssler},K}=K\Delta
t_{\text{R\"{o}ssler}}$ to indicate memory duration in the system time units.
The top and bottom ranges of the error bar mark the maximum and minimum NRMSE
from the 5 runs, and the data point marks their mean. For each test, we used
$N_{\text{train}}=1000$ and $N_{\text{test}}=2000$. For $K<20$
($t_{\text{R\"{o}ssler},K}<2.5$), the unit memory outperforms both the random
and decaying masks. The decaying mask ($\beta=0.8$) experiences an increase in
NRMSE for $K>8$ ($t_{\text{R\"{o}ssler},K}>1$), likely due to a memory-
limiting effect as the phase change encoding past history inputs falls below
the noise floor. The random mask begins to slightly outperform the unit mask
for $K>20$ ($t_{\text{R\"{o}ssler},K}>2.5$), reaching a minimum mean NRMSE of
$1.71\times 10^{-2}$ ($\text{NRMSE}_{y}=1.51\times 10^{-2}$,
$\text{NRMSE}_{z}=1.94\times 10^{-2}$) at $K=30$
($t_{\text{R\"{o}ssler},K}=3.75$). This result improves on the state-of-the-
art for Rössler observers based on physical RC systems
($\text{NRMSE}_{y}=3.6\times 10^{-2},\ \text{NRMSE}_{z}=0.1$) [69]. While all-
digital RCs have reported lower NRMSE for noise-free data, our result is
comparable to the NRMSE obtained in a digital RC with simulated additive noise
at an signal-to-noise ratio (SNR) of 40 dB [55]. (See SI for more information
on photonic NG-RC performance versus SNR). Note that all-digital NG-RCs have
not been applied to the Rössler task, to the best of our knowledge, so we
cannot make a direct performance comparison. We are able to compare more
closely in the next section in which the Lorenz system is investigated.

Figure 4: a) NRMSE versus memory length $K$/duration
$t_{\text{R\"{o}ssler},K}$ and mask type for the Rössler observer. b) NRMSE
for the Rössler system versus number of training points ($N_{\text{train}}$)
using a memory with $K=30$ and a random mask.

In addition to high prediction accuracy, the photonic NG-RC requires a
dramatically shorter warm-up time than traditional physical or digital RCs.
The photonic NG-RC only requires a warm-up time equal to $K$, the number of
memory elements in the input vector. The NG-RC optimized for the Rössler
system requires only 30 time steps to warm up, compared to physical RCs which
often require $10^{3}$ to $10^{5}$ steps [70, 71]. Another major advantage of
NG-RCs is that they require fewer training points than conventional RCs. We
explore this in depth by studying the NRMSE versus the number of training
points $N_{\text{train}}$. The results are shown in Fig. 4(b) for a fixed
memory of length $K=30$ and a random memory mask. We begin to see satisfactory
performance for $N_{\text{train}}$ as low as 200, finding
$\text{NRMSE}=5.3\times 10^{-2}$. The errors drop as $N_{\text{train}}$
increases until the improvement tapers off at $N_{\text{train}}=500$
($\text{NRMSE}=2.6\times 10^{-2}$). At this point the training duration
corresponds to a time interval $t_{\text{R\"{o}ssler}}=62.5$, compared to
$t_{\text{R\"{o}ssler}}=260$ training period for the simulated RC in Ref.
[55].

##  IV Lorenz System

Having found good results the Rössler observer, it is important to verify that
the NG-RC applies more generally to a diverse set of problems. In this
section, we examine the Lorenz system

$\displaystyle x_{t}$ $\displaystyle=a(y-x)$ $\displaystyle y_{t}$
$\displaystyle=x(b-z)-y$ $\displaystyle z_{t}$ $\displaystyle=xy-cz,$ (7)

with the standard parameters $a=10$, $b=28$ and $c=8/3$. This system is also
chaotic and characterized by a double-lobe attractor. We again employ the
photonic NG-RC as an observer that cross-predicts the $y$ and $z$ coordinates
from measurements of $x$. We generate time series by fourth-order Runge-Kutta
integration with a step size of $\Delta t_{\text{Lorenz}}=0.025$. As in the
previous section, we begin with a demonstration of results for a well-tuned
photonic NG-RC and follow up with the memory hyperparameter investigation that
informed the optimization. Finally, we repeat the experiment on time-series
data with doubled step size $\Delta t_{\text{Lorenz}}=0.05$ to separate the
performance impact of changing $K$ and $t_{\text{Lorenz},K}$. Lorenz cross-
prediction results for the parameters $\Delta t_{\text{Lorenz}}=0.025$,
$N_{\text{train}}=1000$, and $N_{\text{test}}=2000$ are shown in Fig. 5. The
memory is defined to have $K=30$ ($t_{\text{Lorenz},K}=0.75$) and a random
mask, and the resulting errors are $\text{NRMSE}_{y}=1.99\times 10^{-2}$ and
$\text{NRMSE}_{z}=2.47\times 10^{-2}$.

We studied the prediction accuracy as a function of memory length and mask
type, as shown in Fig. 6(a). The error bars are generated the same way as in
the previous section. One point of difference with the Rössler system is the
appearance of a larger separation between the error obtained using unit and
random memories. However, the hyperparameter optimization of the Lorenz system
shares key features with the Rössler equations: the unit mask incurred the
smallest error for short memories ($K\leq 20$, $t_{\text{Lorenz},K}<0.5$), the
best results occurred for a random mask with $K=30$
($\text{NRMSE}_{y}=2.04\times 10^{-2}$, $\text{NRMSE}_{z}=2.77\times
10^{-2}$), and the worst performance was found for the decaying mask (with
$\beta=0.8$).

Figure 5: Output of the Photonic NG-RC observer (red dashed) with $K=30$ and a
random mask compared to ground truth (blue solid) for $N=3000$ samples of the
Lorenz system. The grey shaded region marks the training set of
$N_{\text{train}}=1000$ points following a $K-1=29$-step warm-up time. The
errors for the data shown are $\text{NRMSE}_{y}=2.04\times 10^{-2}$,
$\text{NRMSE}_{z}=2.77\times 10^{-2}$.

Figure 6: a) Lorenz obsever NRMSE versus memory length $K$/duration
$t_{\text{Lorenz},K}$ and mask type, with system step size $\Delta
t_{\text{Lorenz}}=0.025$. b) A repeat of part a with doubled step size $\Delta
t_{\text{Lorenz}}=0.05$

.

Next, we double the time series step size to $\Delta t_{\text{Lorenz}}=0.05$.
As seen in Fig. 6(b), the required memory length for peak performance is
halved to $K=15$. Since the optimal memory extent $t_{\text{Lorenz},K}=0.75$
is unchanged, this suggests that the memory should be optimized in terms of
system time. In contrast to Fig. 6(a), the smallest mean NRMSE of $2.00\times
10^{-2}$ ($\text{NRMSE}_{y}=1.61\times 10^{-2},\text{NRMSE}_{z}=2.38\times
10^{-2}$) occurred for a unit memory mask, improving slightly upon the $K=30$
random mask at step size $\Delta t_{\text{Lorenz}}=0.05$. For an overall
performance comparison, the digital NG-RC in Ref. [50] achieved
$\text{NRMSE}_{z}=1.75\times 10^{-2}$ for a Lorenz observer using variables
$x$ and $y$ to predict $z$. The results presented here are comparable despite
using only $x$ values as input. Finally, we note that the errors for a random
memory mask drop below those of the unit mask at $K\geq 25$, indicating that
the total number of points in the input feature vector may dictate the optimal
choice for the mask function. This suggests that systems with higher
dimensionality, even if only short memories are required, may benefit the most
from the application of a random weighting of the input vector elements.

##  V Kuramoto-Sivashinsky system

In the previous two sections, we showed that the photonic NG-RC works well for
low-dimensional systems. However, one of its largest strengths is the ability
to scale to the analysis of high-dimensional systems without a significant
increase in power consumption or computational cost. To illustrate this
capability, we tested the performance of the fiber NG-RC as an observer for
the spatiotemporally chaotic Kuramoto-Sivashinsky (KS) system defined by the
partial differential equation

$y_{t}=-yy_{x}-y_{xx}-y_{xxxx}.$ (8)

The data consist of $P=64$ evenly spaced grid points between spatial
coordinates $x=0$ and $x=22$. Choosing $Q$ equidistant spatial samples to be
continuously available, the observer interpolates the profile by cross-
predicting the other $P-Q$ elements. The ground truth KS data is found by
integrating Eq. 8 with the Crank-Nicholson method using periodic boundary
conditions and a time step $\Delta t_{\text{KS}}=0.25$. In this case, the
NRMSE is defined as [70]

$\text{NRMSE}=\sqrt{\frac{\sum_{s,n}\left(y_{s,n}-\hat{y}_{s,n}\right)^{2}}{\sum_{s,n}y_{s,n}^{2}}},$
(9)

computed over all time ($n$) and space ($s$) indices for cross-predicted
points during the test phase $n>N_{\text{train}}$. All error bars are computed
with the same procedure as in the Rössler and Lorenz systems.

Figure 7: Reconstruction of the spatiotemporal pattern for the KS system using
a photonic NG-RC with $K=10$ and a random mask. The input dimension is $Q=4$,
with the sample points (8,24,40,56) coinciding with the tick marks on the
panel $y$ axes. The top panel shows the reconstruction, the middle contains
the ground truth, and the bottom shows the relative error.

Fig. 7 provides an example of an interpolated profile for $Q=4$ observer
positions, memory of length $K=10$ ($t_{\text{KS},K}=2.5$), training length
$N_{\text{train}}=15000$ ($t_{\text{KS}}=3750$) and testing on the remaining
$N_{\text{test}}=5000$ time steps (up to $t_{\text{KS}}=5000$). The top
subfigure contains the spatiotemporal pattern recovered by the fiber NG-RC,
the middle panel shows the ground truth, and the bottom panel shows the error.
We found that the photonic NG-RC successfully reproduced the general
spatiotemporal pattern with $\text{NRMSE}=0.30$.

Next, we increased the input dimensionality to $Q=8$ and repeated the
experiment with $N_{\text{train}}=15000$ and $N_{\text{test}}=5000$ while
adjusting the memory. As shown in Fig. 8(a), the error initially decreased
with memory length, reaching a minimum $\text{NRMSE}=5.04\times 10^{-2}$ for a
random mask with $K=10$ ($t_{\text{KS},K}=2.5$). This value represents the
state-of-the-art for RC observers, improving upon the all-digital results in
Ref. [55] while reducing training requirements by a factor of 4. In contrast
to the Rössler and Lorenz systems, we find that both the random and decaying
masks lead to lower error than the unit mask for all but the shortest memory
($K=2$, $t_{\text{KS},K}=0.5$). We make sense of this outcome based on the
results in figure Fig. 6(b): the most important factor determining the optimal
memory mask is the total number of points in the input feature vector. For
example, at $Q=8$ and $K=4$, the input feature vector length is 32. For the
Rössler and Lorenz system, the random mask outperformed the unit mask for
feature vectors with $K\geq 25$ ($Q=1$), or a total input feature vector size
of 25. Thus, for high-dimensional systems such as the KS system, the random
mask provides superior performance. We also note that the random mask was
capable of optimal or near-optimal performance for each of the systems studied
here, making it a good initial option when applying the photonic NG-RC to new
tasks.

Figure 8: a) NRMSE versus memory length ($K$)/duration ($t_{\text{KS},K}$) and
mask function for the KS observer with $Q=8$ observation points. b) NRMSE
versus $Q$ for the RC compared to a direct ridge regression on input vectors.

Finally, we measured the NRMSE versus input dimension $Q$ using the optimized
memory of length $K=10$ and a random mask. To evaluate the importance of
generating output feature vectors by a nonlinear transformation, we also
attempted to predict the KS pattern by applying a ridge regression directly to
the input feature vectors, $\mathbf{v}_{\text{in},n}$. As shown in Fig. 8(b),
the NG-RC outperformed this linear optimization for each value of $Q$,
confirming the importance of applying a nonlinear transform to the input data.
In comparison to results reported for digital RCs on the KS task [55], the
photonic NG-RC achieved similar NRMSE for $Q\leq 6$ and improved on the state-
of-the-art for $Q\geq 7$. Again, these results were obtained using
significantly less training data and without the long warm-up time required in
traditional RCs.

##  VI Energy consumption and latency

In addition to providing state-of-the-art predictive performance, the photonic
NG-RC has the potential for significantly lower latency than digital
electronic (software-based) NG-RCs. Accounting for the averaging required to
increase SNR, the photonic NG-RC presented in this work requires $\sim$125
$\mu$s (Rössler, Lorenz) or $\sim$25 $\mu$s (KS) to map an input feature
vector to an output feature vector. We may improve this transformation time by
implementing two changes. The first is to increase the signal power to
eliminate the need for averaging. The second is to increase the encoding rate
from the 100 MHz used in this work to 10 GHz, a typical modulation rate in
optical fiber. As a result, we reduce the time to generate output feature
vectors to 5 ns using a 5 m fiber. To eliminate the latency due to output
vector digitization and software-based weighting, we may apply output weights
instead by the analog scheme from Refs. [14, 67] to complete inference with a
total latency of 5 ns. The photonic NG-RC could also support multiplicative
increases in speed through spatial or wavelength multiplexing techniques,
providing the potential for exceptionally low latency.

We also analyzed the energy consumption required by the photonic NG-RC to map
an input feature vector to an output feature vector. In particular, we
considered the power consumption required by the laser, the modulators, the
photodetector, the digital to analog converter (DAC) used to drive the
modulator and the analog to digital converter (ADC) used to record the RBS
pattern. Details of the energy consumption calculation are provided in the SI.
We then compared the photonic NG-RC to the energy consumption required by a
digital electronic NG-RC [50], which constructs output feature vectors by
applying nonlinear transformations in software. For this comparison, we first
calculated the number of multiplication operations required to generate a
single output feature vector. We compute this number as a function of the
dimensionality of the observed data, $Q$, the memory length, $K$, and the
highest monomial order of the output feature vector, $d$. The number of
multiplication operations required to generate feature vectors of the form
described in Ref. [50] can be expressed as:

$M_{\text{\text{NG-
RC}}}=\sum_{p=2}^{d}\binom{KQ+p-1}{p}-\binom{(K-1)Q+p-1}{p},$ (10)

where the parentheses denote the binomial coefficient operator. The detailed
derivation of $M_{\text{NG-RC}}$ is provided in the SI. Assuming that each
multiplication required $\approx$ 0.5 pJ (which is typical of state-of-the-art
GPUs, as well as other digital electronic computing architectures [72]), the
energy required to generate the output feature vector at each time step can be
expressed as $E_{\text{NG-RC}}=M_{\text{NG-RC}}\times 0.5$ pJ. Note that this
analysis can be considered a conservative estimate because it did not account
for the power required to store and access memory. In contrast, the photonic
NG-RC effectively performs “in-memory” computing and does not need to store
intermediate values to compute the output from a given input vector.

To compare the energy consumption of the photonic NG-RC with that of a
monomial-based digital-electronic NG-RC, we first note that the nonlinear
transformation imposed by the photonic system is different from the explicit
construction of monomials used to form output feature vectors in the digital
approach. However, using the results of this work, we assume that the
performance of a monomial-based NG-RC is equivalent to a photonic NG-RC with
output feature vector length equal to 10 times the input vector length. Figure
9 compares the energy consumed by equivalent-performance photonic and
software-based NG-RCs as a function of the input feature vector length $(KQ)$
for a fixed memory of length $K=10$. In the photonic approach, we found that
the laser is the primary source of energy consumption for $Q<100$. The ADC
dominates at larger $Q$ as the number of required operations increases with
the size of the output feature vectors. For the digital-electronic NG-RC, the
energy requirements increase with $Q$ in all cases. Repeating for digital NG-
RCs with highest order monomial $d=2$, 3, or 4, we find that increasing $d$
increases the slope in energy consumption as a function of $Q$. We find that
the monomial-based NG-RC is more energy efficient than the photonic NG-RC For
low-dimensional systems with modest nonlinearity (which can be captured using
$d=2$ or 3). However, as the input dimensionality increases, particularly in
systems requiring higher-order nonlinearities, the photonic NG-RC provides the
potential for dramatically lower energy consumption. For example, with input
dimension $Q=100$ and $d=3$, the photonic NG-RC could generate an output
feature vector using 1000$\times$ less energy than the software-based
implementation. It is possible to further reduce the energy requirements by
applying analog weights at the prediction phase [14, 67], so that only a
single ADC operation would be required to obtain a prediction.

Figure 9: Energy required to generate an output feature vector,
$\mathbf{v}_{\text{out},n}$, for the photonic NG-RC (blue lines) and the
digital electronic NG-RC described in Ref. [50] (orange lines). The individual
blue curves indicate the energy required by varying components and the dark
blue line labelled “Photonic NG-RC” represents the total energy consumption.
The memory and DAC operations have identical energy contribution curves. The
photonic scheme requires lower energy for systems with large $Q$ or high non-
linearity (requiring large $d$).

##  VII Discussion

In this work, we introduced a photonic NG-RC that uses Rayleigh backscattering
to implement a nonlinear map between input and output feature vectors using a
combination of coherent, interferometric mixing and a quadratic readout. We
showed that this type of nonlinear projection enables excellent performance on
a variety of benchmark tasks while maintaining the advantages of digital
electronic NG-RCs: short warm-up time, reduced training time, and simple
hyper-parameter optimization. The photonic NG-RC also addresses the
scalability limitations of digital electronic NG-RCs. Since the RBS process is
entirely passive, the photonic NG-RC is capable of performing the equivalent
of $KQD$ multiply and accumulate (MAC) operations while requiring only $KQ$
modulations (to encode the input feature vectors). As a result, the photonic
NG-RC has the potential for dramatically lower latency and power consumption
than a digital-electronic NG-RC when analyzing high-dimensional and/or highly
nonlinear systems.

There are also several areas for improvement. First, in this demonstration, we
used a data encoding rate of 100 MHz, leading to latency of 5 $\mu$s in the
generation of output feature vectors. The effective latency was further
increased by a factor of 5 (KS) or 25 (Rössler, Lorenz) due to the averaging
employed to increase SNR. In the future, higher-speed operation could enable
similar computations in as little as 5 ns. We may also use spatial or
wavelength multiplexing to further accelerate the photonic NG-RC. Second, NG-
RC performance depends strongly on the ability of the nonlinear transformation
to capture the physics of the dynamical system being studied. Fortunately,
additional nonlinearities are easily accessible in this platform with minimal
impact on power consumption or latency. Highly nonlinear fiber could be used
to access all-optical nonlinearities in the RBS fiber while maintaining modest
power levels. It is also possible to introduce additional nonlinearity in the
readout stage by choosing a photo-detector with lower saturation power to
clamp the measured RBS signal. We expect that the optimal transformation will
depend strongly on the target task. Finally, in this demonstration, the entire
input feature vector was encoded on an optical carrier at each time step. In
the future, delayed copies of the observed data could be created using a
series of delay lines, requiring the observed data at a given time step to be
encoded only once. Although Fig. 9 showed that the re-encoding cost is
relatively modest (see the energy consumed by the DAC and modulator), such an
approach could enable the photonic NG-RC to operate directly on data
transmitted over fiber in the analog domain without requiring intermediate
digitization and storage.

In summary, this work showed that a photonic platform can efficiently apply
the nonlinear transformation required for NG-RC simply by performing random
nonlinear projections. While the fiber backscattering architecture presented
here has many attractive features, we expect that this work may spur the
development of alternative photonic NG-RC platforms tailored to different
applications that take advantage of the design flexibility afforded by
foregoing the need to construct a physical cavity.

## Supplementary Material

See the supplementary material for a mathematical model of the photonic NG-RC
system, an analysis on the effect of averaging on performance, and a
derivation of the power consumption for photonic and digital electronic NG-
RCs.

###### Acknowledgements.

The authors acknowledge support from the U.S. Naval Research Laboratory (6.1
Base Funding).

## Conflicts of Interest

The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Rumelhart, Hinton, and Williams [1986] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature 323, 533–536 (1986).
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, “Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication,” Science 304, 78–80 (2004).
  * Maass, Natschläger, and Markram [2002] W. Maass, T. Natschläger, and H. Markram, “Real-Time Computing Without Stable States: A New Framework for Neural Computation Based on Perturbations,” Neural Computation 14, 2531–2560 (2002).
  * Fernando and Sojakka [2003] C. Fernando and S. Sojakka, “Pattern Recognition in a Bucket,” in  _Advances in Artificial Life_, Vol. 2801, edited by G. Goos, J. Hartmanis, J. Van Leeuwen, W. Banzhaf, J. Ziegler, T. Christaller, P. Dittrich, and J. T. Kim (Springer Berlin Heidelberg, Berlin, Heidelberg, 2003) pp. 588–597.
  * Schürmann, Meier, and Schemmel [2004] F. Schürmann, K. Meier, and J. Schemmel, “Edge of Chaos Computation in Mixed-Mode VLSI - A Hard Liquid,” in  _Advances in Neural Information Processing Systems_ , Vol. 17, edited by L. Saul, Y. Weiss, and L. Bottou (MIT Press, 2004).
  * Appeltant  _et al._ [2011] L. Appeltant, M. Soriano, G. Van Der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. Mirasso, and I. Fischer, “Information processing using a single dynamical node as complex system,” Nature Communications 2, 468 (2011).
  * Paquot  _et al._ [2012] Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar, “Optoelectronic Reservoir Computing,” Scientific Reports 2, 287 (2012).
  * Larger  _et al._ [2012] L. Larger, M. C. Soriano, D. Brunner, L. Appeltant, J. M. Gutierrez, L. Pesquera, C. R. Mirasso, and I. Fischer, “Photonic information processing beyond Turing: An optoelectronic implementation of reservoir computing,” Optics Express 20, 3241 (2012).
  * Chembo [2020] Y. K. Chembo, “Machine learning based on reservoir computing with time-delayed optoelectronic and photonic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 013111 (2020).
  * Martinenghi  _et al._ [2012] R. Martinenghi, S. Rybalko, M. Jacquot, Y. K. Chembo, and L. Larger, “Photonic Nonlinear Transient Computing with Multiple-Delay Wavelength Dynamics,” Physical Review Letters 108, 244101 (2012).
  * Brunner  _et al._ [2018] D. Brunner, B. Penkovsky, B. A. Marquez, M. Jacquot, I. Fischer, and L. Larger, “Tutorial: Photonic neural networks in delay systems,” Journal of Applied Physics 124, 152004 (2018).
  * Hart  _et al._ [2019] J. D. Hart, L. Larger, T. E. Murphy, and R. Roy, “Delayed dynamical systems: Networks, chimeras and reservoir computing,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20180123 (2019).
  * Duport  _et al._ [2016a] F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, “Virtualization of a Photonic Reservoir Computer,” Journal of Lightwave Technology 34, 2085–2091 (2016a).
  * Duport  _et al._ [2016b] F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, “Fully analogue photonic reservoir computer,” Scientific Reports 6, 22381 (2016b).
  * Dejonckheere  _et al._ [2014] A. Dejonckheere, F. Duport, A. Smerieri, L. Fang, J.-L. Oudar, M. Haelterman, and S. Massar, “All-optical reservoir computer based on saturation of absorption,” Optics Express 22, 10868 (2014).
  * Carroll and Hart [2022] T. L. Carroll and J. D. Hart, “Time shifts to reduce the size of reservoir computers,” Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 083122 (2022).
  * Vatin, Rontani, and Sciamanna [2019] J. Vatin, D. Rontani, and M. Sciamanna, “Experimental reservoir computing using VCSEL polarization dynamics,” Optics Express 27, 18579 (2019).
  * Argyris, Bueno, and Fischer [2018] A. Argyris, J. Bueno, and I. Fischer, “Photonic machine learning implementation for signal recovery in optical communications,” Scientific Reports 8, 8487 (2018).
  * Nguimdo  _et al._ [2017] R. M. Nguimdo, E. Lacot, O. Jacquin, O. Hugon, G. Van Der Sande, and H. Guillet De Chatellus, “Prediction performance of reservoir computing systems based on a diode-pumped erbium-doped microchip laser subject to optical feedback,” Optics Letters 42, 375 (2017).
  * Soriano  _et al._ [2013] M. C. Soriano, S. Ortín, D. Brunner, L. Larger, C. R. Mirasso, I. Fischer, and L. Pesquera, “Optoelectronic reservoir computing: Tackling noise-induced performance degradation,” Optics Express 21, 12 (2013).
  * Ortín  _et al._ [2015] S. Ortín, M. C. Soriano, L. Pesquera, D. Brunner, D. San-Martín, I. Fischer, C. R. Mirasso, and J. M. Gutiérrez, “A Unified Framework for Reservoir Computing and Extreme Learning Machines based on a Single Time-delayed Neuron,” Scientific Reports 5, 14945 (2015).
  * Antonik  _et al._ [2017] P. Antonik, F. Duport, M. Hermans, A. Smerieri, M. Haelterman, and S. Massar, “Online Training of an Opto-Electronic Reservoir Computer Applied to Real-Time Channel Equalization,” IEEE Transactions on Neural Networks and Learning Systems 28, 2686–2698 (2017).
  * Harkhoe  _et al._ [2020] K. Harkhoe, G. Verschaffelt, A. Katumba, P. Bienstman, and G. Van Der Sande, “Demonstrating delay-based reservoir computing using a compact photonic integrated chip,” Optics Express 28, 3086 (2020).
  * Sozos, Mesaritakis, and Bogris [2021] K. Sozos, C. Mesaritakis, and A. Bogris, “Reservoir Computing Based on Mutually Injected Phase Modulated Semiconductor Lasers as a Monolithic Integrated Hardware Accelerator,” IEEE Journal of Quantum Electronics 57, 1–7 (2021).
  * Danilenko  _et al._ [2023] G. O. Danilenko, A. V. Kovalev, E. A. Viktorov, A. Locquet, D. S. Citrin, and D. Rontani, “Impact of filtering on photonic time-delay reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 33, 013116 (2023).
  * Donati  _et al._ [2022] G. Donati, C. R. Mirasso, M. Mancinelli, L. Pavesi, and A. Argyris, “Microring resonators with external optical feedback for time delay reservoir computing,” Optics Express 30, 522 (2022).
  * Hülser  _et al._ [2022] T. Hülser, F. Köster, L. Jaurigue, and K. Lüdge, “Role of delay-times in delay-based photonic reservoir computing [Invited],” Optical Materials Express 12, 1214 (2022).
  * Abdalla  _et al._ [2023] M. Abdalla, C. Zrounba, R. Cardoso, P. Jimenez, G. Ren, A. Boes, A. Mitchell, A. Bosio, I. O’Connor, and F. Pavanello, “Minimum complexity integrated photonic architecture for delay-based reservoir computing,” Optics Express 31, 11610 (2023).
  * Larger  _et al._ [2017] L. Larger, A. Baylón-Fuentes, R. Martinenghi, V. S. Udaltsov, Y. K. Chembo, and M. Jacquot, “High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification,” Physical Review X 7, 011015 (2017).
  * Liang  _et al._ [2023] Z. Liang, M. Zhang, C. Shi, and Z. R. Huang, “Real-time respiratory motion prediction using photonic reservoir computing,” Scientific Reports 13, 5718 (2023).
  * Takano  _et al._ [2018] K. Takano, C. Sugano, M. Inubushi, K. Yoshimura, S. Sunada, K. Kanno, and A. Uchida, “Compact reservoir computing with a photonic integrated circuit,” Optics Express 26, 29424 (2018).
  * Guo  _et al._ [2019] X. X. Guo, S. Y. Xiang, Y. H. Zhang, L. Lin, A. J. Wen, and Y. Hao, “Four-channels reservoir computing based on polarization dynamics in mutually coupled VCSELs system,” Optics Express 27, 23293 (2019).
  * Ortín and Pesquera [2017] S. Ortín and L. Pesquera, “Reservoir Computing with an Ensemble of Time-Delay Reservoirs,” Cognitive Computation 9, 327–336 (2017).
  * Hou  _et al._ [2019] Y.-S. Hou, G.-Q. Xia, E. Jayaprasath, D.-Z. Yue, W.-Y. Yang, and Z.-M. Wu, “Prediction and classification performance of reservoir computing system using mutually delay-coupled semiconductor lasers,” Optics Communications 433, 215–220 (2019).
  * Keuninckx, Danckaert, and Van Der Sande [2017] L. Keuninckx, J. Danckaert, and G. Van Der Sande, “Real-time Audio Processing with a Cascade of Discrete-Time Delay Line-Based Reservoir Computers,” Cognitive Computation 9, 315–326 (2017).
  * Yue  _et al._ [2021] D.-Z. Yue, Z.-M. Wu, Y.-S. Hou, C.-X. Hu, Z.-F. Jiang, and G.-Q. Xia, “Reservoir Computing Based on Two Parallel Reservoirs Under Identical Electrical Message Injection,” IEEE Photonics Journal 13, 1–11 (2021).
  * Yang  _et al._ [2022] Y. Yang, P. Zhou, T. Chen, Y. Huang, and N. Li, “Optical neuromorphic computing based on a large-scale laterally coupled laser array,” Optics Communications 521, 128599 (2022).
  * Bogris  _et al._ [2021] A. Bogris, C. Mesaritakis, S. Deligiannidis, and P. Li, “Fabry-Perot Lasers as Enablers for Parallel Reservoir Computing,” IEEE Journal of Selected Topics in Quantum Electronics 27, 1–7 (2021).
  * Paudel  _et al._ [2020] U. Paudel, M. Luengo-Kovac, J. Pilawa, T. J. Shaw, and G. C. Valley, “Classification of time-domain waveforms using a speckle-based optical reservoir computer,” Optics Express 28, 1225 (2020).
  * Ashner  _et al._ [2021] M. N. Ashner, U. Paudel, M. Luengo-Kovac, J. Pilawa, and G. C. Valley, “Photonic reservoir computer using speckle in multimode waveguide ring resonators,” Optics Express 29, 19262 (2021).
  * Antonik, Marsal, and Rontani [2020] P. Antonik, N. Marsal, and D. Rontani, “Large-Scale Spatiotemporal Photonic Reservoir Computer for Image Classification,” IEEE Journal of Selected Topics in Quantum Electronics 26, 1–12 (2020).
  * Rafayelyan  _et al._ [2020] M. Rafayelyan, J. Dong, Y. Tan, F. Krzakala, and S. Gigan, “Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction,” Physical Review X 10, 041037 (2020).
  * Bueno  _et al._ [2018] J. Bueno, S. Maktoobi, L. Froehly, I. Fischer, M. Jacquot, L. Larger, and D. Brunner, “Reinforcement learning in a large-scale photonic recurrent neural network,” Optica 5, 756 (2018).
  * Sunada, Kanno, and Uchida [2020] S. Sunada, K. Kanno, and A. Uchida, “Using multidimensional speckle dynamics for high-speed, large-scale, parallel photonic computing,” Optics Express 28, 30349 (2020).
  * Vandoorne  _et al._ [2014] K. Vandoorne, P. Mechet, T. Van Vaerenbergh, M. Fiers, G. Morthier, D. Verstraeten, B. Schrauwen, J. Dambre, and P. Bienstman, “Experimental demonstration of reservoir computing on a silicon photonics chip,” Nature Communications 5, 3541 (2014).
  * Vinckier  _et al._ [2015] Q. Vinckier, F. Duport, A. Smerieri, K. Vandoorne, P. Bienstman, M. Haelterman, and S. Massar, “High-performance photonic reservoir computer based on a coherently driven passive cavity,” Optica 2, 438 (2015).
  * Nakajima  _et al._ [2018] M. Nakajima, M. Inubushi, T. Goh, and T. Hashimoto, “Coherently Driven Ultrafast Complex-Valued Photonic Reservoir Computing,” in  _Conference on Lasers and Electro-Optics_ (OSA, San Jose, California, 2018) p. SM1C.4.
  * Nakajima, Tanaka, and Hashimoto [2021] M. Nakajima, K. Tanaka, and T. Hashimoto, “Scalable reservoir computing on coherent linear photonic processor,” Communications Physics 4, 20 (2021).
  * Butschek  _et al._ [2022] L. Butschek, A. Akrout, E. Dimitriadou, A. Lupo, M. Haelterman, and S. Massar, “Photonic reservoir computer based on frequency multiplexing,” Optics Letters 47, 782 (2022).
  * Gauthier  _et al._ [2021] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. S. Barbosa, “Next generation reservoir computing,” Nature Communications 12, 5564 (2021).
  * Barbosa and Gauthier [2022] W. A. S. Barbosa and D. J. Gauthier, “Learning spatiotemporal chaos using next-generation reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 093137 (2022).
  * Gauthier, Fischer, and Röhm [2022] D. J. Gauthier, I. Fischer, and A. Röhm, “Learning unseen coexisting attractors,” Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 113107 (2022).
  * Kent, Barbosa, and Gauthier [2024] R. M. Kent, W. A. S. Barbosa, and D. J. Gauthier, “Controlling chaotic maps using next-generation reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 34, 023102 (2024).
  * Redding  _et al._ [2024] B. Redding, J. B. Murray, J. D. Hart, Z. Zhu, S. S. Pang, and R. Sarma, “Fiber optic computing using distributed feedback,” Communications Physics 7, 75 (2024).
  * Lu  _et al._ [2017] Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, and E. Ott, “Reservoir observers: Model-free inference of unmeasured variables in chaotic systems,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 041102 (2017).
  * Zimmermann and Parlitz [2018] R. S. Zimmermann and U. Parlitz, “Observing spatio-temporal dynamics of excitable media using reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 043118 (2018).
  * Cunillera, Soriano, and Fischer [2019] A. Cunillera, M. C. Soriano, and I. Fischer, “Cross-predicting the dynamics of an optically injected single-mode semiconductor laser using reservoir computing,” Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 113113 (2019).
  * Neofotistos  _et al._ [2019] G. Neofotistos, M. Mattheakis, G. D. Barmparis, J. Hizanidis, G. P. Tsironis, and E. Kaxiras, “Machine Learning With Observers Predicts Complex Spatiotemporal Behavior,” Frontiers in Physics 7, 24 (2019).
  * Barmparis  _et al._ [2020] G. Barmparis, G. Neofotistos, M. Mattheakis, J. Hizanidis, G. Tsironis, and E. Kaxiras, “Robust prediction of complex spatiotemporal states through machine learning with sparse sensing,” Physics Letters A 384, 126300 (2020).
  * Doan, Polifke, and Magri [2020] N. A. K. Doan, W. Polifke, and L. Magri, “Learning Hidden States in a Chaotic System: A Physics-Informed Echo State Network Approach,” in  _Computational Science – ICCS 2020_, Vol. 12142, edited by V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos, and J. Teixeira (Springer International Publishing, Cham, 2020) pp. 117–123.
  * Pammi  _et al._ [2023] V. A. Pammi, M. G. Clerc, S. Coulibaly, and S. Barbay, “Extreme Events Prediction from Nonlocal Partial Information in a Spatiotemporally Chaotic Microcavity Laser,” Physical Review Letters 130, 223801 (2023).
  * Packard  _et al._ [1980] N. H. Packard, J. P. Crutchfield, J. D. Farmer, and R. S. Shaw, “Geometry from a Time Series,” Physical Review Letters 45, 712–716 (1980).
  * Rössler [1976] O. Rössler, “An equation for continuous chaos,” Physics Letters A 57, 397–398 (1976).
  * Lorenz [1963] E. N. Lorenz, “Deterministic Nonperiodic Flow,” Journal of the Atmospheric Sciences 20, 130–141 (1963).
  * Kuramoto [1978] Y. Kuramoto, “Diffusion-Induced Chaos in Reaction Systems,” Progress of Theoretical Physics Supplement 64, 346–367 (1978).
  * Sivashinsky [1977] G. Sivashinsky, “Nonlinear analysis of hydrodynamic instability in laminar flames—I. Derivation of basic equations,” Acta Astronautica 4, 1177–1206 (1977).
  * Smerieri  _et al._ [2012] A. Smerieri, F. Duport, Y. Paquot, B. Schrauwen, M. Haelterman, and S. Massar, “Analog readout for optical reservoir computers,” in  _Advances in Neural Information Processing Systems_ , Vol. 25, edited by F. Pereira, C. J. Burges, L. Bottou, and K. Q. Weinberger (Curran Associates, Inc., 2012).
  * Mermelstein  _et al._ [2001] M. D. Mermelstein, R. Posey, G. A. Johnson, and S. T. Vohra, “Rayleigh scattering optical frequency correlation in a single-mode optical fiber,” Optics Letters 26, 58 (2001).
  * Ma  _et al._ [2022] S. Ma, T. M. Antonsen, S. M. Anlage, and E. Ott, “Short-wavelength reverberant wave systems for physical realization of reservoir computing,” Physical Review Research 4, 023167 (2022).
  * Lu, Hunt, and Ott [2018] Z. Lu, B. R. Hunt, and E. Ott, “Attractor reconstruction by machine learning,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 061104 (2018).
  * Griffith, Pomerance, and Gauthier [2019] A. Griffith, A. Pomerance, and D. J. Gauthier, “Forecasting chaotic systems with very low connectivity reservoir computers,” Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 123108 (2019).
  * Wang  _et al._ [2023] X. Wang, B. Redding, N. Karl, C. Long, Z. Zhu, S. Pang, D. Brady, and R. Sarma, “Integrated Photonic Encoder for Terapixel Image Processing,” (2023), arxiv:2306.04554 [physics] .