Thin liquid film as an optical nonlinear-nonlocal medium and memory element in integrated optofluidic reservoir computer

Abstract. Understanding light–matter interaction lies at the core of our ability to harness physical effects and to translate them into new capabilities realized in modern integrated photonics platforms. Here, we present the design and characterization of optofluidic components in an integrated photonics platform and computationally predict a series of physical effects that rely on thermocapillary-driven interaction between waveguide modes and topography changes of optically thin liquid dielectric film. Our results indicate that this coupling introduces substantial self-induced phase change and transmittance change in a single channel waveguide, transmittance through the Bragg grating waveguide, and nonlocal interaction between adjacent waveguides. We then employ the self-induced effects together with the inherent built-in finite relaxation time of the liquid film, to demonstrate that the light-driven deformation can serve as a reservoir computer capable of performing digital and analog tasks, where the gas–liquid interface operates both as a nonlinear actuator and as an optical memory element.


INTRODUCTION
Light-matter interaction is central to advancing our understanding of the properties of both light and matter. Emergence of robust integrated photonic platforms over the last two decades, characterized by miniaturized cross-section and higher refractive index contrast, were leveraged to enhance the optical intensity thus achieving nonlinear response of various matter degrees of freedom. More recently, nonlinear integrated photonics platforms [1] demonstrated promising capabilities to conduct basic scientific research in technologically attractive applications such as frequency conversion and third-harmonic generation [2], supercontinuum generation [3], emerging quantum photonics applications [4], and is considered as an attractive platform for future non-conventional computation architectures [5]. In particular, efficiency of neuro systems to process computational tasks which are challenging to traditional Turing -von Neumann machines [6] due to sequential 'line-by-line' operation, inspired development of machine learning based Neuromorphic Computing (NC) and Recurrent Neural Network (RNN) computational models [7,8] as well as its subset called Reservoir Computing (RC) [9,10]. Similarly to RNN, the recurrent signal in RC constitutes a memory capable to participate in the parallel computation; however, in contrast to RNN, where computationally complex algorithms are required to tune internal weights in the network, in RC the dynamics occurs in a fixed recurrent network and only external weights are digitally tuned thus leading to a significant reduction of the training computation time and the size of the required memory. Consequently, various physical systems can in principle serve as powerful RC platforms [11] where the complex and nonlinear signal output produced by the physical system is collected and used for supervised learning in order to obtain the desired set of weight coefficients during the digital learning stage. Key properties physical RC systems should satisfy include sufficiently large reservoir dimensionality allowing to separate features of distinct states according to their dynamics, short term (fading) memory, and a balance between sensitivity to separability. Among the different proposed physical mechanisms, optical-based systems [12,13] are particularly attractive due to their inherent parallelism, speed of light data propagation and processing, and relatively low energy consumption leading to both on-chip [14][15][16] and free space [17,18] realizations.
In this work we employ the recently theoretically proposed light-liquid nonlinear and nonlocal interaction mechanism [19,20], where localized light-induced heating invokes thermocapillary (TC) [21][22][23] driven deformation of optically thin gas-liquid interface [24,25], to predict and characterize a set of new nonlinear-nonlocal effects in integrated components, and then leverage these effects to achieve RC capabilities. Fig.1 presents schematic illustration of the typical setup where the optical mode which propagates in the waveguide (WG) covered with thin liquid film, invokes heating in a small gold patch on the top facet of the WG due to optical absorption, leading in turn to TC effect and to liquid deformation. In case the liquid film is sufficiently thin, liquid deformation modifies the overlap of the evanescent tail with the gas phase above the liquid film, leading to self-induced change of optical mode properties such as phase, intensity or coupling coefficient to another WG, which we consider below. For numerical simulations we define silicon (Si) or silicon-nitride (SiN) channel WGs designed to carry single mode at wavelength 1550 nm (we consider TM polarization), which is buried in silica (SiO 2 ) layer in order to provide a flat surface for more straightforward numerical implementation. Our extensive 3D multiphysics [26] numerical analysis which takes into account the full dynamics of various degrees of freedom, allows quantitative analysis beyond previous studies which considered limiting assumptions such as small liquid deformation and effective 2D approximation schemes of surface plasmon polariton [19] and slab WG modes [20]. In particular, our simulations include intricate coupling where optical dissipation in the gold patch leads to Joule heat generation and its transport to the gas-liquid interface; the latter triggers surface tension gradients along the gas-liquid interface which in turn invoke TC flow accompanied by liquid's thickness change above the active WG (Fig.1b) (in our terminology active/passive WG corresponds to presence/absence of gold patch). For the simplest case of a single active WG the corresponding deformation of the dielectric cladding leads to a self-induced change of the accumulated optical phase thus constituting a two-way interaction where the optical mode affects itself through the liquid, provided the gas-liquid interface overlaps with the evanescent tail. Furthermore, substantial in-plane spatial scale of liquid's indentation allows to achieve the so-called optical nonlocality where the active optical mode affects the passive optical mode even if the passive WG is displaced from the region of maximal optical intensity. We then employ the self-induced nonlinear phase change in a single FIG. 1. Schematic description of the key components in the integrated optofluidic system under study describing the underlying mechanism of the light-liquid interaction and the memory used for RC. (a) 3D perspective presenting box-shaped liquid cell of length L = 10 µm bracketed by vertical silicone walls, and an active Si channel WG on silica substrate covered with a gold patch of dimensions x × y × k enabling light-induced heating, subsequent heat transport to the gas-liquid interface which in turn triggers surface tension gradients, represented by blue arrows in (b), leading to TC effect and to liquid film thinning. The latter leads to a self-induced phase change and/or to transmittance change in the single WG setup, and also to a nonlocal effect where the phase in the adjacent passive WG is modified due to liquid deformation. Darker colors on the liquid's surface in (a) correspond to higher temperature whereas the dark red lines denote equal height levels. A sequence of high and low optical power pulses correspond to logic '0' and '1', respectively. The finite relaxation time of the gas-liquid interface allows to employ it as a short term memory where n − 1-th pulse of power P (n − 1) induces h c (n − 1) liquid thickness affecting the subsequent n-th pulse of power P (n). (b) 2D normal section, consisting of Si channel WGs of width w = 500 nm and height u = 220 nm buried in depth b = 80 nm under silica top surface and hosting a gold patch of thickness k = 20 nm and in-plane dimensions x = y = 500 nm; the initial liquid depth is h 0 = 500 nm. WG or alternatively the nonlocal coupling change between adjacent WGs, as a nonlinear actuator which provides internal feedback and also exhibits memory capability due to finite relaxation time of the gas-liquid interface, allowing to perform digital and analog RC tasks.

A. Details of multiphysics simulation
To capture the complex coupling between light propagation, fluid dynamics, heat transport and surface tension effects, we employed finite-element based commercial software COMSOL Multiphysics [26] with Wave Optics, heat transfer and CFD module where changes of thin liquid film geometry are simulated by the moving mesh method. In the second step we extracted the corresponding data describing evolution of h c (t) and performed linear fitting to derive the following reduced with power dependent coefficients α(P ), β(P ) (see SM Fig.S5b). While the first step enabled us to generate sufficiently large dataset in a reasonable time, the second step allowed to reduce computation time of the RC algorithm shortly described below. Finally, we implemented RC scheme [27] by using Matlab (R2019b, Mathworks, Natick, MA) [28] to simulate the optical power at the detectors in MZI circuit, described in Fig.5a, as a function of the exciting signal. In particular, we employed Euler method in order to integrate Eq.1 with time step of 0.1 ms and employed the relation Fig.S5c to derive the corresponding phase shift ∆ϕ T C and the accompanying optical signal at the detectors. To implement RC by employing our system we first obtain training matrix, without tuning internal system properties, and then use this matrix to perform test stage which performs relevant computation. In particular, to accomplish the training stage we first collected the output data in D 1,2 and arranged it in N p w ×2 dimensional matrix where N is the number of bits used for training and p w is the excitation time (which affects the number of rows in the output matrix). For instance, in our XOR simulation for training step we used p w = 25, N = 1000. We then rearrange the output matrix in matrix X = (M, N ), where M = 2p w and each column in the matrix X represents data from time step i (i = 1, ..., N ), and use it to solve the following linear equation for the weight matrix W out where Y is a pre-determined desired output. Here, W out is the M × 1 dimensional vector that is determined by employing Tikhonov (ridge) regression, which is expected to be executed digitally in future experimental realizations, via Matlab's built-in 'ridge' function, which minimizes the following expression with β serves as a regularization parameter and |...| denotes the Euclidean norm. Note that the first term in Eq.3 penalizes large differences between the output vector Y and the desired output Y T , whereas the second term penalizes large weight values which facilitates better performance [27]. The corresponding closed form for W out is then given by where non-zero β multiplies a unity matrix I N of dimension N thus adding nonzero values to the diagonal of matrix X T · X and potentially regularizing it. To implement test stage we feed previously unseen driving sequence, generate new design matrix X test and operate on it with the previously derived W out via To finalize the computation result the output vector Y test is subjected to threshold step with values above (below) 0.5 being rounded to one (zero). The corresponding error rate (Er) is then determined by comparing the computation result after the threshold with the true result, leading to where N test is the total number of bits in the test sequence and γ is the number of errors in the computation. For test stage of the XOR task we employed the values N test = 20 and M = 50, whereas for 'zero'/'one' handwritten digits classification we employed N = 12, 665, N test = 2, 115 and M = 28×28×2 = 1, 568 weights. The corresponding ridge parameter value we used to train all models in our work was 0.01, yet determining its optimal value or region of applicability allowing even better performance would require employing ridge parameter estimation methods [49] which is more data science oriented investigation and thus beyond the scope of this work.

A. Nonlinear self-induced phase change
Consider a single Si channel WG with integrated gold patch, covered with thin liquid film described in Fig.2a Fig.2b presents the gas-liquid deformation as a feedback response (circular arrow) allowing to store information and affect subsequent optical pulses at the same compact spatial region, without the need of dedicated large footprint optical feedback elements, e.g., optical rings [29]. The colormap describes electric field intensity at t = 20 ms and thin liquid-deformation along the x − z plane where x is the propagation direction; note the effect of light reflection from the gold patch leading to distortion of the incoming mode. Fig.2c presents the electric field intensity at two different times along the planes x = −4.5 µm, x = 0 µm and x = 4.5 µm, where x = 0 µm is the center of the gold patch. Specifically, the images at the first row indicate that at t = 0 ms, i.e., prior to liquid film deformation, the incident TM mode experiences reduction in intensity but preserves its TM polarization. However, once the liquid film begins to deform, the second row indicates that after 20 ms the incident mode changes its shape due to superposition with the reflected wave from the gold patch, and the transmitted mode is no longer top-bottom symmetric. Fig.2d,e,f,g present, the self-induced phase change, deformation of the thin liquid film above the gold patch center h c (t) (i.e., at region of minimal thickness just above the center of the gold patch), average temperature of the gold patch T c (t), and the corresponding transmittance as a function of time, respectively, under various optical powers. Naturally, increasingly higher power levels and accompanying temperature gradients, lead to more significant h c (t) and the corresponding self-induced phase change effects. Interestingly, at optical powers 0.1 mW, 0, 07 and 0.05 mW the corresponding liquid film deformation induces also prominent transmittance change leading to a power reduction by approximately a factor of 2/3, which we attribute to impedance mismatch created by the liquid deformation. Notably, owing to the relatively low power levels required to activate TC-driven film thinning as well as the large refractive index contrast across the gas-liquid interface, sufficiently thin liquid film with high overlap of the evanescent optical tail with the gas phase supports few orders of magnitude higher self-induced phase change compared to that of thermo-optical (TO) effect. In particular, Fig.S1 in SM indicates that the system described in Fig.2a admits approximately three orders of magnitude higher self-induced phase change compared to a similar system where instead of liquid film one employs a dielectric solid material which doesn't support TC effect. It is worth mentioning that due to the moving mesh method employed in the multiphysics software [26], the simulation terminates once the gas-liquid interface reaches few elements thick film as measured from the bottom of the liquid cell, leading to phase change ∼ 0.3 rad. In future experimental realizations, which are not subject to such limitation, we expect to obtain even higher values of self-induced phase change and nonlocality scales, as well as self-induced transmittance/reflection modulation, shortly described below.
Furthermore, given the fact that liquid deformation is mainly concentrated around the compact gold patch, based on our study the size of the liquid cell can be in principle reduced to 2 µm, whereas higher values of self-induced phase change can be achieved by placing several gold (or other metal) patches along the WG.

B. Nonlocal effect
Consider the optofluidic components schematically presented in  to the close separation regime, initially low power (which is not sufficient to deform the liquid film) optical TM mode was injected also into the passive WG in order to analyze its evolution under the effect of the higher power mode in the active WG. Notably, the maximal separation between the WGs is much higher than the corresponding evanescent scale and also significantly exceeds other known optical nonlocality scales stemming from alternative optical nonlinearities such as photorefractive [31], TO effect [32] and long range molecular interaction between liquid crystals molecules [33][34][35].

C. Self-induced transmittance and reflection in channel Bragg waveguides
Consider SiN Bragg WG on silica substrate presented in Fig.4 where some of the ribs are covered with gold patch of thickness k. Here, we employ SiN rather than Si because the former has a lower refractive index (2 at wavelength 1550 nm) leading to larger mode volume and thus to higher sensitivity of perturbations in the geometry of the gas-liquid interface. First, we design the periodic WG to admit a stop-band at wavelength 1550 nm when it is covered with a thin liquid film of thickness 1 µm (as periodic structure with slightly different period can be designed to support propagation once the WG is covered with liquid film and stop band once the liquid film is replaced by air. Following the same arguments as above, the self-induced liquid thinning over the periodic structure is expected to shift the propagation conditions and lead to lower transmittance as indeed described by the blue curve in Fig.4d. Interestingly, the presented effects are reminiscent of electromagnetically induced transparency which relies on a very different phenomenon where a quantum interference allows the propagation of light through an otherwise opaque atomic medium [36].

D. Reservoir computing of digital and analog tasks
Consider Fig.5a describing Mach-Zehnder interferometer (MZI) circuit with two directional couplers and with an optofluidic cell in one of the arms, which supports the nonlinear self-induced phase change effect described above, and the circuit described by actuation of the liquid, the corresponding actuation time scale (τ ) must be smaller than pulse width (τ w ), whereas in order to ensure memory, i.e., that the liquid deformation imprinted by the n − 1-th bit will affect the phase/coupling change of the n-th bit, the distance between subsequent pulses (τ r ) must be smaller than τ . Both conditions can be summarized as and these ensure that the optically imprinted light pulse in gas-liquid interface at time moment n−1 affects through the liquid the subsequent pulse at time moment n. Since for an optical signal of given power and operation time, thinner film is expected to introduce more significant nonlinear response (unless it is so thin that the optical signal leads to saturation due to very fast liquid depletion above the gold patch), we choose liquid thickness in the region between 0 − 0.5 µm. Below we employ the nonlinear self-induced phase change and the self-induced coupling change (nonlinear effect) between adjacent WGs and the accompanying memory effect, to perform both digital and analog tasks.
Digital task: Consider the delayed XOR operation, where an arbitrary input time series {x n } N n=1 (n = 1, 2, · · · , N ) of zeros and ones is mapped to output sequence {y n } N n=1+δ according to where ⊕ is the addition modulo two (i.e. XOR operation), δ > 0 is an integer encoding the delay. In particular, δ = 1 corresponds to the simplest case of smallest delay where XOR operation performed on adjacent bits whereas δ = 2 corresponds to the case where it is performed on bits which are separated by one bit. The corresponding dynamics of the thin liquid film, applicable irrespective of nonlocal or nonlinear circuit in Fig.5, admits the following evolution equation explicitly given by To implement XOR operation we can either employ the (nonlinear) self-induced phase change or the (nonlinear-nonlocal) self-induced coupling change with circuits Fig.5(a,c), respectively. The input series {x n } N n=1 is encoded as a sequence of square pulses of power level P n = P or 2P depending on the value of logic '0' or '1' as described by Fig.5b; the corresponding power level operates during a time τ w and followed by relaxation time τ r . Fig.5(d,e) present P − τ r performance diagram with different colors encoding the corresponding prediction error, for the case of nonlinear circuit and nonlocal circuit described by Fig.5(a,c), respectively. Note that performance diagram of nonlocal circuit presents (white) region of vanishing error which has larger area compared to the corresponding area in the performance diagram of nonlinear circuit. Furthermore, the lowest power value of the white region in the nonlinear case is ∼ 40 µW and lower value around ∼ 10 µW for nonlocal case. Higher performance of the nonlocal circuit stems from higher sensitivity of the intensity of the transmitted signal with respect to changes of the liquid thickness, compared to the nonlocal case. Specifically, sensitivity estimates can be manifests also for δ = 2 XOR task as can be seen by comparing the corresponding P − τ r performance diagrams of nonlocal circuit with minimal error ∼ 11%, presented in Fig.5f, to performance diagram of nonlinear circuit presented in Fig.S6c with minimal error ∼ 25%. In general, as can be seen from the graph Fig.S6d, presenting test error for a fixed P and τ r , performance of XOR gate quickly degrades for higher values of the delay δ, which also demonstrates prominent memory of only one time backward in time. The underlying dynamics of the thin liquid film was obtained by constructing a reduced 1D model which was obtained from analyzing numerous 3D multiphysics simulations [26] results where the thin liquid film was subjected to various initial conditions and actuation powers. In our reduced model, the transient dynamics of the thin liquid film is subject to first order ordinary differential equation with power dependent coefficients, allowing simple numerical scheme to predict the accompanying phase/coupling change and the resultant power levels at the detectors (see Methods section and SM). In particular,  Fig.5i presents four main curves with internal substructure, corresponding to '1' or '0' which was preceded by '1' or by '0'. These represent the four dominant combinations which allow efficient solution of the delayed (δ = 1) XOR task, whereas internal structure represents next order memory effects which allows to some extent to solve δ = 2 XOR task. More formally, since the dynamics of our system converges to states which do not depend on the initial conditions, it satisfies the so-called common-signal-induced synchronization (CSIS) (see [38] and references within), also known as echo state [8] which is a necessary property enabling a dynamical system to serve as RC. To determine the required optical energy E needed for training, we assume that the number of '0' and '1' pulses is equal leading to E = (N P/2 + N P )τ w ; for instance N = 1000, P = 10 µW and τ w = 20 ms, yields E = 0.1 J. Similarly, the energy during the test stage E test , required for solution of a string of length N test , is given by E test = (N test P/2 + N test P )τ w . For the case N test = 20 and similar parameters used during the training stage discussed above yields E test = 6 µJ.
Analog task: Next, we consider the classification task of hand-written 'zero' and 'one' digits; Fig.6a  Specifically, the pixels are encoded according to power values cP max where P max is the power needed to encode the brightest pixel of unity brightness (c = 1), whereas c = 0 corresponds to dark pixels which carry no features and other values 0 < c < 1 correspond to pixels of intermediate brightness. Fig.6b presents an example of a hand-written 'one' digit used for classification where the 28 ms long signal encoding the highlighted red horizontal row is given by Fig.6c with left/right pixels corresponding to early/late times. To perform image classification we employ either nonlinear or nonlocal circuits presented in Fig.5a and Fig.5c, respectively, where the top and the bottom arms receive in parallel the 1 ms long bits encoding i-th and the following i + 1-th rows, respectively, thus generating more correlations between different lines and leading to higher accuracy compared to feeding just one arm at a time. Fig.6d presents the classification error during the training and testing stages, with training values naturally admitting lower values. Interestingly, both errors admit fairly low values below 0.5 % across a wide range of power values P max spanned between 0.01 mW to 0.1 mW. The total energy needed to feed one image is determined by 28 × 28 × 1 ms × P max ×c, wherec is the average brightness level of the image laying between 0 <c < 1 and typically admitting value ∼ 0.15. The corresponding energy for training/testing stage is achieved by multiplying the one image energy above, by the total number of corresponding images (see Methods section); for P max = 10 µW the corresponding average power level is given bycP max = 1.5 µW. Importantly, the single liquid cell circuits in Fig.5(a,c) can be extended to form more complex network. For instance, Fig.S9 presents a network consisting of four inputs, two liquid cells and two detectors. While it does not present enhancement in accuracy it allows to reduce the computation time by a factor of two without increasing the total required energy. Analyzing networks with more inputs, allowing to achieve optimization of handwritten digits recognition task and further reduction in operation time is beyond the scope of this work.
To learn more about the role of the reservoir in our modality, we also performed linear regression analysis where the output signal is made equal to the input signal which is followed by ridge regression allowing to obtain the corresponding weights that perform the corresponding task with corresponding error given by 0.61%. Interestingly, employing a reservoir presented in Fig.5a but without liquid response, yields similar 0.61% performance for 'zero'/'one' digits classification, which is around 4.3 times higher compared to the highest 0.14% reservoir performance for this task. It is worth mentioning that the performance described in Fig.5d and Fig.S8 is not a prominent function of the optical power and using larger power doesn't lead to significantly enhanced performance. While in this work the input data used for the analog task was not subject to any preprocessing, we expect that common methods such as edge detection should improve the accuracy.

DISCUSSION
To summarize, we believe that the family of physical effects predicted to take place due to nonlinear-nonlocal interaction between thin liquid film to channel WG modes in integrated chip-scale photonics platforms, and the emergent concept suggesting to employ these effects to realize optofluidic RNN supporting RC applications, is expected to stimulate both pure light-matter interaction studies and development of novel RNN architectures/schemes. In particular, our 3D simulation results which take into account optical propagation with losses, heat transport, fluid dynamics as well as the underlying surface tension physics, indicate that the magnitude of the self-induced phase change due to the TC effect in the liquid film, is approximately three orders of magnitude higher compared to the more common heat-based TO effect. Therefore, the predicted effects of self-induced coupling change and self-induced bandstructure and transmittance change can stimulate experimental and theoretical studies validating the predictions and exploring bandstructure tuning in photonic crystals and its topological properties [40,41], as well as studying nonlocal interaction between large number of WGs even if their separation exceeds the optical evanescent scales. While in our work we employed gold patch due to its optical absorption properties, studying metal structures designed to support resonant absorption due to surface plasmon polariton oscillations is another intriguing possibility capable to induce and modify resonant absorption in channel WGs as we briefly mention in Fig.S11. From ML-based perspective, the employed self-induced phase/coupling change and allowing to achieve nonlocal-nonlinear interaction and memory at the same compact actuation region, and capable to perform digital XOR and analog handwritten images classification task, can lead to new RNN elements and architecture which take advantage of the built-in internal feedback capability. Notably, our work demonstrates for the first time that controlled deformation of gas-liquid interface on sub-micron scale, can operate as RNN which can be extended to a network with larger number of liquid cells (or larger number of WGs in a single liquid cell), and further provides design principles of compact optofluidic and complementary metal-oxide-semiconductor (CMOS) compatible optofluidic RC systems which are capable to induce additional nonlinearity beyond the commonly employed square-law detection, without utilizing optical resonant structures, and furthermore operates under wide range of conditions which do not require phase matching conditions. In particular, our design of optofluidic-based RC system is capable to lead to four-five orders of magnitude reduction in size compared to previous liquid-based RC [42,43] systems, and also to few orders of magnitude faster computation compared to the surface waves dynamics and reaction-diffusion processes which were employed in [42] and [43], respectively. While the presented optofluidic approach for RC requires more energy compared to the present state-of-the-art approaches such as opto-electronic based method, offering ∼ 10 GHz modulation rate with hundreds of fJ per bit [44] (see Table S2), we should keep in mind that Electro-optic approach requires external modulation unit whereas in Optofluidic approach the modulation is internal due to the specificity of the gas-liquid interaction. In particular, internal feedback allows to invoke nonlinearities and memory effects in the interaction region and thus translates into very compact computation region without feedback loops which are required in systems that employ external feedback. It is not clear how to achieve self-induced refractive index response of order O(10 −1 ) with Electro-optic effect, though self-induced modulation can be achieved by employing two photon absorption and other nonlinearities [29] by using significantly higher power levels compared to state-of-the-art externally modulated mirco-ring resonator devices [44]. Furthermore, we should keep in mind that neural brain activity also occurs on a ms time scale and is based on several ionic transport processes which also take place in liquid environment. Therefore, even if our study in its current version cannot suggest computational capabilities comparable to that of a human brain, the general strategy of combining several different physical processes rather than relying on a small number of physical effects, may prove to be as a useful strategy to achieve more efficient computational approaches. It is worth mentioning that while phase-change materials, where optically induced phase change between crystalline and amorphous phases invokes self-induced transmittance change, were recently employed for NC applications [45][46][47], the corresponding material phase change is not able to nonlocally affect nearby WGs. Given the fact that in this work we demonstrated that nonlocality enhances performance of both digital and analog tasks, we believe that our study will motivate future studies with more efficient optofluidic-based architectures which employ nonlocality effects, and in parallel also stimulate material science-oriented studies how nonlocality can be harnessed in more conventional phase-change materials. Noteworthy, we performed preparatory experiments demonstrating controlled deposition of non-volatile liquid into etched trenches above silicon channel WG, bringing future experimental realization of the predicted effects within a closer reach (see Fig.S10).

S.1.2. Heat transport
Heat diffusion transport is governed by the following diffusion equation for the temperature field T , Where ρ m , c p , k th are the corresponding material density, specific heat and heat conductivity, and the term J · E is the Joule heat source term.

S.1.3. Fluid dynamics of liquid and gas
Dynamics of a Newtonian, non-compressible fluid in cartesian coordinates is governed by the Navier-Stokes equations given by where τ ij is the corresponding energy-momentum tensor given by Here, u i , F i , ρ, p and µ are the fluid velocity components, body force components, density, pressure, and dynamic viscosity, respectively. The indices i, j run over the three Cartesian coordinates x, y, z, and summation convention over repeated indices is employed.

S.1.4. Gas-liquid matching conditions
The interfacial Stress Balance Equation (SBE) which holds on the gas-liquid (or liquid-liquid or gas-gas) interface, is given by the following matching conditionŝ which in vector notation takes the following form andt, respectively, and are given bŷ where in the last line we usedt · ∇σ ≡ ∇ t σ and assumed the commonly employed linear dependence of surface tension on temperature σ(T ) = σ 0 − σ T T (σ 0 and σ T are typically positive constants).

S.1.5. Boundary conditions
Heat transport: Dirichlet boundary conditions of fixed temperature 20 o on the boundary.
Liquid: Navier boundary conditions on the vertical walls with built-in factor of minimum element length equal to one. On horizontal walls we employed vanishing slip velocity conditions.
Optical field: PML boundary conditions in order to reduce reflections from the boundaries.

S.1.6. Numerical values used
The table below specifies the numerical values employed in the multiphysics simulations. All parameters are at 293.15 K, 1 atm, and for wavelength 1550 nm.

CHANGE
To compare the self-induced phase change due to TC and TO effects, denoted by described in Fig.2 to a similar simulation where the 500 nm thick liquid film is replaced by a solid film of identical thickness and refractive index, but which is not allowed to deform. Fig.S1 presents comparison between the self-induced phase change of the two cases indicating that the magnitude of the effect is of the order ∆ϕ T C /∆ϕ T O ∼ 500 − 1000 Furthermore, while the ∆ϕ T O presents linear behavior in optical power, ∆ϕ T O presents nonlinear behavior due to nonlinear relation and also to saturation due to depletion of the liquid above the gold patch.

S.4. EVOLUTION OF DYNAMICAL PROCESS AS RNN AND PROOF OF EQ.3
First order partial differential equation governing evolution of physical observable A can be written as where H is a function which depends on A, its spatial derivatives ∇A, and time dependent coefficient P (t); J(t) is an additional source term. Introducing spatial discretization and grid points labelled by vectors α, β (e.g. α = (x α , y α , z α ) is a position vector pointing to the relevant grid point in the 3D space), allows to rewrite Eq.S9 as the following system of first order ordinary differential equations where A α (t) ≡ A( r α , t), J α (t) ≡ J( r α , t) and excitation described in Fig.5 the time step coincides with τ W,r which isn't necessarily small), allows to rewrite Eq.S10 as where we have used the Calculus theorem stating that if derivative exists then its left and right limits exist and are equal. Substituting the expression for A α (t) given by Eq.S11 to express A β (t), A γ (t) in the right hand side of Eq.S11, and furthermore keeping first order terms in ∆t yields which admits a functional form of RNN update equation similar to Eq.3. While Eq.S9 describes dynamics of various degrees of freedom such as heat transport, deformation of thin liquid films, mass transport and quantum wave-function, the approach above can be generalized to practically any system including those described by second order systems in time. For the particular case of thin liquid film of constant viscosity, vanishing body and surface forces, as well as vanishing slip velocity on the boundary, Eq.S9 takes the form ∂h ∂t which can be formally rewritten as Eq.S10 (S14) Consequnetly, following the same steps as above one can readily rewrite Eq.S14 as RNN update equation Eq.S12 (i.e. Eq.3 in the main text with different notation for the time argument).

S.5. CONSTRUCTING RC SIMULATION
To construct RC simulation we employed the following three steps: (i) collected

S.5.2. Reduced 1D model for h c (t) and phase change as a function of liquid thickness
In order to incorporate the numerous data curves presented in Fig.S3 into a compact numerical model which allows prediction of h c (t) under various driving optical powers, we fit the data to the following first order ordinary differential equatioṅ Here, h c (t) is liquid thickness above gold patch center as a function of time t, P is a constant value of optical power carried by the WG, and the fit to simulation data is given in Fig.S5a, and α(P ) and β(P ) are power dependent coefficients explicitly given by α(P ) = 0.4368 · P − 0.0589; β(P ) = −0.5546 · P + 0.0469, (S17) and are presented in Fig.S5b. It is worth mentioning that based on a fitting presented in Fig.S5a, a cubic fit would lead to additional terms in Eq.S16 but to small improvement in the model. We then employ Eq.S17 to predict liquid evolution h c (t) under power levels P and 2P (corresponding to logical 'zero' and 'one') during time window τ w , and zero power during τ r . Note that here P stands for the value of the optical power above the WG and not the power in the original input WG; the first coupler divides the power according The one-to-one relation between liquid thickness h c above gold patch center and the self-induced phase change ∆ϕ T C . The fit (black curve) between the two quantities can be interpolated by

S.5.3. Reservoir dynamics: fading memory and delayed 2-bit XOR
Consider evolution equation Eq.S16 with power dependent coefficients α(P ), β(P ) and assume that the power P maintains constant power levels on some time intervals For any interval where the optical power maintains constant value, the solution to Eq.S16 is given by where h 0 is the initial thickness time (t = 0) and α, β are the corresponding values of the coefficients which are constant along the relevant time interval. Employing the solution given by Eq.S18 for each one of the intervals introduced by the partitioning above, yields h n = a n−1 + b n−1 a n−2 + b n−1 b n−2 a n−3 + ...
where only the last term, highlighted by the curly bracket, depends on the initial thickness value h 0 . Here, the constants a k , b k are defined by where α k , β k are the values of α P and β P , respectively, along the time interval T k .
Importantly, the last term h 0 is multiplied by a series of n numbers; each number b k (k = (0, ..., n − 1)) is smaller than unity because α(P ) < 0. Consequently, for sufficiently large number of steps n the last term will become arbitrary close to zero and therefore negligible compared to a n−1 and other terms. Fig.S6a presents the fading memory effect where the evolving thin liquid film, governed by Eq.S16, indeed 'forgets' the initial conditions as we derived in Eq.S19, and all curves approach towards a common attractor state.
Eq.S21 thus allowing simpler implementation of nonlocal RC.

S.5.5. Comparison against optical-based and liquid-based RC platforms
In the following implies that the total heat energy is given by q = ρ · V · c · ∆T ≈ 7.5 nJ. Here, we assumed parameters corresponding to integrated titanium-tungsten (TiW) heaters with density ρ = 4429 kg/m 3 , specific heat c = 526.3 J/(kg · K), and volume V = 0.2 · 4 · 100 = 80 µm 3 of the 100 µm long metal patch. To achieve smaller modulation by a factor of 20, comparable to modulation of Optofluidic circuit presented in Fig.S7 which is sufficient for XOR implementation, requires shorter metal patch of length 100µm/20 = 5 µm and correspondingly smaller modulation energy 7.5/20 ≈ 0.3 nJ. The All-photonic approach on the other hand, presents an example of chip-scale system which enables to obtain a linear response (in electric field) due to delay lines which provide the capability to form interference between signals injected at different times, but it does not have built-in nonlinear matter response apart the square law detection. While in Phase change and Optofluidic modalities the nonlinearity stems from light-matter interaction between WG modes to phase change substance or liquid, in the Fluid-mechanical approach the nonlinearity stems from inherently nonlinear evolution of the gas-liquid interface; linear evolution occurs in a small amplitude regime. Interestingly, while optical nonlinearity in the phase change approach is limited to a single WG, Optofluidic approach is the only one which allows to harness nonlocality for computation and even demonstrates superior performance compared to nonlinear effect (see Fig.5,6).

S.5.6. Analog task: row by row data injection
In Fig.6 we presented performance of nonlinear and nonlocal circuits, schematically presented in Fig.5(a,c), respectively, to perform the analog task of handwritten simultaneously injects a pair of columns instead of rows. Fig.S8 presents simulation results where the nonlinear and nonlocal circuits perform classification of handwritten 'zero'/'one' digits by employing parallel 'bit-after-bit' injection of two pairs of columns (28 bits long each); one column is injected in each one of the two input arms, similarly to the row injection scheme. We see that overall the performance is comparable to rows injection presented in Fig.6, indicating that the degree of correlation generated between adjacent rows/columns is similar in the type of basic circuits we considered. While optimizing the injection scheme and constructing optimized task-dependent circuits is an important question, it is beyond the scope of this work.

S.5.7. Optofluidic RNN
While the circuits presented in Fig.5 employ single or two input WGs, it should require substantial computation time for task which admit large number of bits, including 'zero'/'one' handwritten classification task presented above. To achieve time-division multiplexing one can naturally employ more input waveguides each carrying independent information. Here, we present a basic example towards this direction by employing four inputs which decrease the computation time by a factor of two compared to the schemes presented in Fig.5(a,c). empty, (d) partially filled, (e) fully filled. Note that due to wetting of the vertical walls, liquid thickness is thinner at the center reaching to sub micron thickness which according to our model is enough to initiate self-induced film deformation.
Note that in our modeling (e.g. Fig.2) we assumed initially flat gas-liquid interface, from cell's center to its boundary), and the WG is located at region where the liquid film has a minima point and hence the local slope is even smaller.
Our previous experiments where thin liquid film was deposited on gold plasmonic grating [25], indicates that thin liquid thinning can be invoked by external light illumination across a wide range of power levels. For instance ∼ 400 mW can initiate large deformation which lead to liquid depletion (holes), whereas ∼ 6 µW leads to subnanometer thickness changes. Importantly, in all cases the liquid film restored to initial configuration under surface tension forces. Since our modeling predicts that power levels of just around 0.1 mW (see Fig.2) are enough in order to achieve sufficiently high temperature increase and liquid deformation due to higher power density in the channel WGs, we believe that the proposed approach is promising in terms of future experimental realization.