Anne-Sophie Van Schelt
email: a.vanschelt@amsterdamumc.nl

PhD Lukas M Gottwald

MSc Nienke P M Wassenaar

MD Jurgen H Runge

PhD Ralph Sinkus

MD Jaap Stoker

PhD Aart J Nederveen

PhD Eric M Schrauben

MSc Anne-Sophie Van Schelt

A.-S Van Schelt

Single Breath-Hold MR Elastography for Fast Biomechanical Probing of Pancreatic Stiffness

36% (for becoming eligible for surgery) and has a high toxicity, reducing quality of life. [START_REF] Perri | Response and survival associated with first-line FOLFIRINOX vs gemcitabine and nab-paclitaxel chemotherapy for localized pancreatic ductal adenocarcinoma[END_REF] Typically, PDAC is characterized by prolific stromal deposition, and recently it has been hypothesized that this stromal deposition affects chemotherapy efficacy. [START_REF] Wang | The extracellular matrix: A key accomplice of cancer stem cell migration, metastasis formation, and drug resistance in PDAC[END_REF] As stromal deposition, tumor, and surrounding healthy tissue may all have differing biomechanical properties, the local tissue stiffness could be a potential biomarker to predict stromal deposition and hence chemotherapy efficacy, [START_REF] Erkan | How fibrosis influences imaging and surgical decisions in pancreatic cancer[END_REF] which has been exploited in, for example, diffusion MRI. [START_REF] Nissan | Modifications of pancreatic diffusion MRI by tissue characteristics: what are we weighting for?[END_REF] As such, prior insight in the underlying mechanical structure of a pancreatic tumor could influence treatment plans.

Imaging by MR elastography (MRE) allows for noninvasive determination of tissue visco-elastic properties. [START_REF] Manduca | MR elastography: Principles, guidelines, and terminology[END_REF] Longitudinal (over time) or local (over space) differences in tissue stiffness might add to the growing knowledge of pancreatic disease pathophysiology. Specifically, MRE research to date has mainly focused on large organs, such as the liver and brain. [START_REF] Guo | MR elastography in cancer[END_REF] Assessment of the pancreas with MRE remains challenging due to its small size, elongated shape, and central abdominal location. Shear wave quality is paramount to obtaining reliable results regarding biomechanics, which rely on the analysis and reconstruction of the propagation of shear waves into and through the pancreas. [START_REF] Manduca | MR elastography: Principles, guidelines, and terminology[END_REF] Consequently, pancreatic MRE requires careful MRI sequence parameter optimization. Moreover, the quantitative accuracy of MRE in small structures depends on the local stability of the tissue being imaged over all wave-offsets and motion-encoding directions. [START_REF] Zhuo | MR artifacts, safety, and quality control[END_REF] Multiple abdominal MRE acquisition strategies have been reported, which include three-dimensional (3D) MRE using spin-echo echo planar imaging (SE-EPI) in six breath-holds, [START_REF] Shi | Feasibility of using 3D MR elastography to determine pancreatic stiffness in healthy volunteers[END_REF] multi-slice multi-frequency SE-EPI in free-breathing, [START_REF] Dittmann | Tomoelastography of the abdomen: Tissue mechanical properties of the liver, spleen, kidney, and pancreas from single MR elastography scans at different hydration states[END_REF][START_REF] Tirkes | Magnetic resonance imaging as a non-invasive method for the assessment of pancreatic fibrosis (MINIMAP): A comprehensive study design from the consortium for the study of chronic pancreatitis, diabetes, and pancreatic cancer[END_REF] gradient-echo (GRE) EPI in three consecutive breath-holds, [START_REF] Steinkohl | T1 relaxation times and MR elastography-derived stiffness: New potential imaging biomarkers for the assessment of chronic pancreatitis[END_REF] and GRE sequences in breath-holding. [START_REF] Kolipaka | Magnetic resonance elastography of the pancreas: Measurement reproducibility and relationship with age[END_REF] Free-breathing acquisitions are more comfortable for patients, however they introduce respiratory-dependent position-based errors in spatial sampling of consecutive wave-phase (WP) offsets throughout the acquisition, when sequence repetition time (TR) and respiratory rhythm are not in sync.

The reference MRE acquisition method at our institution is the Ristretto sequence in combination with a gravitational transducer. [START_REF] Guenthner | Ristretto MRE: A generalized multi-shot GRE-MRE sequence[END_REF][START_REF] Runge | A novel magnetic resonance elastography transducer concept based on a rotational eccentric mass: Preliminary experiences with the gravitational transducer[END_REF] The Ristretto sequence is a generalized multi-shot GRE acquisition with fractional motion encoding, encompassing an intrinsic delay for optimal acquisition. [START_REF] Guenthner | Ristretto MRE: A generalized multi-shot GRE-MRE sequence[END_REF] The gravitational transducer relies on the centrifugal force exerted on the transducer casing by an eccentric rotational mass. This method requires four consecutive ($16 seconds) breath-holds, which may introduce respiratory position-based errors in pancreas location, limit spatial and temporal resolution, and can be particularly uncomfortable for patients. [START_REF] Troelstra | Assessment of imaging modalities against liver biopsy in nonalcoholic fatty liver disease: The Amsterdam NAFLD-NASH cohort[END_REF][START_REF] Sinkus | Rheological determinants for simultaneous staging of hepatic fibrosis and inflammation in patients with chronic liver disease[END_REF] Compressed sensing (CS) is an acceleration technique that uses prospective incoherent undersampling of k-space paired with an iterative reconstruction algorithm using a 10/6.9 10/6.9 10/6.9

12/9.2 10/6.9

15/6.9 15/6.9 15/6.9 15/6.9

Flip angle ( In experiment 2, the TR and flip angle were adjusted to accommodate for the 3D acquisition.

MS, multi-slice; 3D, three-dimensional; CS, compressed sensing, with acceleration factor 12; BH, breath-hold; TR, effective repetition time defined as one image-shot; TE, echo time; MEG, motion encoding gradient. Sense = sensitivity encoding acceleration in the phase encoding direction with an acceleration factor 2 or 2.7.

sparsifying transformation and nonlinear optimization. [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF] This allows for faster image acquisition while maintaining similar image quality, and to date CS has been used in a number of MRI techniques in multiple organs of the body. [START_REF] Blanken | Coronary flow assessment using accelerated 4D flow MRI with respiratory motion correction[END_REF][START_REF] De Jonge | Evaluation of compressed sensing MRI for accelerated bowel motility imaging[END_REF][START_REF] Gottwald | Pseudo-spiral sampling and compressed sensing reconstruction provides flexibility of temporal resolution in accelerated aortic 4D flow MRI: A comparison with k-t principal component analysis[END_REF] While not previously investigated, this suggests that CS could be used to accelerate the MRE acquisition to a single breathhold, whilst still providing reliable and high-quality images.

The aim of this study was twofold: 1) to optimize the rapid interleaved (Ristretto) GRE acquisition for pancreatic MRE, and 2) to develop and test a prospectively undersampled and CS-reconstructed single breath-hold 3D MRE in the pancreas.

Materials and Methods

This prospective study was approved by the local institutional review board and all participants gave written informed consent. A total number of 30 healthy volunteers and five patients diagnosed with PDAC were included. Patients were recruited in the context of a larger project on the assessment of pancreatic stiffness using MRE. Patients were eligible with pathological proven PDAC independent of staging or prior treatment and were >18 years old. Healthy volunteers were excluded when they had a history concerning pancreatic diseases (such as pancreatitis, pancreatic masses, or diabetes mellitus) and were under the age of 18.

All scanning was done at 3.0 T (Ingenia, Philips, Best, the Netherlands) with participants positioned head first supine. All participants fasted hours prior to scanning. A gravitational transducer was used, which uses an eccentric mass to induce shear waves in the abdomen. [START_REF] Runge | A novel magnetic resonance elastography transducer concept based on a rotational eccentric mass: Preliminary experiences with the gravitational transducer[END_REF] The transducer with a 3D-printed curved polyactic-plate (to increase the contact area of the transducer with the participant) was strapped at the approximate level of the pancreas head. A reusable gel pad was placed between the transducer plate and the participant for more comfort and better wave transmission. The axial field of view (FOV) was: 336 Â 192 Â 27 mm and 384 Â 252 Â 27 mm depending on participant size, with (3 mm3 ) isotropic voxels in nine slices. T2-weighted images were acquired for anatomical reference for every participant (voxel size 1.2 Â 1.2 Â 3 mm, 20 slices, TR/echo time [TE] = 378/80 msec).

The aims are accomplished in three consecutive experiments: 1) the application of Ristretto GRE-MRE (hereafter, MS-MRE) in the pancreas is optimized in healthy volunteers; 2) the feasibility of pancreatic CS-MRE is assessed and optimized in healthy volunteers and PDAC patients; and 3) the repeatability of CS-MRE is assessed in healthy volunteers.

Experiment 1: Optimization of Multi Breath-Hold Ristretto GRE-MRE for Pancreas

In this experiment we describe the optimization of fractionally encoded multi breath-hold multi-slice GRE-MRE in the pancreas. This was performed in 10 healthy volunteers (female = 4, mean Consideration for the lower vibration frequency was based on a deeper penetration depth, which potentially is beneficial for application in the pancreas. [START_REF] Steinkohl | Practical and clinical applications of pancreatic magnetic resonance elastography: a systematic review[END_REF] This is a trade-off between penetration depth and longer wavelengths for lower frequencies, which decreases resolution of elastographic details. Furthermore, a higher number of WP offsets can improve accuracy of the reconstruction algorithm, as it increases number of sampling points. Lastly, a longer TE allows for longer motion-encoding gradients, which increases the phase accrual caused by the wave-induced displacement of tissue. This improves the accuracy of shear wave quantification in regions of low wave amplitude due to an increased sensitivity to motion. However, it prolongs breath-hold duration, which should be below a clinically feasible breath-hold time (15-20 seconds [START_REF] Parkes | Defining short and prolonged breath-holds[END_REF] ).

The preferred acquisition strategy was established by determining which acquisition showed the parameters with the highest quality for the pancreas. The highest quality was defined as the method with the highest MRE quality parameters, which are described in postprocessing.

Experiment 2: Implementation of Prospective Undersampling and CS Reconstruction for Single Breath-Hold Ristretto GRE-MRE

Single breath-hold 3D Ristretto GRE-MRE testing was performed on 20 healthy volunteers (female = 11, mean age = 30 AE 8 years). In addition, five PDAC patients were included (male = 4, mean age = 69 AE 5 years). MRE acquisitions were performed using the determined optimized parameters from experiment 1: nine slices with 40 Hz vibration frequency, 5 WP, and TE = 6.9 msec. Three consecutive pancreatic MRE scans were performed: (I) standard multi-slice (MS) in four breath-holds (MS-MRE); (II) 3D four breath-hold acquisition (3D-MRE); and (III) 3D CS-accelerated single breath-hold acquisition (CS-MRE). Scan parameters are reported in Table 1. In patients, placement of the gravitational transducer was guided by the T2-weighted images to be as close as possible to the tumor. Patients underwent scans I and III only, to limit scan time and number of breath-holds.

To incorporate prospective CS undersampling a 3D acquisition is recommended, as it gives more freedom in optimal undersampling pattern design. [START_REF] Ye | Compressed sensing MRI: A review from signal processing perspective[END_REF] Incoherent undersampling was achieved using a variable-density Poisson disc, which was varying in undersampling patterns in the WP offset dimension, repeated for each encoding direction. [START_REF] Uecker | Berkeley Advanced Reconstruction Toolbox[END_REF] Furthermore, a phase-conjugate symmetry method (half-scan) in both ky (67%) and kz (89%) directions was introduced. Acceleration factors were chosen to ensure a clinically feasible breath-hold duration (<18 seconds). For CS acquisitions, reconstruction was performed offline using ReconFrame (version 4.4.4, Gyrotools, Zurich, Switzerland) in MATLAB (version 2021a; The MathWorks Inc., Natick, MA, USA) together with the Berkeley Advanced Reconstruction Toolbox (version 0.5.00, BART; https://mrirecon.github.io/bart/) using a total variation sparsifying transform in time (λ = 0.01, 50 iterations). [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] The sparsifying transform was applied in the WP offset dimension.

Experiment 3: Repeatability of Compressed Sensing

Repeatability of in vivo CS-MRE was tested in 15 of the 20 healthy volunteers who had undergone experiment 2 (female = 10, mean age = 32 AE 8 years). Two consecutive single breath-hold MRE acquisitions were performed without repositioning, to test the within-session repeatability. Prospective CS undersampling and reconstruction were done in the same manner as described earlier.

Postprocessing

Inversion of the data to extract the visco-elastic properties was done using dedicated software (ROOT, as described by Fovargue et al). [START_REF] Fovargue | Robust MR elastography stiffness quantification using a localized divergence free finite element reconstruction[END_REF] Briefly, phase data were unwrapped and smoothed, and a Gaussian filter was applied (σ = 0.75, 3 Â 3 Â 3 pixels) before a 3D Fourier transformation was applied using a 3D Blackman-Harris window. This resulted in a 3D shear wave field. The shear wave displacement was calculated locally per voxel and a finite element method was used to extract the visco-elastic properties, including shear wave speed (SWS), phase angle (φ), elasticity (G d ), and viscosity (G I ). Visco-elastic properties were depicted spatially in an elastogram. Average visco-elastic values were calculated by manually drawing regions of interest (ROIs) on the anatomical magnitude elastography images over the pancreatic head and the tumor for healthy volunteers and patients, respectively. ROIs were drawn by a researcher (AS) with over 4 years of experience in pancreatic MRI. Parameters that were used to compare data quality were the local shear wave amplitude, the nonlinearity, and wave signal-to-noise ratio (SNR). The nonlinearity gives an indication of the deviation of the phase signal from a perfect sinusoidal wave and should be below 50%. [START_REF] Prezzi | Characterization of small renal tumors with magnetic resonance elastography[END_REF] The process of drawing ROIs was done with careful consideration of the nonlinearity. The jqj/jdivj, or wave SNR, is the curl of the displacement vector field over the divergence of the vector field (which is zero up to 10 À6 precision due to the incompressible nature of tissue), with a higher ratio indicating increased quality.

Statistical Analysis

Statistical testing was done using SPSS (version 28; IBM Corp., Armonk, NY, USA). Repeatability was assessed using Bland-Altman analysis [START_REF] Bland | Statistical methods for assessing agreement between two methods of clinical measurement[END_REF] and between-subject and within-subject coefficients of variation (CoVs). Data were tested for normality using the Shapiro-Wilk test. Comparison between healthy volunteers and patients was done using an unpaired t test or Mann-Whitney U test in case of skewed data separately for scans I and III to verify if differences in visco-elastic parameters remain after acceleration. Repeated measures of analysis of variance (ANOVA) were performed with a subsequent multiple comparisons with a Bonferroni correction. Data are expressed as mean AE standard deviation (SD) or median (ranges) as appropriate. An (adjusted) P-value <.05 was considered statistically significant.

Results

Study Population

The clinical population consisted of five patients diagnosed with PDAC, regardless of staging, supplemented by 30 healthy volunteers. Three patients received chemotherapy treatment before scanning. One of the five patients underwent a resection 1 week after inclusion. Baseline characteristics are presented in Table 2.

Experiment 1: Optimization of Multi Breath-Hold Ristretto GRE-MRE for Pancreas

Results of experiment 1 are shown in Fig. 1. The mean shear wave amplitude in the pancreas head was observed to be 65 AE 21, 75 AE 29, 90 AE 28, and 82 AE 27 μm for combinations A, B, C, and D, respectively. Statistical analysis exhibited significant differences between the four methods (F [3, 27] = 11). Combination C displayed a significant increase in shear wave amplitude when compared to combination A and B. While A-B (P = .560), A-D (P = .076), B-D (P = .118), and C-D (P = 1.000) remained constant. Nonlinearity showed no significant differences between all combinations (F[3, 27] = 1.5, P = .232). The average shear wave speed was 1.17 AE 0.23, 0.97 AE 0.2, 0.98 AE 0.2, and 0.96 AE 0.19 m/second for combinations A, B, C, and D, respectively. Multiple comparisons demonstrated a significant difference in SWS between combination A-B, A-C, and A-D.

The preferred approach was identified as combination C (40 Hz vibration frequency, 5 WP offsets, and a TE of 6.9 msec) due to the significantly higher shear wave amplitude, absence of variations in shear wave speed and 3 and Fig. 4.

The results showed significant differences in all viscoelastic parameters between healthy volunteers and PDAC patients (Fig. 4). However, phase angle and viscosity differences were not evident in scan I (φ: P = .083 and G l : P = .060). The ratio for SWS of healthy pancreatic tissue in healthy volunteers and tumor tissue in patients with PDAC were 0.71 and 0.72 for scans I and III, respectively.

Figure 5 summarizes comparisons of MRE quality parameters between MS (I) and 3D MRE (II) with 4 breathholds and the single breath-hold CS MRE (III). A significant decrease was found in nonlinearity in healthy volunteers between scans I and II. Nonetheless, wave SNR and SD within an ROI of the SWS demonstrated no significant differences (F[2, 38] = 2.312, P = .113 and F[2, 38] = 1.435, P = .251, respectively).

Experiment 3: Repeatability of Compressed Sensing

The average visco-elastic parameters for scan 1 and scan 2 both were 0.97 AE 0.07 m/second. Bland-Altman analysis of the average SWS showed limits of agreement (LoA) of [À0.09, 0.10] m/second and a bias of À0.006 m/second (Fig. 6). LoA for the φ were [À0.12, 0.11] rad. For SWS the within-subject CoV was 4.8%.

Discussion

Ristretto GRE-MRE accelerated with CS (CS-MRE) allows a single breath-hold MRE acquisition with comparable MRE parameters as conventional four breath-hold multi-slice Ristretto MRE. In addition, both methods displayed comparable MRE quality parameters, such as nonlinearity and wave SNR. Differences were found in pancreatic SWS and elasticity between healthy volunteers and patients with PDAC. These significant differences also hold for the single breathhold acquisition. We found no significant differences in viscosity nor phase angle in a four breath-hold acquisition between healthy volunteers and patients with PDAC. However, differences were found between healthy volunteers and patients with PDAC for viscosity and phase angle when accelerating to a single breath-hold acquisition.

The increased SWS in PDAC patients could be correlated to the increased stiffness of tumorous tissue. Differences in tumor stiffness have been hypothesized to be related to stromal disposition and have been utilized by surgeons by means of palpation. [START_REF] Narkhede | Diagnosis and management of pancreatic adenocarcinoma in the background of chronic pancreatitis: Core issues[END_REF] Moreover, the phase angle, which is the ratio of elasticity over viscosity, has been associated with the fluidity of soft tissue. Previously published literature that looked at the shear modulus phase angle found a difference between healthy volunteers and PDAC. [START_REF] Zhu | Distinguishing pancreatic cancer and autoimmune pancreatitis with in vivo tomoelastography[END_REF] The increase in viscosity and phase angle found when accelerating up to a single breath hold could be caused by high vascularity and permeable vessels, the existence of necrotic tissue, or a shift in the amount of collagen. [START_REF] Shahryari | Reduction of breathing artifacts in multifrequency magnetic resonance elastography of the abdomen[END_REF] The ratio of SWS between healthy pancreatic tissue of healthy volunteers and tumor tissue in patients remained similar for the multi breath-hold and a single breath-hold acquisition. These findings may indicate that a single breath-hold acquisition can still differentiate between healthy and diseased pancreatic tissue. Furthermore, the SWS for healthy pancreas and PDAC were overall lower when compared to previously published literature performed at a 40-Hz vibration frequency. [START_REF] Shi | Feasibility of using 3D MR elastography to determine pancreatic stiffness in healthy volunteers[END_REF][START_REF] Shi | The use of magnetic resonance elastography in differentiating autoimmune pancreatitis from pancreatic ductal adenocarcinoma: A preliminary study[END_REF] This could be due to differences in inversion algorithms. While a direct inversion is most commonly used, the assumptions such as local homogeneity or incompressibility vary. The differences found in stiffness between healthy volunteers and patients with PDAC were comparable to previously published studies for both the four and single breath-hold MRE acquisitions. [START_REF] Shi | Feasibility of using 3D MR elastography to determine pancreatic stiffness in healthy volunteers[END_REF][START_REF] Zhu | Distinguishing pancreatic cancer and autoimmune pancreatitis with in vivo tomoelastography[END_REF][START_REF] Itoh | Feasibility of magnetic resonance elastography for the pancreas at 3T[END_REF] Imaging with CS has become widely used in a number clinical applications. [START_REF] Smith | Potential of compressed sensing in quantitative MR imaging of cancer[END_REF] The choice for an acceleration factor of 12 was determined through careful balancing of scan time (a clinical feasible single breath-hold time for healthy controls and patients of <18 seconds), the loss of spatial resolution introduced by the undersampling, the 3D volumetric coverage (large enough to cover the entire pancreatic head/tumor), and the use of the optimal scan parameters from experiment 1. While sampling and reconstruction parameters were empirically optimized for these criteria, further research is warranted for single breath-hold CS-MRE in other organs, FOVs, or disease types.

The use of a 3D approach in CS-MRE did not have a negative impact on the acquisition efficacy. The flip angle was optimized for the 3D acquisition, which led to an adjusted flip angle of 10 degrees to accommodate the repeated wholevolume excitation compared to 20 degrees in a multi-slice acquisition. Although the longer TR will decrease signal intensity, this is counteracted by a ffiffi ffi 3 p gain in signal intensity achievable with the 3D approach. Moreover, a 3D acquisition allows for lower motion encoding gradient frequencies at a fixed TE, due to a shorter radio frequency (RF) pulse and slab selective gradients, thereby increases encoding efficiency.

Pancreatic MRE using a gravitational transducer and a Ristretto GRE-MRE sequence with fractional motion encoding is feasible. Wave penetration may be the most important factor in the pancreas and was crucial for successful pancreatic MRE. With decreasing shear wave amplitudes, the shear wave displacement may be lower than the noise level and hence becomes indistinguishable from noise. [START_REF] Papazoglou | Scatter-based magnetic resonance elastography[END_REF] Therefore, a mechanical frequency of 40 Hz was chosen as the optimal frequency for application in the pancreas. Compared with other MRE setups which use multiple transducers, [START_REF] Dittmann | Tomoelastography of the abdomen: Tissue mechanical properties of the liver, spleen, kidney, and pancreas from single MR elastography scans at different hydration states[END_REF][START_REF] Gültekin | Added value of tomoelastography for characterization of pancreatic neuroendocrine tumor aggressiveness based on stiffness[END_REF] the single gravitational transducer here was situated near the pancreatic tumor. Multiple transducers may provide a greater depth of penetration but come at the cost of destructive wave interference. Furthermore, discretization effects can introduce errors if there are not enough pixels sampled per wavelength. It has been reported that the optimal ratio is six to nine voxels per wavelength for direct inversions. [START_REF] Manduca | MR elastography: Principles, guidelines, and terminology[END_REF] The ratio in healthy volunteers was 8 voxels per wavelength for both 40 Hz and 50 Hz, indicating that both methods are within the appropriate range. For tumor tissue this ratio increases up to 10, thus indicating that the wavelength becomes larger than the voxel size and inversion could be highly affected by low SNR.

The results obtained for both the within-subject CoV and LoA showed improved results compared to previously published work. [START_REF] Shi | Feasibility of using 3D MR elastography to determine pancreatic stiffness in healthy volunteers[END_REF][START_REF] Kolipaka | Magnetic resonance elastography of the pancreas: Measurement reproducibility and relationship with age[END_REF][START_REF] An | Test-retest reliability of 3D EPI MR elastography of the pancreas[END_REF][START_REF] Shi | Multi-frequency magnetic resonance elastography of the pancreas: Measurement reproducibility and variance among healthy volunteers[END_REF] Shi et al reported a lower withinsubject CoV of 13.4% at 40 Hz compared to 60 Hz. [START_REF] Shi | Feasibility of using 3D MR elastography to determine pancreatic stiffness in healthy volunteers[END_REF] An et al reported 95% LoA of [À0.51, 0.60] kPa between two scans obtained a day apart. [START_REF] An | Test-retest reliability of 3D EPI MR elastography of the pancreas[END_REF] These improved results support the hypothesis of superior repeatability that would result from the stable positioning of the pancreas within a single breathhold compared to deviations in positioning between multiple breath-holds for the other acquisitions. Nonetheless, in a more recent study by Shi et al, the 95% LoA reported were [À0.072, 0.088] m/second at an approximately 7-day interval using a multi-frequency approach. [START_REF] Shi | Multi-frequency magnetic resonance elastography of the pancreas: Measurement reproducibility and variance among healthy volunteers[END_REF]

Limitations

Only a very small cohort of patients with PDAC was included. To investigate differences in underlying mechanical structures within the patient population a larger cohort is necessary. Furthermore, in this study no validation of measured biomechanical properties through histopathology or ex vivo shear testing was done. Ex vivo validation is warranted to substantiate the findings and ensure accuracy. Moreover, the gravitational transducer allows for high amplitude shear waves with accuracy in inversion due to the precise fundamental vibration frequency. However, this single transducer approach does not easily allow adjustable placement and due to the high intrinsic attenuation of shear waves, caused by gas or scattering in the abdomen, vibration of the entire abdomen is challenging from a single point of shear wave introduction (i.e., a single transducer). Lastly, the nonlinearity was slightly higher in PDAC patients compared to healthy volunteers for the multi and single breath-hold acquisitions, stretching towards 50%. However, accuracy of MRE inversion is only guaranteed for a nonlinearity below 50%, thus our findings could be considered borderline acceptable as the ROI shows an average nonlinearity that is just below 50%. [START_REF] Runge | A novel magnetic resonance elastography transducer concept based on a rotational eccentric mass: Preliminary experiences with the gravitational transducer[END_REF] Moreover, for some individual pixels viscosity values were found below zero, which is in violation of conservation of energy and gives an indication for insufficient data quality for assessment of this parameter. This could be caused by the low mechanical frequency and the high acceleration factor that could hamper precision in the inversion, especially in tissues with high stiffness. However, this increase was not substantial compared to nonlinearity and viscosity in healthy volunteers, indicating that this effect is negligible for acceleration up to a single breath-hold. Further research is necessary to ascertain differences within PDAC patients in quantification of stromal disposition and investigate the possibility of single breath-hold MRE as a biomarker for PDAC chemotherapy efficacy.

Conclusion

The use of Ristretto GRE-MRE accelerated with CS may allow for a single breath-hold MRE acquisition with comparable SWS as four breath-hold multi-slice Ristretto MRE. The differences found in SWS and phase angle between healthy volunteers and patients with PDAC were comparable for both the four and single breath-hold MRE.

FIGURE 1 :

 1 FIGURE 1: (Left) T2-weighted image of a representative volunteer with the pancreas delineated in blue on bottom. The subsequent images for the optimal method of the wave amplitude, nonlinearity, and shear wave speed map delineated in red from top to bottom. On the shear wave speed map the pancreas, liver, and kidney are delineated in blue, yellow, and purple, respectively. (Right) From top to bottom: the average wave amplitude, the average nonlinearity, and the average shear wave speed (SWS) in the pancreas head for each method (*P = .01, ***P ≤ .001, ****P ≤ .0001).

FIGURE 2 :

 2 FIGURE 2: Shear wave speed maps of the pancreas overlaid on T2-weighted images of a 73-year-old male patient with pancreatic ductal adenocarcinoma from left to right: scan I (MS 4BH) and scan III (CS-MRE) with the tumor delineated in purple.

FIGURE 3 :

 3 FIGURE 3: Shear wave speed and phase angle maps of the pancreas overlaid on T2-weighted images of a 63-year-old male patient with pancreatic ductal adenocarcinoma from left to right: scan I (MS 4BH) and scan III (CS-MRE) with the tumor delineated in purple.

  nonlinearity, and the capability of clinical implementation with a feasible breath-hold time. Experiment 2: Implementation of Prospective Undersampling and CS Reconstruction for Single Breath-Hold Ristretto GRE-MRE In Figs. 2 and 3, shear wave speed and phase angle maps are shown of a representative patient for each scan. Comparison of visco-elastic parameters between scans revealed no differences in healthy volunteers (F[2, 38] = 0.462, 0.789, 2.620, 0.054, and P = .633, .462, .086, .948 for SWS, φ, G d , and G l , respectively). Similarly, no significant differences were observed in patients between scans I and III (P = .225, .345, .893, .498 for SWS, φ, G d , and G l , respectively). The findings for both healthy volunteers and PDAC patients are summarized in Table

FIGURE 4 :

 4 FIGURE 4: Visco-elastic parameters: (a) shear wave speed (SWS), (b) phase angle (φ), (c) elasticity (G d ), and (d) viscosity (G I ), for all scans. I) MS 4 breath-hold, II) 3D sense accelerated 4 breath-holds, and III) 3D compressed sensing accelerated single breath-hold for the pancreas head in healthy volunteers and tumor in PDAC patients. Mean values for healthy volunteers and PDAC patients are compared separately for scan I and III (*P ≤ .05, **P ≤ .01, ***P ≤ .001).

FIGURE 5 :

 5 FIGURE 5: (a) Healthy volunteers in blue and PDAC tumor SWS in red; (b) nonlinearity for healthy volunteers (HV) and patients (PDAC); (c) Curl over divergence (wave SNR) for HV and PDAC; (d) STDEV within the ROI for healthy volunteers (HV) and patients (PDAC) (*P ≤ .05).

FIGURE 6 :

 6 FIGURE 6: Bland-Altman analysis of the repeatability of the shear wave speed (SWS). Limits-of-agreement (LoA) were À0.09 and 0.10 m/second and the bias was À0.006 m/second.

Table 1 .

 1 Scan Parameters for the Three Experiments: Pancreatic MRE Optimization, Compressed Sensing, and Repeatability

	2. Compressed Sensing
	1. Pancreatic MRE optimization combinations

Table 2 .

 2 Baseline characteristics of the healthy volunteer cohorts and patient cohort

					PDAC Patients
	Characteristics	Healthy Volunteers	1	2	3	4	5
	Age (years)	30.9 AE 9.2	68	63	73	76	67
	Sex (% male)	33%	M	M	M	M	F
	Weight (kg)	70.8 AE 9.3	70	90	82	81	62
	Height (cm)	176 AE 8	170	178	190	185	169
	BMI (kg/m 2 )	22.9 AE 2	2 4	2 8	2 3	2 4	2 2
	Diagnosis at time of MRI	-	2	2	1	3	3
	Resectable PDAC (1)
	Metastatic PDAC (2)
	LAPC (3)
	Tumor location	-	1	3	2	2	1
	Head (1)
	Corpus (2)
	Tail (3)
	Chemotherapy (sessions)	-	>12	>8	0	0	>12
	Time after diagnosis (months)	-	8	11	1	1	2
	Average tumor size (mm)	-	20	6	7	3	4

BMI, body mass index; PDAC, pancreatic ductal adenocarcinoma; MRI, magnetic resonance imaging, M, male, F, female. age = 27 AE 3 years) using four combinations of vibration frequency, number of measured WP offsets, and TE. Table

1

provides a summary of the different combinations tested and motivations for the optimization parameters are given below.

Table 3 .

 3 Mean AE SD or Median (IQR) Values for Healthy Volunteers for Different Quality and Visco-Elastic Parameters of MRE and Mean Values for Visco-Elastic Gd, elasticity, GI, viscosity; MRE, magnetic resonance elastography; SNR, signal-to-noise ratio; STDEV, standard deviation, SWS, shear wave speed; MS 4BH, multi-slice four breath-hold; 3D 4BH, three-dimensional four breath-hold, CS 1BH, compressed sensing single breath-hold. van Schelt et al.: Accelerated Pancreas MR Elastography 15222586, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/jmri.28773 by Cochrane France, Wiley Online Library on [23/11/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	MRE Quality Parameters
	Parameters of PDAC Patients

van Schelt et al.: Accelerated Pancreas MR Elastography 15222586, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/jmri.28773 by Cochrane France, Wiley Online Library on [23/11/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Magnetic Resonance Imaging


Ivan Boikov
email: ivan.boikov@thalesgroup.com

Daniel Brunner

Alfredo De Rossi

Direct coupling of nonlinear integrated cavities for all-optical reservoir computing

We consider theoretically a network of directly coupled optical microcavities to implement a space-multiplexed optical neural network in an integrated nanophotonic circuit. Nonlinear photonic network integrations based on direct coupling ensures a highly dense integration, reducing the chip footprint by several orders of magnitude compared to other implementations. Different nonlinear effects inherent to such microcavities are studied when used for realizing an all-optical autonomous computing substrate, here based on the reservoir computing concept. We provide an in-depth analysis of the impact of basic microcavity parameters on computational metrics of the system, namely, the dimensionality and the consistency. Importantly, we find that differences between frequencies and bandwidths of supermodes formed by the direct coupling is the determining factor of the reservoir's dimensionality and its scalability. The network's dimensionality can be improved with frequency-shifting nonlinear effects such as the Kerr effect, while two-photon absorption has an opposite effect. Finally, we demonstrate in simulation that the proposed reservoir is capable of solving the Mackey-Glass prediction and the optical signal recovery tasks at GHz timescale.

Introduction

Artificial Neural Networks (ANNs) are a computing approach different from the conventional algorithmic concepts. They are inspired by the most basic principles of human brains, where numerous neurons are connected with synapses to create a network. Currently, ANNs are mostly implemented on the von Neumann architecture, i.e. on digital computers, and, while successful in many fields, this approach has its own drawbacks, mainly a data bottleneck inherent to the architecture as well as significant power requirements. In addition, the latency between an input and a computed output presents a major concern in some applications. One example is optical communications, where an exponentially growing demand for internet connectivity requires increasingly faster processing of optical signals [START_REF] Argyris | Photonic neuromorphic technologies in optical communications[END_REF]. A conventional approach is based on the processing with digital algorithms [START_REF] Ip | Compensation of dispersion and nonlinear impairments using digital backpropagation[END_REF][START_REF] Fernando P Guiomar | Mitigation of intra-channel nonlinearities using a frequency-domain Volterra series equalizer[END_REF][START_REF] Mecozzi | Kramers-Kronig coherent receiver[END_REF] or more recent ANN-inspired solutions [START_REF] Häger | Nonlinear interference mitigation via deep neural networks[END_REF]. However, at higher bandwidths, an efficient processing becomes complicated. The development of ANNs in a physical layer can help overcome this challenge [START_REF] Bhavin | Principles of neuromorphic photonics[END_REF][START_REF] Brunner | Photonic reservoir computing[END_REF]. Besides analog electronics [START_REF] Csaba | Coupled oscillators for computing: A review and perspective[END_REF], the manipulation of optical signals is highly relevant when a large bandwidth operation is important, as was demonstrated in the field of microwave photonics [START_REF] Marpaung | Integrated microwave photonics[END_REF]. For implementing such computing systems, integrated photonics is a promising platform that provides compactness and intrinsic nonlinear effects necessary for efficient and scalable ANN in hardware [START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Feldmann | Parallel convolutional processing using an integrated photonic tensor core[END_REF].

Reservoir Computing (RC) has attracted a lot of interest in the in-hardware computing field as it provides an implementation-friendly ANN topology with competitive performance and a simplistic training procedure [START_REF] Verstraeten | An experimental unification of reservoir computing methods[END_REF]. Reservoir Computing on integrated photonic platforms has been demonstrated in simulations [START_REF] Mesaritakis | Micro ring resonators as building blocks for an all-optical high-speed reservoircomputing bit-pattern-recognition system[END_REF][START_REF] Denis | All-optical reservoir computing on a photonic chip using silicon-based ring resonators[END_REF] and experimentally [START_REF] Vandoorne | Experimental demonstration of reservoir computing on a silicon photonics chip[END_REF][START_REF] Sackesyn | Experimental realization of integrated photonic reservoir computing for nonlinear fiber distortion compensation[END_REF]. These works considered optical microcavities or splitters as neurons with long intra-neural connections, i.e. with a non-negligible coupling delay. This was required to match the system's internal transient timescales to available electronics [START_REF] Vandoorne | Experimental demonstration of reservoir computing on a silicon photonics chip[END_REF]. However, these delay lines are costly in terms of a chip footprint and optical losses and limit scalability. As an example, a 280 ps delay requires 2 cm of silicon waveguide [START_REF] Vandoorne | Experimental demonstration of reservoir computing on a silicon photonics chip[END_REF], which results in approximately one neuron per square millimeter of the chip surface.

In this study, we consider evanescent, i.e. direct, coupling between adjacent cavities as an alternative to create maximally compact integrated all-optical computing systems. We propose a RC consisting of a rectangular grid of directly coupled nonlinear microcavities and study the impact of cavity parameters and reservoir geometry on its basic computational properties, namely the dimensionality [START_REF] Dambre | Information processing capacity of dynamical systems[END_REF][START_REF] Skalli | Computational metrics and parameters of an injectionlocked large area semiconductor laser for neural network computing[END_REF] and the consistency [START_REF] Uchida | Consistency of nonlinear system response to complex drive signals[END_REF]. We found that a longer photon lifetime allows for a better scalability due to more variability between coupling-induced supermodes of the system. Importantly, we found that the nonlinear effects can impact computing performance in opposite directions, where two-photon absorption (TPA) reduces the reservoir's dimensionality, while free carrier dispersion (FCD) and the optical Kerr effect increase the dimensionality. Moreover, unlike TPA, a sufficiently strong FCD and Kerr effect do induce chaotic dynamics in this purely dissipative network without gain. Finally, we evaluate the RC's performance in Mackey-Glass prediction and optical signal recovery tasks. For simulations we used cavity and waveguide parameters accepted as standard for GaAs integrated photonic circuits. We found that such integrated photonic reservoirs show a comparable performance to other implementations found in the literature, which shows a great potential of this maximum integration-density approach.

System of coupled nonlinear microcavities

The backbone of our photonic RC is a N × N ⊥ square grid of directly coupled nonlinear cavities (see Figure 1(a)). One side of the grid is coupled to the input waveguide according to cavity-waveguide coupling coefficient κ. Through this waveguide an optical input u(t) is injected. Assuming intra-cavity and cavity-waveguide coupling strengths are weak, the coupled-mode theory can be applied [START_REF] Manolatou | Coupling of modes analysis of resonant channel add-drop filters[END_REF]. For a purely linear case, the electric field of microcavity modes a(t) can be described with

da dt = iω - Γ o + Γ e 2 ⊙ a + ( M µ + M κ )a + K i s(t), (1)
where ω are the cavities' resonant frequencies, Γ o are intrinsic optical losses, Γ e are optical losses to waveguide, ⊙ is the the Hadamard product, M µ , M κ are direct and waveguide-assisted cavity-cavity coupling matrices, K i is an input waveguide-cavity coupling matrix (see A). Here, [START_REF] Haus | Waves and fields in optoelectronics[END_REF] (Section 7.5) if the k-th and m-th cavities are directly coupled and

M µ km = -( M µ mk ) * = |µ km | exp(iϕ km )
K i k = κ i k exp(iϕ i k ), with k, m ∈ 1 . . . N , where N = N N ⊥ . For simplicity, we assume that Γ o k = Γ o , |µ km | = µ and κ i k = κ i , if the k-th cavity is coupled to the waveguide. While this is a good approximation, Γ o k , |µ km | and κ i
k can vary to some degree, e.g. due to fabrication tolerances. The phase terms ϕ km and ϕ i k are defined by a chip geometry and are considered independent and identically distributed (i.i.d.). The strength of intra-cavity and the cavity-waveguide coupling is weak when fulfilling

|iω ⊙ a| ≫ - 1 2 Γ e ⊙ a + ( M µ + M κ )a + K i s(t) . (2)
Direct coupling can only induce an interaction between cavities when resonances sufficiently aligned. For that reason, unless stated otherwise, we assume ω k = ω 0 . Moreover, the direction of cavity modes should be taken into account. For instance, each microring resonator (MRR) resonance is double-degenerate with counter-propagating modes. If a MRR-based reservoir is implemented exactly as in Figure 1(a), both modes of each cavity can be excited and a would have 2N components. In photonic crystal (PhC) cavities, however, resonances are nondegenerate and therefore a has N components. For simplicity we transform the MRR reservoir such that only one mode of each MRR cavity is excited (see Figure 1(b)).

In a fully integrated system, cavities would also be coupled to readout waveguides. This, however, can be challenging to realize for cavities inside the grid, since the compact cavity arrangement leaves no space for such waveguides. One could consider using N × 1 or 2 × N ⊥ grids, where all cavities are accessible. An alternative could be the use of multiple photonic layers [START_REF] Wesley | Monolithically integrated multilayer silicon nitrideon-silicon waveguide platforms for 3-D photonic circuits and devices[END_REF] or three-dimensional waveguides [START_REF] Moughames | Three-dimensional waveguide interconnects for scalable integration of photonic neural networks[END_REF]. Output signals of these waveguides are z(t) = K o ⊙a(t), where K o is a vector of output waveguide coupling coefficients, for which we assume an uniformity according to K o k = κ o . These waveguides introduce an additional optical loss, which has to be included in Γ e .

Input Reservoir Output W in W out W u n x k,n y n κ o N ⟂ N ∥ ¡ a k κ i e s(t) Input waveguide z k (b) (a) (c)
We consider three nonlinear effects commonly present in integrated photonic microresonators: TPA, FCD and the Kerr effect. TPA can be described by adding

- Γ TPA (a k ) 2 a k = - 1 2 β 2 c 2 n 2 V TPA |a k | 2 a k (3)
to the right side of Eq. 1 [START_REF] Uesugi | Investigation of optical nonlinearities in an ultrahigh-Q Si nanocavity in a two-dimensional photonic crystal slab[END_REF]. Here, β 2 is a TPA coefficient, c the speed of light, n the refractive index and V TPA is a nonlinear effective volume [25]

V TPA = ε 2 0 4 ε r χ r 2|U | 4 + |U • U | 2 3 dV -1 , (4)
where

ε r = ε(r)/ε max , χ r = χ (3) (r)/χ (3)
max and U (r) is a cavity's normal mode, normalized such that 0.5

ε 0 ε r |U (r)| 2 dr = 1.
FCD causes a microcavity's resonance to shift by ∆ω k as a function of free electron density N k [START_REF] Bennett | Carrierinduced change in refractive index of InP, GaAs and InGaAsP[END_REF]. This effect adds a term i(dω/dN )N k a k to the right side of Eq. 1 [START_REF] Xu | Carrier-induced optical bistability in silicon ring resonators[END_REF]. Such free electrons could be generated via a linear photon absorption or TPA. The former is possible in direct-bandgap semiconductors, but causes high optical loss independent of the input power, unlike with TPA, where it linearly increases with the optical energy stored in a cavity. Considerable losses in systems without gain make large networks unrealistic. We therefore exclude linear photon absorption as sustainable effect producing nonlinearity and only consider TPA for the free electron generation. The density of TPA-generated electrons is computed by [START_REF] Uesugi | Investigation of optical nonlinearities in an ultrahigh-Q Si nanocavity in a two-dimensional photonic crystal slab[END_REF] dN

k dt = -Γ c N k + Γ TPA (a k ) 2 ωV c I k , (5)
where Γ c is a free electron recombination rate, V c is a volume across which electrons spread through a diffusion, and

I k = |a k | 2
is the optical intensity inside the k-th cavity.

The Kerr effect causes a resonance frequency shift in response to a stronger electric field inside a cavity and is described by adding i(dω/dI)I k a k to the right side of Eq. 1. While the impact of the Kerr effect is typically weaker than FCD, it is still useful to be considered as an isolated case, as, unlike FCD with TPA, it provides a resonance frequency shift without an additional optical losses.

A Photonic Reservoir Computer based on directly coupled nonlinear cavities

A RC is often introduced as an Echo State Network (ESN) [START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF] (see Figure 1(c)):

x n+1 = (1 -α)x n + α tanh( Ŵ x n + Ŵ in [1; u n+1 ]), (6)
where x n and u n are vectors of neuron activations and inputs at a time step n ∈ [1, 2, . . . T ], α is a leaking rate, Ŵ and Ŵ in are a recurrent and an input weight matrices. The ESN's output is defined as y n = Ŵ out x n , where Ŵ out is an output weight matrix. The readout training for software-based models is typically performed using ridge regression [START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF], and here we use the same approach for our simulation. Assume an ESN is being trained to solve a task with input signals u n and the ideal target outputs y tgt n . First, compute update equations for T steps to obtain a set of ESN state vectors x n , which are consequently appended to create a matrix X. Similarly, from y tgt n one creates Ŷ tgt . The optimal output weight matrix is then given by

Ŵ out = Ŷ tgt XT X XT + ζ |X ij | 2 ij Î -1 , (7)
where . . . ij is an average over all matrix indices ij, ζ is a regularization constant and Î is the identity matrix, superscript (. . . ) T is the matrix transpose operation.

In the previous section, we proposed a photonic system that receives optical signal s(t) as an input and provides a set of optical signals z(t). We suggest to treat this system as a reservoir computer. In the integrated photonic hardware, readout weights could be realized based on tunable attenuators with transmission w k , phase-shifters according to ϕ k and a combiner tree to implement a linear combination of z k (t) that creates the output optical signal y(t):

y(t) = k w k exp(iϕ k )z k (t) = W out z(t). (8)
This way, the conversion of optical input s(t) to optical output y(t) is done fully optically, as electronics are only controlling attenuators and phase shifters that are constant during a computation.

Linear regime

We first consider the case when the input signal is not sufficiently strong to induce a nonlinear response by the microcavities. An important value for the following discussion is the input signal's bandwidth BW that we define as the span between its spectral half-power positions. As we will show, computational properties of such photonic reservoirs depend mostly on the microcavities' spectral characteristics. In that regard, an important consequence when directly coupling microcavities is a formation of supermodes with split resonances (see Figure 2). Consider the linear part of Eq. 1:

diag iω - Γ o + Γ e 2 + Mµ + Mκ , (9)
eigenvalues of which are λ E k . Then the average supermode bandwidth is Γ S /2π = -2Re λ E k k /2π. Importantly, the coupling-induced supermode splitting effectively increases the reservoir's response bandwidth significantly beyond the bandwidth of a single resonator. We define the reservoir bandwidth (RBW) as the maximum distance between the network's supermodes measured at half intensity. For example, without the input waveguide, two identical independent cavities with resonances frequencies at ω 0 and coupled with rate µ become split ω 0 ± |µ|, leading to RBW ≈ (2|µ| + Γ S )/2π. For multiple indentical cavities in a chain, Eq. 9 is a tridiagonal Toeplitz matrix, eigenvalues of which are given by [START_REF] Kulkarni | Eigenvalues of tridiagonal pseudo-Toeplitz matrices[END_REF]

λ E k = iω 0 + Γ o 2 + 2 -|µ| 2 cos kπ N + 1 , k = 1 . . . N, (10)
meaning that RBW ≈ (4|µ| + Γ S )/2π. For non-trivial connectivity architectures an analytical derivation of RBW becomes non-trivial. However, in numerical simulations we found that supermode resonances stay within a ω 0 ± π|µ| interval (see Figure 2). We therefore assume that RBW ≈ (2π|µ| + Γ S )/2π for such directly coupled microcavity system. Furthermore, as will be shown later, in the cases we are interested in Γ S ≪ |µ|, and hence RBW ≈ |µ|.

An important computing metric of reservoirs is dimensionality D, which corresponds to the systems' degrees of freedom. Since reservoirs, and ANNs in general, perform a high-dimensional expansion of input signals [START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF], a higher dimensionality usually leads to a better computing capacity. To determine the dimensionality we generate a random sequence s0

(t n ) = √ I 0 (t n ) exp(iϕ 0 (t n ))
, where √ I 0 (t n ) and ϕ 0 (t n ) are randomly sampled according to a white noise distribution. Then, an 8th-order Butterworth low-pass filter with a given bandwidth BW is applied on s0 (t n ) to obtain s 0 (t n ), which is then linearly interpolated to s 0 (t). This signal is used to modulate an optical carrier that is injected into the reservoir as an input. Here and throughout the article the integration of differential equations is carried out using the DifferentialEquations.jl library [START_REF] Rackauckas | Differentialequations.jl-a performant and feature-rich ecosystem for solving differential equations in julia[END_REF] in the Julia programming language [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF]. The photonic reservoir's dimensionality is then determined by applying the Principal Component Analysis (PCA) [START_REF] Pearson | LIII. On lines and planes of closest fit to systems of points in space[END_REF] on z n . Importantly, for the PCA we split electric fields into a real and an imaginary part, which implies that the maximum dimensionality our reservoir can have is 2N . It is therefore convenient to introduce a normalized dimensionality D = D/2N .

At this stage, it is beneficial to consider the supermode concept in more detail. A supermode is an independent, i.e. an orthogonal oscillation of the coupled system's electric fields. Each supermode can be considered as a virtual microcavity, that performs a pass-band filtering at its resonant frequency and bandwidth. Crucially, the coupling-induced splitting of supermode resonances makes their frequencies vary, promising high dimensionality.

To start, a supermode needs to receive an input. For that, we consider the impact of the injection BW on the reservoir's dimensionality. As we see in Figure 3(a), when BW ≪ RBW the dimensionality is low. In this case, a large number of supermodes are off-resonance relative to the injected field, and hence they cannot be excited. With an increasing BW, more and more supermodes are able to interact with the injected information. As a consequence, the reservoir's dimensionality increases until BW = RBW, when all supermodes become excited. Accordingly, the physics of supermodes in coupled arrays result in a necessary bandwidth matching condition for maximizing the system's dimensionality: BW ≥ RBW. Importantly, here the network was driven with white noise, i.e. by a signal with a flat spectrum. However, in realistic and applicationrelevant situations such a signal is unlikely. Certain supermodes will be excited less than others and accordingly contribute less to the networks dimensionality.

The microcavity network's dimensionality relies upon individual supermodes being driven by different aspects of an input signal, which requires supermodes to be sufficiently separated. Therefore, to quantify this metric, we consider an average supermode spacing RBW/N , normalized by the average supermode bandwidth

F Γ = RBW/N Γ S /2π , (11)
which can be considered as a measure of spectral span available for individual supermodes. In Figure 3(b), we see that the dimensionality is low for F Γ ≪ 1 and for F Γ > 1 dimensionality is maximized. This threshold may change depending on system parameters, however, we find that it is typically comparable to unity. Therefore a second requirement for a high dimensionality is F Γ 1. Again, this is not a sufficient condition, since it does not ensure that all supermodes are separated. Spacing between supermodes is typically inhomogeneous, and supermodes can potentially overlap, even if F Γ is large. This effect can be mitigated through a stronger coupling to the input waveguide κ i , which can shift resonances or increase the bandwidth of some supermodes, as can be seen in Figure 2 for a 4 × 3 reservoir with and without the input waveguide. In Figure 3(c), we see that this approach can be effective when N ⊥ is low, i.e. when overall coupling to the input waveguide is relatively stronger for the supermode splitting. From F Γ a scaling behaviour of dimensionality can be determined. Since F Γ ∝ N -1 , a larger number of cavities leads to a lower F Γ , which leads to a lower D. This, however, can be compensated by lowering optical losses, as F Γ ∝ (Γ S ) -1 . Indeed, in Figure 4(a) we see that with a high µ/Γ o (which is analogous to F Γ ), the dimensionality can reach the theoretical limit. Furthermore, dimensionality varies for a given number of cavities depending on their geometric arrangement. A minor fluctuation of dimensionality is due to a seemingly random nature of supermode formation -supermode frequencies of a 24 × 1 grid are different from a 12 × 2 grid. However, there are also outliers with a considerably lower dimensionality. In these cases N ⊥ is high, which allows supermodes to localize spatially further from the input waveguide. As a result, these supermodes become too weakly coupled to the input and as a consequence contribute less to the system's dimensionality (see Figure 4(b)). In addition to that, a lower N ⊥ allows for a higher dimensionality even at a low F Γ , likely due to the impact of input waveguide discussed above.

Nonlinear regime

In order to induce a nonlinear response from the photonic system, an optical power of at least P NL needs to be injected into the cavities. This power depends on material, cavity parameters and which particular nonlinearity is to be leveraged. We define P NL as an input power at which a nonlinear effect becomes equal to a related linear effect. Two-photon absorption increases losses inside the cavities according to Γ TPA , and it is natural to place P NL where these nonlinear losses start matching linear optical losses, i.e. Γ TPA ≈ Γ S . FCD and the Kerr effect cause resonance shift ∆ω, which, in turn, is natural to be compared to an average bandwidth of supermodes, hence also Γ S . Therefore, an optical loss rate is the normalization factor determining the cavity array's nonlinearity power threshold for all considered nonlinear effects.

For a single cavity injected with a monochromatic wave one can determine P NL analytically (see supplementary material). The result is shown in Figure 5(a) for various integrated platforms for MRR and PhC cavities based on a data obtained from literature (see Table 2), and the horizontal axis shows an operating frequency that is BW/2. Assume a cavity with a bandwidth Γ S and an electron relaxation rate Γ c . This cavity can receive most optical power of an injected signal if BW ≤ Γ S /2π. However, a low Γ S corresponds to higher loaded Q-factor and allows for a stronger nonlinearity and, consequently, a lower P NL . We therefore assume that BW ≈ Γ S /2π. Additionally the scale of BW relative to Γ c is also important. In the case when BW ≪ Γ c , the generated free electrons recombine too fast and FCD becomes negligible, whereas TPA becomes the dominant nonlinearity. Then, P NL is computed w.r.t. Γ TPA , see dashed line in Figure 5(a). When BW is comparable to Γ c FCD becomes stronger than TPA. Then, P NL is computed w.r.t. ∆ω and this case is shown as a solid line in Figure 5(a). Finally, when BW ≫ Γ c , the electron density can't react to the changes to the input signal in time and FCD becomes irrelevant, see dotted line in Figure 5(a). Compared to MRRs, PhC cavities require less power and allow a faster operation on the same material, as a smaller mode volume corresponds to a stronger light confinement and a stronger nonlinearity, while a higher surface area enhances Γ c [START_REF] Tanabe | All-optical switches on a silicon chip realized using photonic crystal nanocavities[END_REF].

For a network comprised of a large number of cavities, analytical treatment is not immediately tractable. We therefore resort to a numerical approach. Consider a reservoir with TPA as a nonlinearity. Since its supermodes are split due to direct coupling and, as we found in the previous section, Γ S ≪ RBW, a low-bandwidth excitation of this nonlinearity is not suitable. Therefore, the input changes in time and modulates Γ TPA k on a potentially individual resonator level. In order to then compare it to the scalar Γ S , we consider its standard deviation across many input samples σ t [Γ TPA k ]. Another natural option could be an average over samples, however, consider a case when TPA is dominant and

Γ TPA k (a k (t)) = Γ TPA k t +δΓ TPA k (a k (t)) with δΓ TPA k (a k (t)) ≪ Γ TPA k t ≈ Γ S ,
then the loss terms in Eq. 1 for the k-th cavity is

da k dt = - Γ o k -Γ e k 2 -Γ TPA k t -δΓ TPA k (a k (t)) a k ≈ - Γ o k -Γ e k 2 -Γ TPA k t a k . ( 12
(a k (t)) t ≈ Γ TPA k t ≈ Γ S , whereas with the standard deviation σ t Γ TPA k (a k (t)) = σ t δΓ TPA k (a k (t)) ≪ Γ S . We then compute σ t [Γ TPA k ] k to
obtain a scalar TPA rate for a given input power P . Repeating this process for a range of P we obtain a largely monotonic curve P ( σ t [Γ TPA k ] k ), through which we can determine the input power necessary to induce a particular level of TPAinduced nonlinearity. In that case P NL = P (Γ S ). For FCD and the Kerr effect the process to obtain a scalar characteristic capturing the network's relevant response is identical, except that ∆ω k is used instead of Γ TPA k . In Figure 5(b) we see that P NL scales linearly with the number of cavities, provided N ⊥ is limited. Indeed, for large N ⊥ some supermodes can localize far from the input waveguide and couple to input much weaker than those near it. This creates a nonlinearity imbalance in the reservoir and to reach σ t [Γ TPA k ] k = Γ S , more optical power is needed. Nonlinearity can result in a reservoir operating in a chaotic regime, which for most purposes is harmful to computation. Chaos corresponds to the situation where a network's state is sensitive to infinitesimal changes in the input, and hence the unavoidable presence of noise in hardware would render computations not reproducible for even perfectly identical inputs. For characterizing this effect we carry out a consistency analysis [START_REF] Skalli | Computational metrics and parameters of an injectionlocked large area semiconductor laser for neural network computing[END_REF]. Consistency is a property describing the reproducibility of a dynamical system's (here the reservoir's) responses under an injection of slightly different input signals [START_REF] Uchida | Consistency of nonlinear system response to complex drive signals[END_REF]. We generate a set of inputs in the form s n (t) = s 0 (t) + ∆s n (t) where s 0 (t) and ∆s n (t) are values sampled from a complex-valued white noise distribution, which filtered with a Butterworth 8th order pass-band filter. s 0 (t) is the "true input", while ∆s n (t) emulates an impact of noise in an experiment and we attain a signal to noise ratio of 20 dB by scaling the amplitude of ∆s n (t). Each s n (t) is injected into a separate, but identical reservoir. The absolute value of correlation between trajectories of

γ ij k = E z i k (z j k ) * -E z i k E (z j k ) * E |z i k | 2 E |z j k | 2 , ( 13
)
where E[A] is expected value of A. The consistency for this pair of reservoirs is a root mean square (RMS) across all their cavities γ ij = (γ ij k ) 2 k and equally the consistency of system γ is the RMS of γ ij of all combinations of i and j. It is typically desirable that the consistency is high [START_REF] Lymburn | Consistency in echo-state networks[END_REF].

In Figure 6 we show the consistency and dimensionality of reservoirs with different nonlinearities as a function of input power. Noteworthy, we find that a stronger TPA reduces the dimensionality, while the consistency remains constant and close to unity. The drop of dimensionality is a consequence of the F Γ reduction, as increasing TPA-induced losses effectively increase Γ S . A high consistency could be explained by the negative-feedback nature of TPA, where the system's response to a stronger input is an increased damping that often stabilizes the system. According to the data in Figure 6 FCD and the Kerr effect increase the dimensionality, but eventually the system loses consistency for too high input power. Generally, consistency decreases when σ t [∆ω k ] k approaches Γ S , i.e. when input power surpasses P NL . The effect causing the dimensionality increase is the nonlinear shift of individual supermodes, which otherwise might overlap in the linear regime. We found that generally other system parameters have a relatively small impact compared to P NL .

Interestingly, the nonlinear shift of supermode frequencies also compensates to some degree the negative impact on dimensionality caused by a mismatch between BW and RBW (see Figure 6(d)). This could be explained by the fact that a cavity might belong to multiple supermodes. When a cavity resonance shifts, it causes multiple supermodes to shift as well, even those that are not excited in a linear regime. If the shift is strong enough, these supermodes can end up inside the input bandwidth and receive input. Moreover, as BW reduces, a stronger nonlinearity is needed to move resonances into the input bandwidth, and at some point the concept of supermodes stops being applicable. 1. (d) each point corresponds to an input power at which γ = 0.98.

To conclude, nonlinearity in general has a non-trivial effect on a reservoir. Importantly, this sensitively depends on the particular nonlinear effect that is leveraged. On one hand, dimensionality is reduced by TPA, but increased by FCD and the Kerr effect. On the other hand, using TPA the reservoir stays consistent, while FCD and the Kerr effect are harmful to consistency. However, in the case of FCD and Kerr nonlinearities, the increase of dimensionality comes before the loss of consistency. Hence, there is an interval of optimal input power that allows for a better dimensionality with a high consistency.

Computing performance

We now turn to evaluating the performance when applying a direct coupled cavity array reservoir to computing tests. Here, we consider the Mackey-Glass prediction task, a commonly used benchmark in the RC literature, and the optical signal recovery task that demonstrates a relevant practical application of such a system [START_REF] Argyris | Photonic neuromorphic technologies in optical communications[END_REF].

Mackey-Glass prediction task

The Mackey-Glass equation is a first-order time-delayed differential equation:

dξ dt = αξ(t -τ ) 1 + ξ(t -τ ) g -γξ(t), (14)
where we choose τ = 17, α = 0.2, g = 10 and γ = 0.1, with which ξ(t) exhibits a moderately chaotic behaviour [START_REF] Michael | Oscillation and chaos in physiological control systems[END_REF]. Integrating Eq. 14 with the Euler method with a timestep of 0.17 and downsampling the result with a ratio of 3/0.17, we obtain a discrete series ξ(t n ). The task is to predict y tgt n = ξ(t n+δ ) with ξ(t n ) as an input and a given δ > 0. Both, memory and nonlinearity are required for a successful prediction, which made this task a popular RC benchmark. The prediction performance is measured with the Normalized Root Mean Square Error (NRMSE) [START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF]:

NRMSE y, y tgt = 1 σ n [y tgt n ] 1 T T n=1 y n -y tgt n 2 . ( 15
)
Nevertheless, the reservoir used is time-continuous, we then convert ξ(t n ) into a continuous signal by a zero-hold interpolation. This signal is used to modulate an amplitude of optical carrier that is injected into reservoir as an input. To find the output z k (t) are sampled at t n . Due to the nature of ξ(t n ), the Fourier spectrum of optical input is uneven (see Figure 7). As a result, the excitation of supermodes is not uniform with a non-trivial impact on the reservoir. To visualize the impact of nonlinearities, we consider a linearized Figure 7: The normalized Fourier transform of input signal (blue) with reservoir resonances overlaid (orange). Here BW < RBW, which goes against the conclusions of Section 4. High peaks in the input spectra correspond to a dominating harmonic (ξ(t) largely oscillates, see in Figure 8(c)), however, for a close prediction, it's important to also capture high-frequency irregularities, despite their weak amplitude.

Jacobian of the reservoir

J km = diag iω - Γ o + Γ e 2 km + Mµ km + Mκ km + - 1 2 β 2 c 2 n 2 V TPA I k + i dω dN N k + i dω dI I k δ km , (16)
where last three terms correspond to TPA, FCD and the Kerr effect, with δ km as the Dirac delta. In Figure 8(d) we show trajectories of Ĵ eigenvalues of the system while driven by the input. As expected, TPA increases eigenloss and modifies eigenfrequencies, most likely due to modifications to the spatial shape of the supermodes. With FCD these modifications become notably stronger, and in the case of Kerr effect they are the dominant feature. This behaviour is notably different from ESNs [START_REF] Mustafa C Ozturk | Analysis and design of echo state networks[END_REF]. The inequality of nonlinearity strength is also evident. A supermode near a Fourier peak with a strong connection to an input waveguide is subject to several times stronger nonlinearity than others. This poses a challenge for FCD and the Kerr effect as a nonlinearity concentration in a few supermodes leads the reservoir to chaos before nonlinearity in other supermodes becomes strong enough.

Besides the obvious relevance of the NRMSE for performance, the output signal power also becomes a concern in noisy systems. Even though a noise was not present in simulations, we also consider a power penalty as an additional performance metric: L p = 10 log( |y(t)| 2 t / |s(t)| 2 t ), where . . . t is an average over time.

In Figure 8(a) a considerable improvement with all three nonlinearities is seen. However, with frequency-shifting nonlinearities the performance is lost after the normalized average input power approaches unity. This is consistent with Section 5 as with a loss of consistency the training becomes ineffective. However, in Figure 8(b) we see that before the consistency is lost, the output power is higher compared to TPA, which is likely due to an increase of dimensionality (see Figure 6). For a 8 × 3 reservoir with TPA and FCD a NRMSE of 1.5 • 10 -1 was achieved (see Figure 8(c)), while for TPA 8 • 10 -2 , considering power penalty restriction. For comparison, a performance baseline is NRMSE(u n , u n+3 ) ≈ 1.1, an all-optical time-multiplexed reservoir with 330 virtual neurons based on a semiconductor laser experimentally achived NRMSE of 0.23 [START_REF] Bueno | Conditions for reservoir computing performance using semiconductor lasers with delayed optical feedback[END_REF], a simulated time-multiplexed reservoir based on a nonlinear MRR with 25 virtual nodes had NRMSE of 7.2 • 10 -2 [START_REF] Donati | Microring resonators with external optical feedback for time delay reservoir computing[END_REF]. However, an ESN with 400 neurons outperforms by a large margin, predicting approximately 20 steps ahead with a NRMSE of 1.2 • 10 -4 [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks-with an erratum note[END_REF]. Nevertheless, the proposed reservoir performs such computation fully-optically in real time.

Simulations have shown that prediction performance improves with the number of cavities (see Figure 9) for N ⊥ = 1 or 3 and for all three nonlinearities. With FCD and the Kerr effect N ⊥ = 1 and N ⊥ = 3 reservoirs have shown similar performance, while for TPA N ⊥ = 3 was better than N ⊥ = 1. The lowest error with an NRMSE of 0.037 and L p of -25 dB is achieved for 48 cavities with N ⊥ = 3. We also note the importance of the direct coupling by comparing our system with a set of independent cavities, a linear regime of which was considered by [START_REF] Hermans | Memory in linear recurrent neural networks in continuous time[END_REF]. For a fair comparison, resonances of cavities were set to resonances of supermodes of a corresponding N × 1 reservoir. With all nonlinearities coupled cavities provided considerably better performance, likely due to a nonlinearity-induced supermode mixing.

Figure 8: Impact of nonlinearities on (a) prediction error and (b) power penalty on a 3-step Mackey-Glass prediction. For all cases GaAs MRR parameters are used (see Table 1), but for TPA-only dω/dN = 0, while for Kerr effect β 2 = 0 and dω/dI < 0 is chosen arbitrarily. Regularization constant was increased until power penalty reached -25 dB. (c) 3-step prediction with TPA+FCD case at an optimal input power, here NRMSE is 0.15 and L p ≈ -22 dB. (d) Trajectories of linearized Jacobian J km for (c).

Optical signal recovery

An optical signal propagating through an optical fiber is subject to the chromatic dispersion and the Kerr nonlinearity, among other effects. Electric field E(z, t) in the fiber can be described with a nonlinear Schrödinger equation [START_REF]Nonlinear fiber optics[END_REF]:

i ∂E ∂z + i α loss 2 E - β 2 2 ∂ 2 E ∂t 2 + γ|E| 2 E = 0, (17)
where we choose α loss = 0.2 dB/km, β 2 = 21.7 ps 2 /km, γ = 1.3 W -1 /km, E(0, t) = s 0 (t) and s 0 (t) is the transmitted signal. Because of distortion, after propagating through a fiber of length L the signal at the receiver side E(L, t) = s(t) differs from s 0 (t). Examples of the distortion are shown in Figure 10(a,b). If s 0 (t) is a carrier of encoded symbols, the reservoir's task is to recover them from s(t).

We generate a random sequence s 0 (t n ) of symbols encoded with on-off keying at 25 GHz with 8 samples per bit and filter with a first order low-pass filter with 0.6 • 25 GHz cutoff frequency (see Figure 10(a)). The launch power is 10 dBm. Then Eq. 17 is solved to obtain distorted signal samples s(t n ), which are linearly interpolated to s(t) and injected into the reservoir. The reservoir parameters are the same as in the Mackey-Glass task (see Table 1) except for µ = 40 GHz, chosen according to the s(t) spectrum (see Figure 10(b)) and nonlinear effects disabled. As in the previous section, µ is increased for N × 1 reservoirs to equalize reservoir bandwidth. For RC training we set y tgt (t) = s 0 (tτ ) where τ is an output delay and sample s 0 (t), s(t) and z k (t) once per bit in its middle (see dots in Figure 10(a)). For RC testing y(t) and y tgt (t) are similarly sampled to produce y n and y tgt n and a linear classifier C is used. Performance is measured with the Symbol Error Rate (SER): SER(y,

y tgt ) = n 1 [C(y n ) = C(y tgt n )] /N t ,
where N t is a number of transmitted symbols. However, SER cannot be lower than 1/ min c { n 1[C(y tgt n ) = c]}, which in our case is approximately 2/N t . The presence of an output delay is motivated by group velocity dispersion that causes different input signal frequencies to arrive to a receiver at different times. By delaying y tgt (t) we allow the reservoir to accumulate more information about the current symbol before producing an output. Simulations show that it leads to a better separation of output samples at a given power penalty. For used fiber parameters a delay of 4 bits is adequate and is used throughout the section (see Figure 11).

A larger reservoir is capable of handling a distortion from a longer fiber (see Figure 12). The performance of a N × 1 grid of cavities was comparable to a N × 3 grid and outperformed independent cavities.

Discussion

The proposed setup can successfully be used as a RC. However, there are notable differences that need discussion. First of all, an ESN connection matrix is typically random and sparse [START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF], and two neurons are unlikely to have a connection equally strong in both directions. In the photonic case if two cavities are directly coupled, the connection is reciprocal and coupling coefficients are negative complex conjugate of each other (see Figure 13). Moreover in ESNs a nonlinear function encompasses the input and the neuron connection terms (see Eq. 6). In our reservoir, however, a nonlinearity depends only on the current reservoir state (see Eq. 1). Nonlinearities themselves are differrent as well. Typically, in ESNs a sigmoid or hyperbolic tangent are used. They add a saturation of neuron activation, which is, effectively, an increase of loss. During a typical computation, ESN's linearized Jacobian eigenvalues largely move in the direction of the z-plane origin (corresponding to an increase of loss) and go back with, perhaps, slight frequency deviations [START_REF] Mustafa C Ozturk | Analysis and design of echo state networks[END_REF]. In the photonic case with TPA as the only nonlinearity, the behaviour is similar, but with FCD and the Kerr effect the trajectories of eigenvalues change considerably (see Figure 8(d)). Nonetheless, all three nonlinearities positively contribute to performance (see Figure 8(a)).

In Section 4 a relation between cavity parameters and their effect on reservoir dimensionality was proposed. Whether these requirements could be met with the available technology needs to be discussed. First, consider that Γ o < Γ S , F Γ ≈ 1 and RBW ≤ BW

N Γ o 2π ≤ N Γ S 2π ≈ N F Γ Γ S 2π = RBW ≤ BW ⇒ Γ o 2π N BW. ( 18
)
Another requirement demands of the reservoir to remember previous inputs: an optical timescale should be larger than an input signal timescale, i.e.

1 BW < 2 Γ S < 2 Γ o ⇒ Γ o < 2BW, (19)
which is already satisfied by Eq. 18. For example, consider a reservoir with 20 cavities. To process a signal with BW = 20 GHz, an approximate lower bound for Γ o would be 2π GHz, which would correspond to an intrinsic Q-factor of ω 0 /Γ o ≈ 0.2 • 10 6 . This is rather high, but not impossible. An intrinsic Qfactor of 0.7 • 10 6 was demonstrated in passive GaAs PhCs with 10 6 considered achievable [START_REF] Combrié | GaAs photonic crystal cavity with ultrahigh Q: microwatt nonlinearity at 1.55 µm[END_REF]. In an AlGaAs MRR a Q-factor of 1.5•10 6 was demonstrated [START_REF] Chang | Ultra-efficient frequency comb generation in AlGaAson-insulator microresonators[END_REF] and in a silicon PhC a Q-factor of more than 11 • 10 6 [START_REF] Asano | Photonic crystal nanocavity with a Q factor exceeding eleven million[END_REF]. The low-loss silicon nitride allowed for even higher Q-factor of 37 • 10 6 in a MRR [START_REF] Ji | Ultra-low-loss on-chip resonators with sub-milliwatt parametric oscillation threshold[END_REF].

Next, consider a nonlinearity timescale. TPA and the Kerr effect respond almost instantaneously to a change of electric field, similar to an ESN. However, FCD is also affected by the free electron recombination rate (see Figure 5(a)), which should be comparable to a rate at which of optical signal in a cavity changes, i.e. half of its bandwidth. In our case each cavity belongs to multiple supermodes and thus its bandwidth approaches the RBW, which is close to BW. The electron recombination rate of silicon MRRs was shown to be approximately 2 GHz [START_REF] Stefan F Preble | Ultrafast all-optical modulation on a silicon chip[END_REF], i.e. BW 4 GHz would be supported. The recombination can be accelerated to almost 20 GHz with a reverse-biased p-i-n junction embedded in a cavity [START_REF] Stefan F Preble | Ultrafast all-optical modulation on a silicon chip[END_REF]. Due to a higher surface area, in PhCs the recombination can reach 10 GHz for silicon [START_REF] Sw | Ultrafast band-edge tuning of a two-dimensional silicon photonic crystal via free-carrier injection[END_REF] and 100 GHz for GaAs [START_REF] Combrié | GaAs photonic crystal cavity with ultrahigh Q: microwatt nonlinearity at 1.55 µm[END_REF].

Until this point we assumed that cavities's resonances are aligned ∀k = 1..N ω k = ω 0 . However, this might not be the case due to fabrication tolerances. In order for directly coupled cavities to interact, resonances should be sufficiently close. A weak intercavity interaction leads to two issues, the first being a dimensionality reduction when some cavities are not coupled to the input waveguide directly (e.g. for N ⊥ > 1). Since such cavities can only receive signals through other cavities, without an interaction they stop receiving signals and hence cannot contribute to dimensionality. This is demonstrated in Figure 14 for 8 × 3 and 24 × 1 linear MRR reservoirs. As a resonance disorder is increased, at first D increases as the supermode overlap is alleviated, but decreases afterwards. The 24 × 1 reservoir is affected less than the 8 × 3 one, since in the former 12 cavities are coupled to the waveguide (due to MRR mode direction matching, see Figure 1(b)), but only 4 in the latter. If MRR mode direction was not considered and both MRR modes were excited, all cavities of a 24 × 1 MRR reservoir are coupled to a waveguide. In that case a further increase of resonance disorder simply increases F Γ and D will eventually reach unity. A similar outcome is expected for a PhC reservoir. The second issue is a weakening of supermode mixing, i.e. delocalization. When cavities interact well, supermodes span multiple cavities, otherwise, they isolate in individual cavities. Because of that, a nonlinearity of an individual cavity can only affect a single supermode. As was shown in Figure 8(d), in this case nonlinearity could not contribute to computing performance as much. Here, N ⊥ > 1 cases could be affected to a lesser degree than N ⊥ = 1, since for N ⊥ > 1 each cavity has 2-4 neighbouring cavities, while only 1-2 for N ⊥ = 1 and the probability of resonance mismatch with all neighbours is lower with more neighbours. The weakening of supermode mixing is shown in inset of Figure 14. Each heatmap represents a matrix Λ where Λ km is an intensity of m-th supermode in the k-th cavity. For clarity, each m-th column of Λ is normalized to unity. We see that when σ k [ω k /2π] < 0.5µ, supermodes span multiple cavities, but only a few otherwise. The resonance disorder can be mitigated with the thermo-optic effect, but for more cavities it becomes challenging. In addition to that, a close proximity of cavities results in a strong thermal crosstalk. However, the resonance disorder can be less severe when cavities are spatially close [START_REF] Dodane | Fully embedded photonic crystal cavity with Q= 0.6 million fabricated within a full-process CMOS multiproject wafer[END_REF].

Conclusion

We have applied a network of directly coupled nonlinear nanophotonic cavities for a reservoir computing. Such coupling allows for a lower chip footprint by several orders of magnitude, excluding readout. We have derived general conditions under which such a system attains a high dimensionality, an important property for RC. The role of nonlinear effects present in integrated microcavities, namely TPA, FCD and the Kerr effect were studied. While their impact on the reservoir dimensionality and the consistency is different, all three nonlinearities positively contributed to the performance in Mackey-Glass prediction task. The reservoir was also successful in recovery of OOK-encoded optical signal. In both cases the reservoir has shown a fully optical autonomous real-time computation in the GHz range.

Several potential improvements of this design could be proposed. For example, using multiple cavity resonances to operate on multiple channels at the same time. An another option is to overcome dimensionality limitation due to a finite Q-factor. If a reservoir has an upper limit of N cavities, one could consider using two such reservoirs with one receiving a delayed copy of input. This way the dimensionality of the setup could be expected to increase two-fold.

A Derivation of waveguide-related differential equation matrices

A mode amplitude reduces exponentially outside a resonator, therefore we can assume that each cavity-waveguide interaction happens in a restricted region between two reference planes separated by d k (dashed lines in Figure 15) [START_REF] Manolatou | Coupling of modes analysis of resonant channel add-drop filters[END_REF].

For a single cavity coupled to a waveguide the equations are known [START_REF] Manolatou | Coupling of modes analysis of resonant channel add-drop filters[END_REF]:

Omitting waveguide-induced loss and generalizing for a row of n cavities coupled to one waveguide d

a n dt = M n a n + R n 1 s 1+ + R n 2 s 2+ , (22)
  s 1- s 2-   =   0 T n 1 T n 2 0     s 1+ s 2+   +   Q n 1 Q n 2   a n , (23)
where M n , R n 1 and R n 2 are assumed to be known. Adding a new cavity to the row da dt = κ 1 s 3+ + κ 2 s 4+ ,

         s 3-= e -iβd (s 4+ -κ * 2 a) s 4-= e -iβd (s 3+ -κ * 1 a) s 3+ = e -iβL s 2- s 2+ = e -iβL s 3- , (24)
The goal is to find the matrix representation of a system with n + 1 cavities based on a system with n cavities. We first eliminate the flux between the cavities s 2+ = e -iβL e -iβd (s 4+κ * 2 a) s 3+ = e -iβL T n 2 s 1+ + Q n 2 a n .

Then, cavity dynamics can be rewritten using only new inputs s 1+ and s 4+ :

d an dt = M n a n + R n 1 s 1+ + e -iβL e -iβd R n 2 s 4+ -e -iβL e -iβd R n 2 κ * 2 a, da dt = e -iβL T n 2 κ 1 s 1+ + e -iβL Q n 2 κ 1 a n + κ 2 s 4+ , (27)
Consequently, new outputs are:

s 1-= T n 1 s 2+ + Q n 1 a n = = T n 1 e -iβL e -iβd (s 4+ -κ * 2 a) + Q n 1 a n = = T n 1 e -iβL e -iβd s 4+ + Q n 1 a n -T n 1 e -iβL e -iβd κ * 2 a, (28)

B Comparing evanescent and dissipative coupling

Consider a linear case of two directly coupled lossless cavities. Their state will be described with Eq. [START_REF] Haus | Waves and fields in optoelectronics[END_REF]:

d dt a 1 a 2 = 0 µ -µ * 0 a 1 a 2 , ( 32
)
where µ is a coupling coefficient and * is a complex conjugation. Then a

depending on ϕ 12 and ϕ 21 there could be an interaction ranging from the one similar to direct coupling to its total absence, i.e. light escapes from cavities to the input waveguide. In practice, however, phases are effectively random. As a result, we choose the direct coupling as a neuron connectivity method.

C Estimation of reference power for a single cavity

Assuming a monochromatic excitation with an input power P NL and a steadystate, a power that the cavity receives should be equal to a power lost:

P NL = - d|a| 2 dt = - da dt a * -c.c. = Γ S + β 2 c 2 n 2 V TPA |a| 2 |a| 2 . ( 34
)
For a TPA-dominated case the reference power is found when Γ TPA (a) = Γ S , then

|a| 2 = Γ S n 2 V TPA β 2 c 2 . ( 35
)
For a FCD-dominated case in a steady state (assuming that a detuning does not impact the power absorbed by the cavity):

∆ω = ∂ω ∂N N = ∂ω ∂N β 2 c 2 n 2 V TPA 1 2 ω 0 V c Γ c |a| 4 = Γ S , (36)
|a| 2 = Γ S Γ c n 2 V TPA 2 ω 0 V c Γ c (∂ω/∂N )β 2 c 2 1/2 . (37)

Figure 1 :

 1 Figure 1: (a) Proposed photonic Reservoir Computer (b) MRR mode direction matching and (c) Echo State Network.

Figure 2 :

 2 Figure 2: Supermodes of directly coupled microcavities for various grid geometries. Bottom figure shows a histogram of supermode frequencies.

Figure 3 :

 3 Figure 3: Dimensionality of linear photonic reservoir w.r.t. (a) input signal bandwidth, (b) average supermode spacing and (c) input waveguide coupling strength. Reservoir parameters are given in Table 1 and by default BW > RBW, but in (a) µ increased to 80 GHz to increase upper limit of dimensionality, in (b) each point corresponds to specific µ and Γ o . (c) different grid geometries, see legend.

Figure 4 :

 4 Figure 4: (a) Scalability of linear reservoir dimensionality w.r.t. optical losses. Dimensionalities of all possible combinations of N and N ⊥ for given N are shown as histograms. Reservoir parameters are shown in Table 1, with µ = 80 GHz, reduced k o and BW > RBW. (b) Impact of F Γ and grid geometry on dimensionality. Both figures use the same dataset.

  ) Here, Γ TPA k t effectively increases the linear loss and hence the k-th cavity is close to linear. However, averaging would still classify the cavity as a nonlinear, since Γ TPA

  k

Figure 5 :

 5 Figure 5: (a) Comparison of material platforms w.r.t. P NL and an operating frequency of a single MRR or PhC cavity, details in text. Here, ALD is atomic layer deposition. For references seeTable 2. (b) Scaling of P NL w.r.t. the number of cavities. Reservoir parameters are given in Table 1, BW matches RBW.

Figure 6 :

 6 Figure 6: Effect of (a) TPA (b) TPA and FCD (c) the Kerr effect on reservoir dimensionality and consistency. Reservoir parameters are given in Table 1. (d) each point corresponds to an input power at which γ = 0.98.

Figure 9 :

 9 Figure 9: Scalability of 3-step prediction performance with various geometries and (a) TPA (b) TPA with FCD and (c) the Kerr effect. For each point the regularization constant was increased until power penalty reached -25 dB. Minimum NRMSE (a) 0.037 (b) 0.070 (c) 0.087.

Figure 10 :

 10 Figure 10: (a) Timetraces of real part of original signal (dashed line) and distorted signal after 150 km of optical fiber (solid line). Black dots show bit middle and act as training target, blue and orange are only for visualization and show distorted signal samples with color corresponding to each symbol. (b) Normalized Fourier spectrum of distorted signal after 150 km of fiber (blue) with 8 × 3 reservoir resonances overlaid (orange). (c) Constellation of distorted signal samples after 150 km of fiber, colors correspond to (a).

Figure 11 :

 11 Figure 11: (a) Effect of output delay (shown on legend) on output symbol separability measured with NRMSE at a given power penalty, controlled with regularization parameter. Here L = 150 km and signal is processed with a 8 × 3 linear reservoir. (b)(c) Constellations of output symbols at L p ≈ -20 dB with 0 and 4 bit delay respectively. Figure 10c shows constellation of input signal samples.

Figure 12 :

 12 Figure 12: Performance scaling of optical signal recovery w.r.t. reservoir size. The reservoir is (a) N independent MRRs (b) N × 1 grid (c) N × 3 grid of linear MRRs. For each computation the regularization parameter is increased until the power penalty reached -25 dB.

Figure 13 :

 13 Figure 13: Neuron connectivity graph of (a) a typical ESN and (b) a photonic RC with a grid of cavities. Connection strength is represented with arrow thickness.

Figure 14 :

 14 Figure 14: Impact of cavity resonance frequency disorder on dimensionality. Reservoir parameters given in Table 1, BW > RBW. Each dot corresponds to a N × N ⊥ (shown on legend) linear MRR reservoir with resonances ω k ∈ N (ω 0 , σ ω ) with various σ ω . Insets represent supermode intensity distribution over cavities Λ.

Figure 15 :+ κ 1 s

 151 Figure 15: Scheme of waveguide-cavity interaction. A rectangle with a n represents a row of n cavities coupled to the waveguide.

s 4 - 1 e 2 e -iβL T n 2 κ 1 κ 2  1 e 1 -* 2 e 2 -

 41221122 = e -iβd (e -iβL s 2-κ * 1 a) = = e -iβd e -iβL T n 2 s 1+ + Q n 2 a nκ * 1 a = = e -iβd e -iβL T n 2 s 1+ + e -iβd e -iβL Q n 2 a ne -iβd κ * 1 a.(29)Therefore, the differential equation and the output matrices can be build incre--iβL e -iβd R n e -iβL T n -iβd e -iβL T n e -iβL e -iβd T n 1 κ -iβd e -iβL Q n

Table 1 :

 1 Default reservoir parameters. Most parameters are omitted as they only impact the required input power.

	Parameter	Value	Units	Comments
	Geometry 8 × 3 grid		GaAs MRRs
	µ	25 GHz	defines RBW
	Γ c	28.5 GHz	[33]
	Q-factor	3 • 10 5		provide a near-optimal Γ S /2π ≈ 0.9 GHz ≈ RBW/N
	k 2 i	4π GHz	every other cavity next to input waveguide is coupled
	k 2 o	π GHz	provides sufficiently strong z k (t)

Table 2

 2

. (b) Scaling of P NL w.r.t. the number of cavities. Reservoir parameters are given in Table 1, BW matches RBW.

Table 2 :

 2 References used in Figure5.

	Key	Cavity	Reference
	GaAs
	Ibrahim2002 MRR	[35]
	Combrie2008 L 5 PhC	[36]
	Husko2009	H 0 PhC	[37]
	Moille2015	2× H 0 PhC + ALD [38]
	Silicon
	Xu2005	MRR	[39]
	Tanabe2005	L 3 and L 4 PhC	[34]
	Implanted silicon
	Waldow2008 MRR	[40]
	Tanabe2007	H 1 PhC	[41]
	InP
	Ibrahim2003 MRR	[33]
	Heuck2013	H 0 PhC	[42]

k-th cavity of i-th and j-th reservoirs z i k (t) and z j k (t) is then computed

[START_REF] Kun | Fundamentals of probability and stochastic processes with applications to communications[END_REF]

  1,2 ∝ exp(±i|µ|t), i.e. no loss was introduced regardless of µ. Comparing to a similar case of waveguide coupling of two PhCs in simplified form (see A):

	d dt	a 1 a 2	=	-|κ| 2 |κ| 2 exp(iϕ 21 )	|κ| 2 exp(iϕ 12 ) -|κ| 2	a 1 a 2	,

Acknowledgments

This work was supported by European Union through the Marie Sk lodowska-Curie Innovative Training Networks action, POST-DIGITAL project number 860360.