Dissipation Alters Modes of Information Encoding in Small Quantum Reservoirs near Criticality

Abstract

Quantum reservoir computing (QRC) has emerged as a promising paradigm for harnessing near-term quantum devices to tackle temporal machine learning tasks. Yet, identifying the mechanisms that underlie enhanced performance remains challenging, particularly in many-body open systems where nonlinear interactions and dissipation intertwine in complex ways. Here, we investigate a minimal model of a driven-dissipative quantum reservoir described by two coupled Kerr-nonlinear oscillators, an experimentally realizable platform that features controllable coupling, intrinsic nonlinearity, and tunable photon loss. Using Partial Information Decomposition (PID), we examine how different dynamical regimes encode input drive signals in terms of redundancy (information shared by each oscillator) and synergy (information accessible only through their joint observation). Our key results show that, near a critical point marking a dynamical bifurcation, the system transitions from predominantly redundant to synergistic encoding. We further demonstrate that synergy amplifies short-term responsiveness, thereby enhancing immediate memory retention, whereas strong dissipation leads to more redundant encoding that supports long-term memory retention. These findings elucidate how the interplay of instability and dissipation shapes information processing in small quantum systems, providing a fine-grained, information-theoretic perspective for analyzing and designing QRC platforms.

1. Motivation and Introduction

Reservoir computing (RC) is a computational paradigm that harnesses the intrinsic dynamics of complex systems to process time-dependent inputs efficiently [1,2,3]. Unlike conventional recurrent neural networks (RNNs), RC requires training only at the readout layer, circumventing expensive weight-update procedures on internal nodes [4]. Quantum reservoir computing (QRC) extends these ideas to quantum platforms, leveraging quantum superposition and entanglement to amplify the dimensionality of the feature space and potentially enhance computational capabilities [3,5]. Early demonstrations of QRC have shown promise in tasks like time-series prediction, classification, and memory capacity estimation, and ongoing efforts explore a range of theoretical and experimental strategies for improving performance [6,7,8,9,10,11,12,13,14,15].

Recent QRC research has primarily focused on many-body quantum systems, where quantum phase transitions are suspected to boost computational expressivity [16,17]. While numerical studies reveal intriguing heuristics, such as enhanced memory capacity near critical points, designing optimal quantum reservoirs remains an open question, partly due to the complexity of analyzing large quantum systems. Here, we adopt a complementary approach by studying a pair of coupled Kerr-nonlinear oscillators, a minimal yet experimentally realizable quantum platform [18]. This system exhibits rich dynamical behaviors, including dynamical instability (bifurcation) and dissipation due to photon loss that can be precisely tuned. As such, it offers a tractable yet nontrivial testbed for exploring how quantum correlation, instability, and dissipation govern quantum information processing.

To dissect how these coupled oscillators encode incoming signals, we draw on information-theoretic concepts from neuroscience, where measures of synergy and redundancy helped analyze how neural networks collectively encode stimuli [19,20,21,22]. As the traditional mutual information metric fails to separate out redundant and synergistic contributions to information encoding, we employ Partial Information Decomposition (PID) [23], which partitions the total information into three components: redundancy, capturing information that both oscillators share; unique information, provided by each oscillator individually; and synergy, arising only when both oscillators are observed together. This perspective provides a fine-grained view of the internal encoding structure of the reservoir.

To connect the system’s dynamics to its information-encoding strategy, we combine numerical simulations with non-equilibrium mean-field theory based on the Keldysh formalism [24], focusing on how small external perturbations propagate through the system. In particular, we study how the coupling strength, frequency detuning, and photon loss rate influence the system’s response, and then show how these distinct dynamical regimes lead to different synergy and redundancy profiles in the oscillators’ outputs.

Our main findings reveal that near a critical coupling strength leading to dynamical bifurcation, the system transitions from predominantly redundant encoding to a regime featuring significant synergistic information. This synergistic behavior arises from the interplay between fast collective oscillations and overdamped soft modes. We show that increasing dissipation suppresses quantum correlations and promotes highly redundant encoding modes. In contrast, near the onset of dynamical instability, synergy is amplified and enriches short-term responsiveness, improving short-term memory retention. Taken together, these results highlight how dissipation and dynamic instability in a minimal system can steer a quantum reservoir toward redundant or synergistic processing, each regime benefiting different computational tasks.

In relation to recent proposals using two coupled Kerr-nonlinear oscillators as quantum reservoirs [25,26], which highlight the roles of dissipation for fading memory and moderate coupling for richer dynamics, our work differs by focusing on PID to examine how synergy and redundancy influences the reservoir’s memory capacity. This complements prior findings to show that critical points in Kerr dynamics can shift encoding from redundant to synergistic regimes. Meanwhile, single Kerr oscillators with large Hilbert spaces [13,27] highlight how dimensionality alone can serve as a computational resource, but our key question of whether the whole can exceed the sum of its parts requires at least two coupled oscillators for emergent synergistic encoding. Lastly, although [28] studies larger arrays of Kerr oscillators and demonstrates the near-bifurcation enhancements of nonlinear memory and employ higher-order cumulant expansions to handle higher photon numbers, we limit ourselves to at most a second-order expansion in a parameter regime where it remains accurate, focusing on the PID-based insights rather than reservoir benchmarks and thus complementing other findings in the literature.

To guide the reader, this paper is organized as follows. Section 2 introduces the coupled Kerr-oscillator model and the relevant information-theoretic measures, including PID and quantum mutual information. Section 3 then presents our core numerical findings on synergy and redundancy, comparing fully quantum dynamics with both mean-field and cumulant expansion analyses. We also discuss the mechanisms driving redundant and synergistic encoding and examine how dissipation influences these encoding modes. In Section 3.4, we connect these insights to the quantum reservoir’s memory capacity. Finally, we conclude in Section 4 by summarizing our results and outlining directions for future research. A pedagogical overview of PID can be found in Appendix A and Appendix B, and the details of Keldysh approach to linear response analysis is provided in Appendix C.

2. Quantum Model, Relevant Information Measures, and Performance Metrics

We first introduce our quantum system to study the interplay of instability and dissipation and their influence on modes of information encoding. Then, we outline the key goals of our study, and introduce information measures and metrics to characterize our quantum systems, and finally describe the numerical methods to study them.

2.1. Model

We consider a minimal model that can exhibit dynamical transitions from simple to more complex dynamics, and also support both redundant and synergistic modes of information encoding: a pair of coupled Kerr-nonlinear oscillators. Such systems are well studied in cavity quantum electrodynamics and nonlinear optics, where the interplay of nonlinearities and external driving yields rich dynamical behavior [29,30,31,32]. By focusing on just two coupled cavities, each supporting a single mode at a particular resonance frequency, we avoid the complexity of many-body phases that arise in the thermodynamic limit, thus isolating the essential ingredients needed to explore the onset of coordinated information encoding behaviors in quantum systems. The schematic diagram of this setup is presented in Figure 1.

We work in a frame rotating at the driving frequency ${\omega }_F$, and the corresponding Hamiltonian is given by

(1)$\widehat{H}\left(t\right)=J({\widehat{a}}_1^{†}{\widehat{a}}_2+{\widehat{a}}_2^{†}{\widehat{a}}_1)+\sum _{i=1,2}\left({\Delta }_i{\widehat{a}}_i^{†}{\widehat{a}}_i+{\displaystyle \frac{1}{2}}{U}_i{\widehat{a}}_i^{†2}{\widehat{a}}_i^2+F\left(t\right)({\widehat{a}}_i^{†}+{\widehat{a}}_i)\right),$

where ${\widehat{a}}_i$ and ${\widehat{a}}_i^{†}$ are the annihilation and creation operators for the $i\mathrm{th}$ cavity mode, respectively. The parameter J governs the coupling strength between the two cavities, enabling photon tunneling and collective mode formation [30]. The detuning ${\Delta }_i={\omega }_i-{\omega }_F$ measures the offset of the $i\mathrm{th}$ cavity’s resonance frequency ${\omega }_i$ from the driving frequency. The nonlinear Kerr coefficient $U_i$ characterizes the anharmonicity of each mode which introduces photon-photon interactions essential for generating nonclassical states. And $F\left(t\right)$ is a common external drive that contains time-dependent information applied to both cavities.

To observe transitions between redundant and synergistic information encoding, the two cavities must differ in their nonlinear properties. In particular, having distinct Kerr coefficients ($U_1\ne U_2$) breaks symmetry and enables nontrivial interactions between the modes. With these minimal ingredients (a coherent drive, a tunable coupling, and carefully chosen nonlinearities), this system provides a controlled setting to study the fundamental mechanisms underlying both redundant and synergistic coding in interacting quantum systems.

To incorporate noise and dissipation, we consider the time evolution of an open quantum system weakly coupled to a Markovian bath. Specifically, we model the dynamics of the system’s density matrix $\widehat{\rho }\left(t\right)$ using the Lindblad master equation [33]:

(2)${\displaystyle \frac{d}{dt}}\widehat{\rho }\left(t\right)=\mathcal{L}\widehat{\rho }\left(t\right)=-i[\widehat{H}\left(t\right),\widehat{\rho }\left(t\right)]+\sum _{i=1,2}2{\gamma }_i\mathcal{D}\left[{\widehat{a}}_i\right]\widehat{\rho }\left(t\right),$

where ${\gamma }_i$ is the photon decay rate of the ith cavity associated with the Lindblad superoperator describing single-photon loss, $\mathcal{D}\left[{\widehat{a}}_i\right]$, which acts on the density matrix as

(3)$\mathcal{D}\left[{\widehat{a}}_i\right]\widehat{\rho }={\widehat{a}}_i\widehat{\rho }{\widehat{a}}_i^{†}-\frac{1}{2}\{{\widehat{a}}_i^{†}{\widehat{a}}_i,\widehat{\rho }\}.$

In our simulations, we take ${\gamma }_1={\gamma }_2=\gamma$ for simplicity.

For the common time-dependent external driving field $F\left(t\right)$, we choose $F\left(t\right)=s\left(t\right)F$, where $s\left(t\right)$ is a dimensionless, time-dependent signal, and F is a characteristic strength of the drive. To ensure that the system’s intrinsic dynamics dominate, we select F to be small or comparable to other energy scales, and regard $F\left(t\right)$ as a perturbation. In this work, $s\left(t\right)$ is taken to be a symmetric telegraph process with $s\left(t\right)\in \{-1,1\}$ [34]. While telegraph noise may not be a directly implementable noise model in all bosonic systems (an open quantum system can couple to telegraph noise if it interacts with a fermionic bath; see [35,36] for related studies), we employ it here as a convenient testing ground. Its well-characterized statistical properties [34] and ease of numerical simulation make it a useful drive model for probing how redundant and synergistic information encoding emerges in quantum dynamics. In Section 3.2, we also compare the results with those obtained using different noise models to assess their generality.

2.2. From Quantum to Semiclassical (Mean-Field) Dynamics

Directly simulating the Lindblad equation is practical only when the average photon number is small, as the Hilbert space dimension grows rapidly with photon occupation. In this low-photon regime, we simulate the full quantum dynamics in Equation (2) to capture all quantum correlations, compute PID and QMI, and analyze how nonclassical effects influence information encoding.

As we increase the driving strength or adjust parameters to reach higher photon-number regimes, the full quantum simulation becomes computationally expensive. In this regime, quantum fluctuations often play a less significant role, and a semiclassical approximation becomes suitable. By factorizing expectation values as $⟨{\widehat{a}}_i{\widehat{a}}_j⟩\approx ⟨{\widehat{a}}_i⟩⟨{\widehat{a}}_j⟩$, the dynamics reduce to coupled nonlinear ordinary differential equations (ODEs) for ${\alpha }_i\left(t\right)=⟨{\widehat{a}}_i⟩$:

(4)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\alpha }_1& =-(\gamma +i\Delta ){\alpha }_1-iJ{\alpha }_2-i{U}_1{\alpha }_1{\left|{\alpha }_1\right|}^2-iF\left(t\right),\hfill \\ \hfill {\displaystyle \frac{d}{dt}}{\alpha }_2& =-(\gamma +i\Delta ){\alpha }_2-iJ{\alpha }_1-i{U}_2{\alpha }_2{\left|{\alpha }_2\right|}^2-iF\left(t\right).\hfill \end{array}\end{array}$

These ODEs are much easier to solve, allowing us to examine information encoding under conditions where the photon number is large and quantum correlations are negligible.

To bridge the gap between the fully quantum and semiclassical treatments, we also employ a second-order cumulant expansion (see Appendix D). This approach partially restores some quantum correlations while remaining more tractable than the full density-matrix simulation. We expect that in parameter regimes where quantum correlations matter, the cumulant expansion will improve upon the semiclassical approximation but still remain simpler than the full quantum approach.

2.3. Characterizing Information Processing

To characterize how our system of coupled Kerr-nonlinear oscillators processes and encodes the input signal $s\left(t\right)$ into the output readouts taken to be

(5)$X_i\left(t\right)\equiv \mathrm{Re}⟨{\widehat{a}}_i\left(t\right)⟩,$

we analyze three complementary figures of merit. First, we use the Partial Information Decomposition (PID) to separate the total information that the output observables encode about the input into redundant and synergistic components. Second, we employ the quantum mutual information (QMI) to quantify the role of quantum correlations in shaping these encoding modes. Finally, we consider the memory capacity in a quantum reservoir computing (QRC) context to assess how information is retained over time. While PID and QMI directly characterize the system’s response to external inputs without any training procedure, the memory capacity inherently involves a training step to quantify how well past inputs can be reconstructed from the system’s outputs. In this way, all three measures together provide a comprehensive view of the system’s information processing capabilities.

2.3.1. Partial Information Decomposition (PID)

Let s, $X_1$, and $X_2$ be random variables representing the input signal and the measured observables from the two oscillators, respectively. The mutual information $I(s:X_1,X_2)$ can be decomposed into redundant, synergistic, and unique components as

(6)$I(s:X_1,X_2)=\mathrm{Rdn}+\mathrm{Syn}+\mathrm{Unq}\left(X_1\right)+\mathrm{Unq}\left(X_2\right),$

where $\mathrm{Rdn}$ is the redundant information present in both $X_1$ and $X_2$, $\mathrm{Unq}\left(X_i\right)$ is the unique information contributed solely by $X_i$, and $\mathrm{Syn}$ is the synergistic information accessible only through the joint knowledge of $X_1$ and $X_2$.

By constructing the empirical joint distribution $P(s,X_1,X_2)$ from simulation data and applying the BROJA-2PID algorithm [37], we isolate $\mathrm{Rdn}$ and $\mathrm{Syn}$. This allows us to determine whether the oscillators encode input information redundantly or synergistically, thereby shedding light on their cooperative information processing strategies across different dynamical phases of the system. More detailed discussions and example calculations of PID can be found in Appendix A and Appendix B.

2.3.2. Quantum Mutual Information (QMI)

Although our primary Partial Information Decomposition (PID) analysis focuses on classical correlations between input-output variables, the underlying reservoir dynamics remain fundamentally quantum. To probe the quantum correlations inherent in our system, we compute the quantum mutual information (QMI) between the two oscillators.

Let ${\rho }_{12}$ be the density matrix describing the joint state of the two coupled Kerr oscillators. We partition the system into subsystems 1 and 2, each corresponding to one oscillator. The QMI is then given by

(7)$I(1:2)=S\left({\rho }_1\right)+S\left({\rho }_2\right)-S\left({\rho }_{12}\right),$

where ${\rho }_1={\mathrm{Tr}}_2\left({\rho }_{12}\right)$ and ${\rho }_2={\mathrm{Tr}}_1\left({\rho }_{12}\right)$ are the reduced density matrices of each oscillator, and $S\left(\rho \right)=-\mathrm{Tr}\left(\rho ln\rho \right)$ is the von Neumann entropy.

Because each oscillator’s Hilbert space is, in principle, infinite dimensional, we employ a finite photon-number basis up to a cutoff $N_{\mathrm{cutoff}}$ to ensure computational tractability. Namely,

(8)${\rho }_{12}\approx \sum _{n_1=0}^{N_{\mathrm{cutoff}}}\sum _{n_2=0}^{N_{\mathrm{cutoff}}}{c}_{n_1,n_2}|n_1⟩⟨n_1|\otimes |n_2⟩⟨n_2|,$

where $|n_i⟩$ is the Fock state with $n_i$ photons in oscillator i, and $c_{n_1,n_2}$ are the matrix elements obtained by time-averaging the steady-state solution of the master equation. We verify that increasing $N_{\mathrm{cutoff}}$ beyond ∼10 does not significantly alter the results within the quantum regime studied in this work, suggesting that the computed QMI converges to a within numerical precision.

To evaluate Equation (7), we first diagonalize the full two-oscillator density matrix ${\rho }_{12}$

(9)${\rho }_{12}=\sum _k{\lambda }_k\left|{\psi }_k\right.⟩\left.⟨{\psi }_k\right|,S\left({\rho }_{12}\right)=-\sum _k{\lambda }_{k}ln{\lambda }_k,$

where $\left\{{\lambda }_k\right\}$ and $\left\{|{\psi }_k⟩\right\}$ are the eigenvalues and eigenbases of ${\rho }_{12}$, respectively. We then obtain the reduced density matrices for each oscillator by tracing out the other as ${\rho }_1={\mathrm{Tr}}_2\left({\rho }_{12}\right),$ and ${\rho }_2={\mathrm{Tr}}_1\left({\rho }_{12}\right).$ Both ${\rho }_1$ and ${\rho }_2$ are similarly diagonalized to compute their von Neumann entropies, $S\left({\rho }_1\right)$ and $S\left({\rho }_2\right)$. Substituting these into Equation (7) gives the QMI. The resulting $I(1:2)$ quantifies the total correlations between the two oscillators, including both classical and quantum components, and thus complements a classical PID-based analysis.

2.3.3. Memory Capacity (MC)

In addition to instantaneous encoding, we are interested in how the system stores information over time, as is relevant in quantum reservoir computing (QRC). The memory capacity quantifies how well the current state of the reservoir (the outputs $X_1\left(t\right)$, $X_2\left(t\right)$) retains information about past inputs $s(t-\tau )$. By analyzing how the reconstruction error of past inputs varies with the delay $\tau$, we derive a memory measure that complements the PID and QMI analyses.

A high memory capacity suggests that the reservoir not only encodes the input at a given instant but also preserves information over extended periods. Comparing the memory capacity with $\mathrm{Rdn}$ and $\mathrm{Syn}$ reveals how different dynamical regimes influence both the instantaneous and temporal dimensions of information processing in this system. Details of the MC calculation are provided in Section 3.4.

3. Results and Discussion

We now present numerical evidence that coupled Kerr oscillators can encode input signals in either a redundant or synergistic fashion, depending on J, $\Delta$, and $\gamma$.

3.1. Emergence of Synergistic Encoding

We begin by examining how the joint mutual information (MI), $I\left(s:(X_1,X_2)\right)$, compares to the individual MIs $I(s:X_1)$ and $I(s:X_2)$ when probing our driven-dissipative system. We focus on two representative parameter sets: a mean-field regime with larger drive and smaller Kerr nonlinearities, and a quantum regime with smaller drive and larger Kerr nonlinearities. In both cases, we fix the detuning and damping at $\Delta =-2$ and $\gamma =0.5$. Concretely, in the mean-field case, we take $F=2.0$, $U_1=6.25\times {10}^{-3}$, and $U_2=2U_1$, while in the quantum regime, we take $F=0.2$, $U_1=4.0$, and $U_2=2U_1$.

3.1.1. Synergy from Total Mutual Information Consideration

Figure 2 compares $I\left(s:(X_1,X_2)\right)$ with $I(s:X_1)$ and $I(s:X_2)$, revealing that

$I\left(s:(X_1,X_2)\right)>I(s:X_1)+I(s:X_2),$

in the mean-field dynamics regime. This information excess indicates that measuring both $X_1$ and $X_2$ jointly can reveal strictly more information about the external drive signal $s\left(t\right)$ than either observable alone, suggesting a potential synergistic encoding mechanism.

3.1.2. Transition from Redundant to Synergistic Encoding

To further investigate whether this information surplus really originates from synergistic effects (rather than unique information in each oscillator), we perform Partial Information Decomposition (PID) according to Equation (6). For clarity, we normalize synergy and redundancy by $I\left(s:(X_1,X_2)\right)$, respectively,

${\mathrm{Syn}}_{\mathrm{norm}}={\displaystyle \frac{\mathrm{Syn}}{I\left(s:(X_1,X_2)\right)}},{\mathrm{Rdn}}_{\mathrm{norm}}={\displaystyle \frac{\mathrm{Rdn}}{I\left(s:(X_1,X_2)\right)}}.$

In Figure 3, we compare the normalized synergy (left) and redundancy (right) across three regimes (mean-field, second-order cumulant, and fully quantum) at fixed detuning and dissipation. As we sweep the coupling strength J from small to large, a pronounced synergy peak emerges near $J\approx |\Delta |=2$, marking a transition from predominantly redundant encoding to notably higher synergy (near $J\simeq |\Delta |$). In the quantum regime, stronger quantum correlations bias the encoding scheme slightly toward redundancy even near the peak. In contrast, the second-order cumulant approach captures the partial quantum correlations that lie between the mean-field regime (where correlations are suppressed) and the fully quantum regime (where all orders of correlations may appear). This second-order cumulant dynamics provides an approximate interpolation for regimes where quantum correlations are significant but not too strong (we set $\Delta =-2,F=0.5,\gamma =0.5,U_1=0.2,U_2=2U_1$ to represent a dynamical regime with non-negligible quantum correlations, motivating the use of second-order cumulants); see Appendix D.

3.2. Underlying Mechanisms Enhancing Synergistic Coding: The Role of Soft and Fast Modes

In this section, we explain the sharp increase in synergy observed near $J\simeq |\Delta |$, and attribute this behavior to the interplay between soft and fast modes in the coupled Kerr oscillators. Specifically, we demonstrate that the dominance of fast modes, enabled by the overdamping of soft modes, enhances coherent collective dynamics, leading to an increase in synergistic information.

3.2.1. Soft Modes and Potential Landscape Flatness

In the mean-field approximation without external drive ($F=0$), the coupled Kerr oscillators evolve in the following effective potential (see Appendix C), which captures the interplay of detuning, Kerr nonlinearities, and coupling

(10)$V\left({\overrightarrow{\alpha }}_c,{\overrightarrow{\alpha }}_c^{\ast }\right)=\sum _{j=1,2}\left(\Delta {\left|{\alpha }_{j,c}\right|}^2+\frac{1}{4}{U}_j{\left|{\alpha }_{j,c}\right|}^4\right)+J\left({\alpha }_{1,c}{\alpha }_{2,c}^{\ast }+{\alpha }_{2,c}{\alpha }_{1,c}^{\ast }\right),$

where ${\Delta }_1={\Delta }_2=\Delta$ and ${\overrightarrow{\alpha }}_c=[{\alpha }_{1,c},{\alpha }_{2,c}]$. At $J=|\Delta |$, the Hessian of V evaluated at the steady-state solution ${\alpha }_{1,c}={\alpha }_{2,c}=0$ develops zero eigenvalues, corresponding to flat or marginal directions. These flat directions represent soft modes, characterized by near-zero oscillation frequencies. Specifically, when evaluated at the steady state ${\alpha }_{1,c}={\alpha }_{2,c}=0$, the Hessian matrix of this effective potential (10) calculated in terms of the vector $(\mathrm{Re}\left({\alpha }_{1,c}\right),\mathrm{Im}\left({\alpha }_{1,c}\right),\mathrm{Re}\left({\alpha }_{2,c}\right),\mathrm{Im}\left({\alpha }_{2,c}\right))$ is

(11)$\mathbf{H}\left[V(\overrightarrow{\alpha },{\overrightarrow{\alpha }}^{\ast })\right]{|}_{\overrightarrow{\alpha },{\overrightarrow{\alpha }}^{\ast }=0}=2\left(\begin{array}{cccc}\Delta & 0& J& 0\\ 0& \Delta & 0& J\\ J& 0& \Delta & 0\\ 0& J& 0& \Delta \end{array}\right).$

The eigenvalues of this Hessian are

$\left\{2\right(\Delta +J),2(\Delta +J),2(\Delta -J),2(\Delta -J\left)\right\}.$

For the relevant parameter regime in which $\Delta <0$ and $J>0$, exactly when $J=|\Delta |$ that two flat directions (zero eigenvalues) emerge, and the other two directions are unstable at the second order (negative eigenvalues); see Figure 4.

3.2.2. Overdamping of Soft Modes at $J=|\Delta |$

Including dissipation with a rate $\gamma$ can transform the nearly flat directions into overdamped dynamics as follows. Following the Keldysh formalism [38,39,40] (see Appendix C), the retarded Green’s function shows poles of the form

(12)${\omega }_s=\pm \left|J-|\Delta |\right|-i\gamma ,{\omega }_f=\pm \left|J+|\Delta |\right|-i\gamma ,$

where ${\omega }_s$ (slow) and ${\omega }_f$ (fast) label the respective branches. Exactly at $J=|\Delta |$, the real part of ${\omega }_s$ vanishes, leaving only $-i\gamma$, indicating an overdamped relaxation to the steady state. In contrast, the fast modes remain oscillatory with frequencies $\mathrm{Re}\left({\omega }_f\right)=|J+|\Delta \left|\right|$. Consequently, at $J=|\Delta |$, the dynamics are dominated by the coherent oscillations of the fast modes, as the soft modes contribute only non-oscillatory relaxation.

3.2.3. Coherence-Driven Synergy Enhancement

When the dissipation rate is comparable to the oscillatory frequency of the fast modes ($\gamma \sim \mathrm{Re}\left({\omega }_f\right)$), the dominance of fast modes at $J=|\Delta |$ reduces competition between oscillatory frequencies (soft modes disappear and only fast modes persist) and enhances the coherence of the system’s collective response. More specifically, following the discussion in Appendix C, one can consider the relaxation dynamics of the small perturbation $\delta {\alpha }_c\left(t\right)$ around the steady state. It follows that the relaxation dynamics of each observable at site j can be expressed as

(13)$\mathrm{Re}\left(\delta {\alpha }_{j,c}\left(t\right)\right)=c_{j,s}{e}^{-\gamma t}cos\left(\mathrm{Re}\left({\omega }_s\right)t\right)+c_{j,f}{e}^{-\gamma t}cos\left(\mathrm{Re}\left({\omega }_f\right)t\right),$

where $c_{j,s}$ and $c_{j,f}$ are initial-perturbation-dependent mode amplitudes. When $\mathrm{Re}\left({\omega }_s\right)\to 0$, the soft mode contribution simplifies to pure exponential decay, and the dissipative dynamics are dominated by the single oscillatory frequency ${\omega }_f$ of the fast modes.

This coherence relaxation eliminates competing oscillations and minimizes overlapping contributions between modes, reducing redundancy. Also, the collective dynamics driven by the fast modes encode information that cannot be captured by any single oscillator alone, thereby enhancing synergistic information. This explains the peak in synergy observed near $J\simeq |\Delta |$.

Figure 5 illustrates how the slow-mode poles (orange) move to the imaginary axis at $J=|\Delta |$, marking the disappearance of competing oscillatory timescales. In this near-critical regime, the relaxation dynamics become dominated by the underdamped (oscillatory) contributions of the fast modes, resulting in coherent dissipation. It is important to note that this result pertains to the regime where $\gamma \sim |\mathrm{Re}\left({\omega }_f\right)|\gg |\mathrm{Re}\left({\Omega }_s\right)|$.

3.2.4. Generality of the Transition to Synergistic Behavior

The observed synergy peak at $J\simeq |\Delta |$ is not specific to the type of input signal driving the system. To demonstrate this, we perform numerical simulations of the master equation describing quantum dynamics with the input signal $s\left(t\right)$ sampled from a uniform distribution in the interval $[-1,1]$, and uncorrelated in time. The results, shown in Figure 6, confirm that the transition in redundant/synergistic behavior persists at $J=|\Delta |$, regardless of the statistical properties of the input.

These results emphasize that the transition to synergistic encoding at $J\simeq |\Delta |$ is an intrinsic feature of the system’s response dynamics, driven by the dominance of underdamped fast modes, and not much by the input signal properties. In the following section, we explore how increasing $\gamma$ impacts the system’s encoding behavior, showing that fast relaxation towards the steady state progressively shifts the system from synergistic to redundant encoding.

3.3. Large Dissipation Leads to Redundant Encoding

As the damping rate $\gamma$ increases, the dynamics progressively shift toward overdamped relaxation for both slow and fast modes. This transition leads to a steady-state regime where the two cavities become nearly identical, resulting in redundant coding of the input information. While the emergence of soft modes and the dominance of fast modes at $J\simeq |\Delta |$ enhance synergy, increasing $\gamma$ gradually suppresses this effect, shifting the system toward a more redundant encoding regime dominated by rapid overdamped relaxation towards the steady state.

This behavior is particularly evident in the quantum regime, where the system approaches a product state at large $\gamma$, rendering the two oscillators effectively independent. As shown in Figure 7, increasing $\gamma$ from the baseline value of $\gamma =0.5$ (as seen in Figure 2 and Figure 3) results in a clear reduction in quantum correlations. Panel (a) of Figure 7 quantifies this through the time-averaged quantum mutual information (QMI), which decreases as $\gamma$ increases. Notably, peaks in QMI at $J\simeq |\Delta |$ coincide with peaks in classical mutual information (MI) as shown in panel (b). This observation aligns with the intuition that higher quantum correlations often translate to enhanced classical information encoding near the critical coupling.

Figure 7c,d further illustrate this transition by showing how absolute synergy diminishes and absolute redundancy grows with increasing $\gamma$. At low $\gamma$, the dynamics favor (weakly) synergistic encoding. At high $\gamma$, however, the system becomes dominated by (highly) redundant encoding, with both cavities responding similarly and independently of each other.

3.4. Memory Capacity of Synergistic and Redundant Encoding

We close our discussion by analyzing the performance of the coupled Kerr oscillators as a quantum reservoir, focusing on their capacity to retain and process temporal information. This memory capacity benchmark highlights how the synergistic and redundant behaviors identified earlier influence practical tasks such as time-series memorization.

3.4.1. Short-Term Memory Task

To quantify memory capacity, we train the system to recall past input signals using a short-term memory task. The input signal $s\left(t\right)$ is sampled from a uniform distribution in the interval $[-1,1]$ and is uncorrelated in time, and the target time series ${\overline{y}}_n\left(t\right)$ corresponds to the input signal at a previous time step:

(14)${\overline{y}}_n\left(t\right)=s(t-n\Delta t),$

where $\Delta t=0.01$ is the time step. The output observables $X_i\left(t\right)=\mathrm{Re}\left(⟨{\widehat{a}}_i\left(t\right)⟩\right)$ and $Y_i\left(t\right)=\mathrm{Im}\left(⟨{\widehat{a}}_i\left(t\right)⟩\right)$ are used as feature vectors (here, we use more output readouts than in the previous sections since this input signal is more difficult to fit with fewer feature vectors). We construct an output vector $\overrightarrow{X}\left(t\right)$:

(15)$\overrightarrow{X}\left(t\right)={(X_1\left(t\right),X_2\left(t\right),Y_1\left(t\right),Y_2\left(t\right),1)}^{\top },$

and fit the target ${\overline{y}}_n\left(t\right)$ using a weight matrix $\mathbf{W}$ via standard linear regression with Tikhonov regularization:

(16)$\widehat{y}\left(t\right)={\mathbf{W}}^{\ast }\overrightarrow{X}\left(t\right),$

where the optimal weights ${\mathbf{W}}^{\ast }$ minimize the mean square error (MSE) during training:

(17)$\mathrm{MSE}(\mathbf{W},\lambda )=\left|\right|\overline{\mathbf{y}}-\widehat{\mathbf{y}}{\left|\right|}^2+{\lambda \left|\right|\mathbf{W}\left|\right|}^2.$

For simplicity, we set $\lambda =0$. After training, the memory capacity [1] for delay step n is evaluated as

(18)$\mathrm{MC}\left(n\right)={\displaystyle \frac{{\mathrm{cov}}^2({\overline{\mathbf{y}}}_n,{\widehat{\mathbf{y}}}_{\mathrm{pred}})}{{\sigma }^2\left({\overline{\mathbf{y}}}_n\right){\sigma }^2\left({\widehat{\mathbf{y}}}_{\mathrm{pred}}\right)}},$

where $0\le \mathrm{MC}\left(n\right)\le 1$, with $\mathrm{MC}\left(n\right)=1$ indicating perfect recall.

3.4.2. Tuning J Towards the Critical Coupling Strength

Figure 8 shows the memory capacity as J approaches the critical coupling $J=|\Delta |=2$. At smaller values of J, the memory capacity remains low across the delay steps. However, as J approaches $|\Delta |$, the system exhibits a notable change in behavior: the memory capacity for short delays (small n) increases significantly, while the memory capacity for longer delays (large n) decays more rapidly. This behavior reflects the trade-off between short-term and long-term memory as the system transitions to the synergistic regime dominated by fast modes.

We note that, for low dissipation (e.g., $\gamma =0.5$), the system may appear to retain information about initial conditions at long delays n, indicating that the fading memory property might not be fully realized. While a more comprehensive benchmark such as the information processing capacity (IPC) [41] could more definitively confirm or refute fading memory in this low-dissipation regime, our current goal is to investigate synergy and redundancy in a small quantum reservoir rather than to perform an exhaustive reservoir-computing benchmark. Instead, we emphasize qualitative trends in short- vs. long-term memory. A more detailed, IPC-based study would be an interesting direction for future work to fully study the computational potential of near-critical quantum reservoirs in this system.

3.4.3. Impact of Dissipation $\gamma$

Figure 9 illustrates how dissipation ($\gamma$) affects memory capacity at the critical coupling $J=|\Delta |=2$. As $\gamma$ increases, the memory capacity can be attributed to the enhanced stability (more rapid relaxation towards) in the reservoir’s steady-state dynamics. Dissipation suppresses oscillatory behavior and stabilizes the reservoir’s response. This stabilization corresponds to a regime of highly redundant encoding, where information is stored near the steady state across somewhat identical subsystems.

Interestingly, as shown in Figure 8, while higher $\gamma$ leads to improved total memory capacity, the decay rate of memory capacity in this redundant regime, $\mathrm{MC}\left(n\right)\sim exp(-\Gamma n)$, remains approximately constant across different dissipation rates. This suggests that dissipation uniformly governs the loss of correlations over time. Differences in memory capacity at high $\gamma$ primarily arise from the proportionality factor in Equation (18). As dissipation increases, the variance in the reservoir’s output prediction decreases, reflecting rapid stabilization to the steady state. This reduced variance amplifies the initial memory capacity but does not alter the exponential decay rate of correlations.

These results highlight a subtle interplay between dissipation, memory capacity, and encoding modes. Near criticality at $J=|\Delta |$, the system exhibits synergistic behavior, where collective dynamical response dominates, enabling the reservoir to respond sensitively to recent input signals. This synergistic encoding boosts short-term memory retention but comes at the expense of long-term storage, as the system’s ability to retain correlations with far past inputs diminishes rapidly due to sensitivity to perturbation.

At higher dissipation rates $\gamma$, the system tends to encode inputs more redundantly and relaxes more quickly to a steady state. Interestingly, as illustrated in Figure 9, both short- and long-term memory capacities, at criticality, increase with $\gamma$ for a simple uniformly random input memorization task, even though one might initially expect a trade-off. The underlying reason could be that the system has not yet fully achieved the fading memory regime, so further investigations, potentially through more comprehensive benchmarks such as IPC analysis, are needed to confirm how dissipation precisely balances short-term responsiveness and long-term stability right at a critical point.

Also from the perspective of quantum reservoir design, dissipation plays a critical role in ensuring the reservoir exhibits the fading memory property, where the influence of past inputs gradually diminishes, which is a necessary condition for designing operational reservoir computing [9,10].

4. Conclusions and Outlook

In this work, we explored coupled Kerr-nonlinear oscillators as a model open quantum system for studying modes of information processing in a small quantum reservoir computing (QRC) platform. By employing Partial Information Decomposition (PID) and analyzing memory retention tasks, we investigated how the interplay of dynamical instability and dissipation governs the encoding of input information into redundant or synergistic modes. These encoding schemes play a crucial role in determining the reservoir’s performance in processing and retaining temporal data.

Our findings reveal several key insights. First, near the critical coupling strength $J=|\Delta |$, the system transitions from predominantly redundant to synergistic encoding. This transition is driven by the dynamics of coherent oscillation modes that dominate as slow modes (soft modes) become overdamped. These collective dynamics enable the reservoir to process information synergistically, boosting short-term memory retention. Importantly, this synergy is robust across different input signals, including telegraph processes and uniform uncorrelated noise, suggesting that the observed transition is an intrinsic feature of the system’s bifurcation near $J=|\Delta |$. This dynamic instability was elucidated through the non-equilibrium (Keldysh) field-theoretic analysis, which highlighted how the disappearance of soft modes amplifies the dominance of fast, coherent modes.

Dissipation ($\gamma$) plays an important role in shaping information encoding and memory capacity for our reservoir near criticality. At low $\gamma$, synergistic encoding enables collective processing and enhances short-term memory retention, as the reservoir’s dynamics are sensitive to recent input signals. However, as $\gamma$ increases, dissipation suppresses coherent dynamics, rapidly driving the system toward redundant encoding. In this regime, large dissipation stabilizes encoding by enforcing redundant representations near the steady state, enabling the reservoir to retain information about far-past inputs at the expense of sensitivity to recent input information.

The connection between encoding modes and memory capacity reveals a trade-off: synergistic encoding favors short-term memory retention but is less suited for long-term storage due to its sensitivity to perturbations. Redundant encoding, on the other hand, sacrifices responsiveness to recent inputs but improves long-term retention. From a quantum reservoir design perspective, dissipation ensures the reservoir exhibits the fading memory property necessary for effective reservoir design [9,42]. By carefully tuning dissipation, one can hope to optimize the balance between short-term responsiveness and long-term stability and tailor the reservoir to specific computational tasks, aligning with the recent work on engineered dissipation as computational resources in quantum systems [10,43,44].

Coupled Kerr-nonlinear oscillators provide a minimal yet insightful testbed for analyzing transitions between encoding modes driven by dissipation and dynamical instability. Extending this framework to other quantum systems with dynamical phase transitions, such as Bose–Hubbard models [40], driven-dissipative platforms [45], and spin systems [13,16], could deepen our understanding of how dissipation and instability shape encoding dynamics in systems with more complex phase spaces. Another interesting direction is to develop a rigorous definition of quantum synergy. Although this work uses Partial Information Decomposition (PID) to analyze classical information derived from quantum observables, incorporating quantum correlations could provide a more comprehensive view of how coherence and other quantum effects either enhance or constrain information encoding, in line with [7,8,13]. Such development would refine our understanding of synergy and redundancy in small quantum systems from a more quantum information theoretic viewpoint.

In practice, realistic quantum devices inevitably face various noise sources and control imperfections. While these perturbations may shift the crossover point between redundant and synergistic encoding or relocate the system’s critical point, our results show that this crossover in the modes of encoding persists over a wide parameter range. With advanced experimental techniques, such as dynamical decoupling and error mitigation in quantum optics, we expect that the essential physics behind this crossover can still be realized in real experimental systems.

Lastly, it is important to recognize that the results presented here, like much of the prior work, assume perfect readout without the presence of shot noise. As system sizes scale, however, readout processes may suffer from exponential concentration phenomena, requiring exponentially many measurement shots to estimate input-dependent readouts accurately as discussed in [11]. This limitation presents a significant barrier to scalability. In larger quantum reservoir systems, future studies must address how synergy and redundancy behave under realistic measurement constraints. Incorporating the effects of finite measurement precision into the framework of encoding dynamics could lead to more scalable and experimentally feasible designs for quantum reservoirs.

Author Contributions

Conceptualization, methodology, investigation, formal analysis, writing, K.C. and T.C.; software, K.C.; funding acquisition, T.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The codes and datasets used in this work are available from the corresponding author upon reasonable request.

Acknowledgments

T.C. acknowledges insightful discussions with Wave Ngampruetikorn on synergistic and redundant coding.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Schematic of two coupled Kerr-nonlinear oscillators. Each cavity i features a Kerr nonlinearity [Forumla omitted. See PDF.] and a photon loss rate [Forumla omitted. See PDF.]. A time-dependent drive [Forumla omitted. See PDF.] (green arrows) injects identical signals into both cavities, while the coherent tunneling of strength J (violet arrow) couples the two modes. We measure the mean fields [Forumla omitted. See PDF.] to probe the system’s response. The Hamiltonian is specified by Equation (1), and the Lindblad Equation (2) governs this driven-dissipative dynamics. We assume both cavities have the same detuning [Forumla omitted. See PDF.] from the drive frequency. This work investigates how the readouts [Forumla omitted. See PDF.] encode the time-dependent drive [Forumla omitted. See PDF.] across different dynamical regimes of the coupled Kerr oscillators.

Figure 2. Classical mutual information [Forumla omitted. See PDF.], compared to [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] in (Left) a mean-field regime and (Right) a quantum regime. In the mean-field regime, [Forumla omitted. See PDF.] exceeds [Forumla omitted. See PDF.] or [Forumla omitted. See PDF.] alone near [Forumla omitted. See PDF.], hinting at synergy. On the other hand, in the quantum regime, [Forumla omitted. See PDF.] is comparable to, but not always exceeding, the sum of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 3. (Left) Normalized synergy vs. the coupling J. (Right) Normalized redundancy vs. the coupling J. A pronounced peak near [Forumla omitted. See PDF.] marks the crossover from predominantly redundant to more synergistic encoding. In the fully quantum description, enhanced quantum correlations can favor redundancy even at the transition, whereas second-order cumulants interpolate between mean-field and quantum descriptions.

Figure 4. Effective potential around the steady state [Forumla omitted. See PDF.] (red dot) in the [Forumla omitted. See PDF.] regime, projected onto [Forumla omitted. See PDF.]. (Left) When [Forumla omitted. See PDF.], the steady state is weakly unstable in all directions, with no soft modes present, and the system predominantly encodes inputs redundantly. (Center) At the critical point [Forumla omitted. See PDF.], flat directions appear, marking the onset of soft modes. In this near-critical regime, collective oscillations enhance synergistic encoding. (Right) For [Forumla omitted. See PDF.], the potential deforms into a saddle, with two stable and two unstable directions. Here, the soft modes again disappear, and the system encodes inputs redundantly.

Figure 5. Retarded Green’s function poles (A16) in the complex-frequency plane as J increases, illustrating the evolution of slow modes [Forumla omitted. See PDF.] (orange dots) and fast modes [Forumla omitted. See PDF.] (blue dots). For [Forumla omitted. See PDF.], both slow and fast modes coexist, and the system tends to encode inputs more redundantly. Near the critical point [Forumla omitted. See PDF.], the real part of [Forumla omitted. See PDF.] approaches zero, indicating a disappearance of slow collective oscillations and an overdamped decay. In this near-critical regime, collective oscillations due to underdamped fast modes dominate and enhance synergistic encoding. For [Forumla omitted. See PDF.], the slow modes shift away from zero frequency, reducing synergy and transitioning the system back toward redundant encoding. This interplay between slow and fast modes governs how the system transitions from predominantly redundant to synergistic processing and back again as J crosses the critical point.

Figure 6. Impact of uniform, uncorrelated noise input on information encoding at the quantum regime ([Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]). (Left) Normalized PID with uniform noise input. (Right) Comparison of total MI and partial MI contributions. The left panel shows a clear transition to synergistic encoding near [Forumla omitted. See PDF.]. The right panel compares [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], highlighting the emergence of synergy near the critical coupling.

Figure 7. Impact of increasing [Forumla omitted. See PDF.] on information metrics in the quantum regime ([Forumla omitted. See PDF.], [Forumla omitted. See PDF.]). (a) Time-averaged quantum mutual information (QMI) between the two oscillators as a function of J for different [Forumla omitted. See PDF.]. (b) Classical mutual information (MI) between the input signal and the outputs vs. J. (c) Absolute synergy vs. J. (d) Absolute redundancy vs. J. These plots illustrate the transition from low-synergistic to high-redundant encoding with increasing [Forumla omitted. See PDF.]. At large dissipation, the two oscillators approach a product state, indicated by low QMI. Interestingly, despite the redundancy dominating at higher [Forumla omitted. See PDF.], the total mutual information near [Forumla omitted. See PDF.] remains approximately constant.

Figure 8. (Left) Memory capacity [Forumla omitted. See PDF.] for n = 1–10 and [Forumla omitted. See PDF.], showing an increase in short-term memory capacity as [Forumla omitted. See PDF.]. (Right) Memory capacity for n = 1–20 and [Forumla omitted. See PDF.], showing that the long-term capacity drops as [Forumla omitted. See PDF.]. Parameters: [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.]. These results are averaged over 50 input realizations.

Figure 9. Memory capacity at the critical coupling [Forumla omitted. See PDF.] as [Forumla omitted. See PDF.] increases. (Left) Linear scale plot: [Forumla omitted. See PDF.] vs. n. (Right) Log scale plot of [Forumla omitted. See PDF.]. The left panel shows the average memory capacity (100 realizations). The right panel illustrates an exponential decay in [Forumla omitted. See PDF.]. The exponential decay rate of the two-time correlation in the memory capacity, [Forumla omitted. See PDF.], remains approximately constant as [Forumla omitted. See PDF.] increases, indicating that dissipation uniformly governs the loss of correlations across different dissipation rates. Differences in memory capacity primarily arise from the proportionality factor, with larger [Forumla omitted. See PDF.] leading to a smaller variance in the output prediction in the denominator of Equation (18), as the output rapidly stabilizes to the steady state. This rapid stabilization corresponds to highly redundant encoding, and in turn enhances the total memory capacity for memorizing uniformly random input time series. Notably, in this redundant coding regime, the highly dissipative dynamics improve the quantum reservoir’s memory capacity.

Appendix A. Partial Information Decomposition (PID)

In classical (Shannon) information theory, we run into conceptual difficulties as soon as there are more than two random variables to handle because there is no single, universally accepted way to break down the total information among multiple variables following Shannon’s original recipe [46]. Quantifying the amount of information that is shared (redundant), exclusive (unique), or complementary (synergistic) among three or more variables is highly nontrivial and remains an active area of research; see, for example, [47,48,49,50,51,52]. In this work, we adopt one particular framework that meets three criteria: (1) It has a relatively straightforward formalism. (2) It captures whether and how our system is redundant or synergistic. (3) It does not exhibit major interpretational drawbacks for our purposes.

This framework, called Partial Information Decomposition (PID) [53], attempts to disentangle the multivariate information into non-overlapping, non-negative parts with clear interpretations. Concretely, consider a set of n random variables $\{S,R_1,R_2,\dots ,R_{n-1}\}$, where S (the target) is the variable whose information we want to capture, and $\mathbf{R}=\{R_1,\dots ,R_{n-1}\}$ (the sources) is the combined set of variables providing that information. The total (multivariate) mutual information between the source and the target is (for discrete distributions)$$I(S:\mathbf{R})=\sum _{s,r_1,\dots ,r_{n-1}}{P}_{(S,\mathbf{R})}(s,r_1,\dots ,r_{n-1})log\left({\displaystyle \frac{P_{(S,\mathbf{R})}(s,r_1,\dots ,r_{n-1})}{P_S\left(s\right)P_{\mathbf{R}}(r_1,\dots ,r_{n-1})}}\right).$$

As outlined in the original PID paper [53], this methodology can in principle be generalized to any number of variables in terms of the PID lattice, but its complexity increases rapidly once $n>3$. (Further developments can be found in [54,55,56,57,58].) In this work, we focus on the simplest nontrivial case $n=3$, for which the PID formalism is most concretely developed and comparatively well understood.

Appendix A.1. PID for Three Variables

Let $X,Y,Z$ be three random variables, and consider $I(X:(Y,Z\left)\right)$, the mutual information between X and the pair $(Y,Z)$. The three-variable PID proposes the following decomposition of $I(X:(Y,Z\left)\right)$ into four distinct parts:(A1)$$I(X:(Y,Z\left)\right)=\mathrm{Rdn}(X:Y;Z)+\mathrm{Syn}(X:Y;Z)+\mathrm{Unq}(X:Y\setminus Z)+\mathrm{Unq}(X:Z\setminus Y),$$ along with corresponding decompositions of the pairwise mutual informations: (A2)$$\begin{array}{cc}\hfill I(X:Y)=& \mathrm{Rdn}(X:Y;Z)+\mathrm{Unq}(X:Y\setminus Z),\hfill \end{array}$$(A3)$$\begin{array}{cc}\hfill I(X:Z)=& \mathrm{Rdn}(X:Y;Z)+\mathrm{Unq}(X:Z\setminus Y).\hfill \end{array}$$

Here, the four partial information (PI) terms have the following interpretations:

• $\mathrm{Rdn}(X:Y;Z)$ (redundancy/shared information): the amount of information about X that is found in common in Y and Z.

• $\mathrm{Syn}(X:Y;Z)$ (synergy/complementary information): the information about X that is only accessible when considering Y and Z jointly, i.e., the “whole is greater than the sum of its parts”, which is typically a signature of the cooperation of constituents in complex systems.

• $\mathrm{Unq}(X:Y\setminus Z),\mathrm{Unq}(X:Z\setminus Y)$ (unique information): the information about X that can be acquired solely from Y (or Z) and not from Z (or Y).

Because there are four unknown PI quantities but only three equations, Equations (A1)–(A3), one cannot solve for the PIs without additional theoretical constraints. Different PID axioms or definitions fix one of these quantities first (e.g., by prescribing a formula or an inequality), allowing the remaining quantities to be determined from the joint distribution $P(X,Y,Z)$. Approaches in the literature include (1) defining redundancy first [49,53,57,59], (2) defining synergy first [50,60], and (3) defining unique information first [23].

Appendix A.2. Co-Information

From Equations (A1)–(A3), one obtains the following notable combination of mutual information:(A4)$$I(X:Y)+I(X:Z)-I\left(X:(Y,Z)\right)=\mathrm{Rdn}(X:Y;Z)-\mathrm{Syn}(X:Y;Z).$$

The left-hand side is often called the co-information [47] (or equivalently, interaction information [61] or sometimes just the mutual information [62]):(A5)$$\mathrm{CoI}(X;Y;Z)=I(X:Y)+I(X:Z)-I\left(X:(Y,Z)\right).$$

Although this is the simplest and most direct extension of mutual information to three variables, it has two important drawbacks. First, it can take positive or negative values (making some interpretations more subtle than the non-negative mutual information). Secondly, it cannot really distinguish redundancy and synergy because the right-hand side of (A4) just reflects their differences.

Nevertheless, for three-variable systems, CoI can sign whether redundancy or synergy dominates. Even though the co-information cannot quantify the amount of redundancy and synergy separately, it can describe whether the system is in a redundant coding or synergistic coding regime [21]. When $\mathrm{CoI}>0$, the system is said to be redundancy dominated, and when $\mathrm{CoI}<0$, the system is said to be synergy dominated [51,63]. This measure also plays a role in certain PID definitions and algorithms [23].

Appendix B. Explicit Examples for PID Calculation

In this appendix, we provide simple but useful examples to show how one can compute Partial Information Decomposition (PID) terms analytically under the BROJA approach [23]. The BROJA method obtains each partial information by solving a constrained optimization problem over all joint distributions $Q(X,Y,Z)$, consistent with certain marginal constraints derived from the true distribution $P(X,Y,Z)$.

Appendix B.1. General Framework

Let $X,Y,Z$ be three discrete random variables with a true joint distribution $P(X,Y,Z)$. The BROJA definitions of the four partial information terms ($\mathrm{Unq},\mathrm{Rdn},\mathrm{Syn}$) are as follows [23]:(A6)$$\begin{array}{cc}\hfill \mathrm{Unq}(X:Y\setminus Z)& =\underset{Q\in {\Delta }_P}{min}{I}_Q(X:Y\mid Z),\hfill \\ \hfill \mathrm{Unq}(X:Z\setminus Y)& =\underset{Q\in {\Delta }_P}{min}{I}_Q(X:Z\mid Y),\hfill \\ \hfill \mathrm{Rdn}(X:Y;Z)& =\underset{Q\in {\Delta }_P}{max}{\mathrm{CoI}}_Q(X;Y;Z),\hfill \\ \hfill \mathrm{Syn}(X:Y;Z)& =I(X:(Y,Z))-\underset{Q\in {\Delta }_P}{min}{I}_Q(X:(Y,Z)),\hfill \end{array}$$ where ${\mathrm{CoI}}_Q(X;Y;Z)$ is the co-information under the joint distribution Q, and$$I(X:(Y,Z))=I_P(X:(Y,Z))$$ is the mutual information computed from the actual distribution P. Any subscripted quantity such as $I_Q(·)$ is computed from a candidate joint distribution $Q\in {\Delta }_P$, not from the true P.

Appendix B.2. Definition of Δ P

The space ${\Delta }_P\subseteq \Delta$ consists of all joint distributions $Q(X,Y,Z)$ whose $(X,Y)$ and $(X,Z)$ marginals match those of P. Formally, if $\Delta$ is the set of all distributions on $\mathcal{X}\times \mathcal{Y}\times \mathcal{Z}$, then(A7)$${\Delta }_P=\left\{Q\in \Delta \left|Q\right(x,y)=P(x,y\left)\mathrm{and}Q\right(x,z)=P(x,z)\right\},$$ for all $(x,y,z)\in \mathcal{X}\times \mathcal{Y}\times \mathcal{Z}$. In other words, we fix the two-dimensional marginals of $(X,Y)$ and $(X,Z)$ to match the true data but allow $Q(Y,Z|X)$ to vary.

Appendix B.3. Foliation into Slices

To handle the high dimensionality of ${\Delta }_P$, a convenient approach parameterizes ${\Delta }_P$ via a foliation; see Appendix A of [23]. Concretely, for each x with $P\left(x\right)>0$, we define(A8)$${\Delta }_{P,x}=\left\{\tilde{Q}(Y,Z)\in \Delta (\mathcal{Y}\times \mathcal{Z})|\tilde{Q}\left(y\right)=P\left(y\right|x),\tilde{Q}\left(z\right)=P\left(z\right|x)\right\}.$$

A joint distribution $Q\in {\Delta }_P$ can then be specified by choosing, for each x, a $\tilde{Q}(Y,Z)\in {\Delta }_{P,x}$, and combining via(A9)$$Q(x,y,z)=P\left(x\right){\tilde{Q}}_x(y,z)=P\left(x\right)Q(y,z|x).$$

Hence, ${\Delta }_P$ can be viewed as the product of slice spaces ${\Delta }_{P,x}$ for all x. Once a suitable parameterization of these slices is chosen, the optimization in (A6) becomes a finite-dimensional convex optimization problem. Next, we show how the procedure is performed in the simple cases following [23].

Appendix B.4. Example 1: AND Gate

Consider three binary variables $\mathcal{X}=\mathcal{Y}=\mathcal{Z}=\{0,1\}$ with$$X=\left(Y\mathrm{AND}Z\right)$$ and $Y,Z$ independent and identically distributed $\mathrm{Bernoulli}(1/2)$. The true joint distribution P is uniform on the events $(X,Y,Z)\in \left\{\right(0,0,0),(0,0,1),(0,1,0),(1,1,1\left)\right\}$ with probability $1/4$ each, and zero otherwise. From this P, one derives the marginals $P(X,Y)$ and $P(X,Z)$.

Appendix B.5. Parametrizing Δ P

Applying the slice formalism, one obtains

• For $x=0$ (which has $P(X=0)=3/4$), the slice ${\Delta }_{P,0}$ is a 1-parameter family$${\tilde{Q}}_0(0,0)=\frac{1}{3}+{\alpha }^{\prime },{\tilde{Q}}_0(0,1)=\frac{1}{3}-{\alpha }^{\prime },{\tilde{Q}}_0(1,0)=\frac{1}{3}-{\alpha }^{\prime },{\tilde{Q}}_0(1,1)={\alpha }^{\prime },$$ with $0\le {\alpha }^{\prime }\le 1/3$.

• For $x=1$ (which has $P(X=1)=1/4$), the slice ${\Delta }_{P,1}$ is trivial because the only consistent distribution is ${\tilde{Q}}_1(1,1)=1$.

Combining ${\tilde{Q}}_0$ and ${\tilde{Q}}_1$ as in (A9) and reparameterizing yields the 1-parameter family$$Q_{\alpha }(x,y,z)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{4}}+\alpha ,\hfill & (x,y,z)=(0,0,0),\hfill \\ {\displaystyle \frac{1}{4}}-\alpha ,\hfill & (x,y,z)=(0,0,1),(0,1,0),\hfill \\ \alpha ,\hfill & (x,y,z)=(0,1,1),\hfill \\ {\displaystyle \frac{1}{4}},\hfill & (x,y,z)=(1,1,1),\hfill \end{array}\right.0\le \alpha \le {\displaystyle \frac{1}{4}}.$$

Under each candidate $Q_{\alpha }$, one can compute ${\mathrm{CoI}}_{Q_{\alpha }}(X;Y;Z)$ and $I_{Q_{\alpha }}(X:(Y,Z))$, and thus solve the optimization problems over $\alpha$ in (A6). Figure A1 shows ${\mathrm{CoI}}_{Q_{\alpha }}$ and $I_{Q_{\alpha }}$ vs. $\alpha$.

Figure A1. Co-information [Forumla omitted. See PDF.] and mutual information [Forumla omitted. See PDF.] for the AND-gate example, plotted as functions of [Forumla omitted. See PDF.]. The optimum for redundancy (respectively, synergy) occurs at [Forumla omitted. See PDF.].

Solving these optimizations shows that $\alpha =\frac{1}{4}$ is the critical point. Hence, each partial information is$$\begin{array}{cc}\hfill \mathrm{Unq}(X:Y\setminus Z)& =\mathrm{Unq}(X:Z\setminus Y)=0\hfill \\ \hfill \mathrm{Rdn}(X:Y;Z)& ={\displaystyle \frac{3}{4}}log{\displaystyle \frac{4}{3}}=0.311log2=0.311\mathrm{bit}\hfill \\ \hfill \mathrm{Syn}(X:Y;Z)& ={\displaystyle \frac{1}{2}}log2=0.5\mathrm{bit}\hfill \end{array}$$

Appendix B.6. Example 2: XOR Gate

Next, consider the binary XOR relation$$X=\left(Y\mathrm{XOR}Z\right),$$ where $Y,Z$ are again independent and identically distributed $\mathrm{Bernoulli}(1/2)$. The true distribution P is uniform on $\left\{\right(0,0,0),(1,0,1),(1,1,0),(0,1,1\left)\right\}$. One again sets up the slices ${\Delta }_{P,0}$ and ${\Delta }_{P,1}$, each giving a family of distributions parameterized by ${\alpha }^{\prime }$ and ${\beta }^{\prime }$.

• ${\Delta }_{P,0}$: the most general distribution that satisfies the constraint (A8) is$${\tilde{Q}}_0(y,z)=\left\{\begin{array}{cc}{\alpha }^{\prime }\hfill & ,(y,z)=(0,0)\hfill \\ {\displaystyle \frac{1}{2}}-{\alpha }^{\prime }\hfill & ,(y,z)=(0,1)\hfill \\ {\displaystyle \frac{1}{2}}-{\alpha }^{\prime }\hfill & ,(y,z)=(1,0)\hfill \\ {\alpha }^{\prime }\hfill & ,(y,z)=(1,1)\hfill \end{array}\right.$$ with $0\le {\alpha }^{\prime }\le {\displaystyle \frac{1}{2}}$.

• ${\Delta }_{P,1}$: since the distribution for XOR operation is symmetric under swapping Y and Z, then it turns out that$${\tilde{Q}}_1(y,z)=\left\{\begin{array}{cc}{\beta }^{\prime }\hfill & ,(y,z)=(0,0)\hfill \\ {\displaystyle \frac{1}{2}}-{\beta }^{\prime }\hfill & ,(y,z)=(0,1)\hfill \\ {\displaystyle \frac{1}{2}}-{\beta }^{\prime }\hfill & ,(y,z)=(1,0)\hfill \\ {\beta }^{\prime }\hfill & ,(y,z)=(1,1)\hfill \end{array}\right.$$ with another free parameter, $0\le {\beta }^{\prime }\le {\displaystyle \frac{1}{2}}$.

Combining them with $P(X=0)=P(X=1)=\frac{1}{2}$ yields a two-parameter family$$Q_{\alpha ,\beta }(x,y,z)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{8}}+\alpha \hfill & ,(x,y,z)=(0,0,0),(0,1,1)\hfill \\ {\displaystyle \frac{1}{8}}-\alpha \hfill & ,(x,y,z)=(0,0,1),(0,1,0)\hfill \\ {\displaystyle \frac{1}{8}}+\beta \hfill & ,(x,y,z)=(1,0,0),(1,1,1)\hfill \\ {\displaystyle \frac{1}{8}}-\beta \hfill & ,(x,y,z)=(1,0,1),(1,1,0)\hfill \end{array}\right.$$ where $-{\displaystyle \frac{1}{8}}\le \alpha ,\beta \le {\displaystyle \frac{1}{8}}$, governing all $Q\in {\Delta }_P$. As before, one computes ${\mathrm{CoI}}_{Q_{\alpha ,\beta }}(X;Y;Z)$ and $I_{Q_{\alpha ,\beta }}(X:(Y,Z))$ to solve the optimization in (A6). Figure A2 shows the surfaces of CoI and I vs. $\alpha ,\beta$.

Figure A2. (Left) [Forumla omitted. See PDF.]; (Right) [Forumla omitted. See PDF.]; for the XOR-gate example. The optimum occurs at [Forumla omitted. See PDF.].

The optimum occurs at $(\alpha ,\beta )=(0,0)$, giving$$\begin{array}{cc}\hfill \mathrm{Unq}(X:Y\setminus Z)& =\mathrm{Unq}(X:Z\setminus Y)=0\hfill \\ \hfill \mathrm{Rdn}(X:Y;Z)& =0\hfill \\ \hfill \mathrm{Syn}(X:Y;Z)& =log2=1\mathrm{bit}.\hfill \end{array}$$

This perfectly aligns with the well-known fact that XOR is synergistic. Consider $X=Y\mathrm{XOR}Z$, and we happen to have only the knowledge of Y, where either we know that $Y=0$ or $Y=1$. Without the knowledge of Z, we cannot infer any useful information about the value of X at all, and the probability $P\left(X\right|Y)$ is always equal to $1/2$ for any $(x,y)$. Only when both Y and Z are presented together do we know exactly what the value of X should be, and it is exactly one bit of information (i.e., one yes/no question) that we need to know in order to eliminate all uncertainty about the value of X.

These two classic gates (AND and XOR) show how the BROJA optimization can be computed analytically in simple discrete cases. In practice, for larger or more complex systems, numerical methods are necessary, but the underlying principle is the same. One restricts to ${\Delta }_P$ to preserve certain marginals and then solves the convex optimization problems in (A6).

Appendix C. Keldysh Action

In this appendix, we outline how to obtain the Keldysh action for the coupled Kerr oscillators described by Equations (1) and (2), and how the system’s linear response to external perturbation can be obtained from the retarded Green’s function. Readers seeking broader context on Keldysh formalism may consult [24,64,65].

From the Linblad Equation (1), one can reformulate this open-system evolution in a path-integral language by introducing “classical” fields ${\alpha }_{j,c}$ and “quantum” fields ${\alpha }_{j,q}$, capturing forward and backward time contours, respectively [24]. After performing the usual Keldysh rotation, the action $S_K$ takes the following form:(A10)$$\begin{array}{cc}\hfill S_K=\int dt[& \sum _{j=1,2}\left({\alpha }_{j,c}^{\ast }i{\partial }_t{\alpha }_{j,q}+{\alpha }_{j,q}^{\ast }i{\partial }_t{\alpha }_{j,c}-{\Delta }_j\left({\alpha }_{j,c}^{\ast }{\alpha }_{j,q}+{\alpha }_{j,q}^{\ast }{\alpha }_{j,c}\right)\right)\hfill \\ \hfill & -J\left({\alpha }_{1,q}^{\ast }{\alpha }_{2,c}+{\alpha }_{1,c}^{\ast }{\alpha }_{2,q}+{\alpha }_{2,q}^{\ast }{\alpha }_{1,c}+{\alpha }_{2,c}^{\ast }{\alpha }_{1,q}\right)\hfill \\ \hfill & -\frac{1}{2}\sum _{j=1,2}{U}_j\left(|{\alpha }_{j,c}{|}^2+{\left|{\alpha }_{j,q}\right|}^2\right)\left({\alpha }_{j,c}^{\ast }{\alpha }_{j,q}+{\alpha }_{j,c}{\alpha }_{j,q}^{\ast }\right)\hfill \\ \hfill & -\sqrt{2}F\left(t\right)\sum _{j=1,2}\left({\alpha }_{j,q}^{\ast }+{\alpha }_{j,q}\right)+i\sum _{j=1,2}{\gamma }_j\left(2|{\alpha }_{j,q}{|}^2+{\alpha }_{j,c}{\alpha }_{j,q}^{\ast }-{\alpha }_{j,c}^{\ast }{\alpha }_{j,q}\right)].\hfill \end{array}$$

Here, ${\Delta }_j$ is the detuning frequency, $U_j$ is the Kerr nonlinearity, J is the coupling rate, and $F\left(t\right)$ is an external drive. Varying $S_K$ with respect to ${\alpha }_{j,c}^{\ast }$ and ${\alpha }_{j,q}^{\ast }$ yields the semiclassical equations of motion that incorporate both Hamiltonian and dissipative dynamics.

Appendix C.1. Mean-Field Equations and the Effective Potential

Mean-field or semiclassical (or saddle-point) dynamics are found by setting ${\displaystyle \frac{\delta S_K}{\delta {\alpha }_{j,c}}}={\displaystyle \frac{\delta S_K}{\delta {\alpha }_{j,q}}}=0$ and similarly for complex conjugates. One obtains ${\alpha }_{j,q}=0$ as a trivial solution [24], and the mean-field dynamics for ${\alpha }_{j,c}$ then follow from ${\displaystyle \frac{\delta S_K}{\delta {\alpha }_{j,q}^{\ast }}}{|}_{{\alpha }_{j,q}=0}=0$ giving(A11)$${\partial }_t{\alpha }_{j,c}=-\left({\gamma }_j+i{\Delta }_j\right){\alpha }_{j,c}-iJ{\alpha }_{j^{\prime },c}-\frac{i}{2}{U}_j{\alpha }_{j,c}{\left|{\alpha }_{j,c}\right|}^2-i\sqrt{2}F\left(t\right),$$ where $j^{\prime }\ne j$. This can be recast as a potential dynamics of the form$$i{\partial }_t{\alpha }_{j,c}={\partial }_{{\alpha }_{j,c}^{\ast }}V({\overrightarrow{\alpha }}_c,{\overrightarrow{\alpha }}_c^{\ast })-i{\gamma }_j{\alpha }_{j,c}+\sqrt{2}F\left(t\right),$$ with the effective potential$$V({\overrightarrow{\alpha }}_c,{\overrightarrow{\alpha }}_c^{\ast })=\sum _{j=1,2}\left({\Delta }_j|{\alpha }_{j,c}{|}^2+\frac{1}{4}{U}_j{\left|{\alpha }_{j,c}\right|}^4\right)+J\left({\alpha }_{1,c}{\alpha }_{2,c}^{\ast }+{\alpha }_{2,c}{\alpha }_{1,c}^{\ast }\right),$$ which is the potential landscape discussed in Section 3.2.

Appendix C.2. Fluctuations and the Inverse Green’s Function

To analyze the system’s linear response to small perturbation about a mean-field or saddle-point solution, we expand the Keldysh action to the second order in the fluctuating fields $\delta {\alpha }_{j,c/q}$. Namely, we set$${\alpha }_{j,c/q}\left(t\right)\to {\alpha }_{j,c/q}\left(t\right)+\delta {\alpha }_{j,c/q}\left(t\right),$$ expand $S_K$ in (A10) to the quadratic order in $\delta {\alpha }_{j,c/q}$, and then consider the Fourier decomposition of the fluctuations$$\delta {\alpha }_{c/q}\left(t\right)={\displaystyle \frac{1}{\sqrt{2}}}\int {\displaystyle \frac{d\omega }{2\pi }}\left[\delta {\alpha }_{c/q}\left(\omega \right)e^{-i\omega t}+\delta {\alpha }_{c/q}(-\omega )e^{+i\omega t}\right].$$

In block-matrix form, the resulting quadratic action reads(A12)$$S_{K,2}=\frac{1}{2}\int {\displaystyle \frac{d\omega }{2\pi }}\delta {\Phi }^{†}\left(\omega \right)\left(\begin{array}{cc}0& {\left[G^A\left(\omega \right)\right]}^{-1}\\ {\left[G^R\left(\omega \right)\right]}^{-1}& D^K\end{array}\right)\delta \Phi \left(\omega \right),$$ where $\delta \Phi \left(\omega \right)$ groups the fluctuations $\{\delta {\alpha }_{j,c/q}\left(\omega \right),\delta {\alpha }_{j,c/q}^{\ast }(-\omega )\}$. Explicitly, in the vector component such that $$\delta \Phi \left(\omega \right)={\left(\delta {\alpha }_{1,c}\left(\omega \right),\delta {\alpha }_{1,c}^{\ast }(-\omega ),\delta {\alpha }_{2,c}\left(\omega \right),\delta {\alpha }_{2,c}^{\ast }(-\omega ),\delta {\alpha }_{1,q}\left(\omega \right),\delta {\alpha }_{1,q}^{\ast }(-\omega ),\delta {\alpha }_{2,q}\left(\omega \right),\delta {\alpha }_{2,q}^{\ast }(-\omega )\right)}^T,$$ the inverse of retarded/advanced Green’s function ${\left[G^{R/A}\left(\omega \right)\right]}^{-1}$ and the Keldysh component of the inverse Green’s function $D^K$ read, respectively,(A13)$$\begin{array}{cc}\hfill {\left[G^R\left(\omega \right)\right]}^{-1}& =\left(\begin{array}{cccc}\omega -g_1+i{\gamma }_1& -{\displaystyle \frac{1}{2}}{U}_1{\alpha }_{1,c}^2& -J& 0\\ -{\displaystyle \frac{1}{2}}{U}_1{\left({\alpha }_{1,c}^{\ast }\right)}^2& -\omega -g_1-i{\gamma }_1& 0& -J\\ -J& 0& \omega -g_2+i{\gamma }_2& -{\displaystyle \frac{1}{2}}{U}_2{\alpha }_{2,c}^2\\ 0& -J& -{\displaystyle \frac{1}{2}}{U}_2{\left({\alpha }_{2,c}^{\ast }\right)}^2& -\omega -g_2-i{\gamma }_2\end{array}\right),\hfill \end{array}$$(A14)$$\begin{array}{cc}\hfill G^A{\left(\omega \right)]}^{-1}& ={\left[G^R{\left(\omega \right)}^{†}\right]}^{-1},\mathrm{and}\hfill \end{array}$$(A15)$$\begin{array}{cc}\hfill D^K& =2i\mathrm{diag}({\gamma }_1,{\gamma }_1,{\gamma }_2,{\gamma }_2),\hfill \end{array}$$ where $g_j={\Delta }_j+U_j{\left|{\alpha }_{j,c}\right|}^2$ for $j=1,2$.

Appendix C.3. Classical Field Response

We are interested in the fluctuations of the classical field variables $\delta {\alpha }_{j,c}\left(\omega \right)$, which is precisely encoded the inverse retarded block, ${\left[G^R\left(\omega \right)\right]}^{-1}$ of (A12), that determines how those fluctuations grow or decay. Let $\delta F\left(t\right)$ be a small external drive coupling linearly to the oscillator modes as in (A11). At the level of fluctuations, one can write$$\delta {\alpha }_{j,c}\left(\omega \right)=G_j^R\left(\omega \right)\delta F\left(\omega \right),$$ where $G_j^R\left(\omega \right)$ is the relevant component (or linear combination of components) of the retarded Green’s function. In the time domain,$$\delta {\alpha }_{j,c}\left(t\right)={\int }_{-\infty }^{\infty }d{t}^{\prime }{G}_j^R(t-t^{\prime })\delta F\left(t^{\prime }\right),$$ with $G_j^R\left(t\right)=0$ for $t<0$. This causality condition sets the retarded (not advanced) nature of $G^R$; the system cannot respond before the perturbation arises.

Appendix C.4. Pole Structure and Oscillation vs. Decay

To analyze the dynamics near a stationary solution, e.g., ${\alpha }_{1,c}={\alpha }_{2,c}=0$, we linearize around that point and compute $G^R\left(\omega \right)$ by inverting the block matrix. There, the inverse retarded Green’s function of (A12) becomes(A16)$${\left[G^R\left(\omega \right)\right]}^{-1}=\left(\begin{array}{cccc}\omega -{\Delta }_1+i{\gamma }_1& 0& -J& 0\\ 0& -\omega -{\Delta }_1-i{\gamma }_1& 0& -J\\ -J& 0& \omega -{\Delta }_2+i{\gamma }_2& 0\\ 0& -J& 0& -\omega -{\Delta }_2-i{\gamma }_2\end{array}\right).$$

The poles of $G^R\left(\omega \right)$ appear where ${\left[G^R\left(\omega \right)\right]}^{-1}\to 0,$ i.e., where the determinant of the retarded block vanishes. Writing such a pole as ${\omega }_{★}=\Omega \pm i\Gamma$ clarifies the physical nature of the fluctuation mode:

• $\Omega =\mathrm{Re}\left({\omega }_{★}\right)$ is the frequency of oscillation.

• $\Gamma =-\mathrm{Im}\left({\omega }_{★}\right)$ is the exponential decay rate if $\Gamma >0$, signifying a stable, dissipative mode.

By setting ${\Delta }_1={\Delta }_2=\Delta$ and ${\gamma }_1={\gamma }_2=\gamma >0$ as discussed in Section 3.2, one obtains the poles describing the fast and slow relaxation modes in (12), that is$${\omega }_s=\pm \left|J-|\Delta |\right|-i\gamma ,{\omega }_f=\pm \left|J+|\Delta |\right|-i\gamma .$$

The real part of a slow-mode pole starts to disappear when the effective potential becomes marginally flat (e.g., $J\simeq |\Delta |$). Although one might expect a zero-frequency “soft” oscillation in a conservative dynamics setting, dissipation (encoded in $-\mathrm{i}\gamma$ terms of ${\left[G^R\left(\omega \right)\right]}^{-1}$) shifts that would-be neutral oscillatory mode into an overdamped decay channel.

In summary, the Keldysh action formalism provides a powerful lens to derive both the mean-field equations of motion and the fluctuation response in an open quantum system. Retarded Green’s functions, in particular, capture the causality of how a driven perturbation modifies the system at later times, thereby revealing the presence of overdamped or oscillatory collective modes. These results underpin the discussion in Section 3.2 on how the disappearance of slow modes leads to synergistic encoding at $J\approx |\Delta |$.

Appendix D. Second-Order Cumulant Expansion

In this appendix, we derive the second-order (2nd-order) cumulant expansion used to obtain semiclassical equations of motion for eight complex-valued expectation variables:$$\left\{⟨{\widehat{a}}_1⟩,⟨{\widehat{a}}_2⟩,⟨{\widehat{n}}_1⟩,⟨{\widehat{n}}_2⟩,⟨{\widehat{a}}_1^2⟩,⟨{\widehat{a}}_2^2⟩,⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩,⟨{\widehat{a}}_1{\widehat{a}}_2⟩\right\}.$$

This second-order cumulant expansion serves as an interpolation between simpler mean-field dynamics, where second-order cumulants factorize and yield no correlations, and the full quantum description, in which higher-order cumulants can be nonzero when quantum correlations are sufficiently strong. We follow the standard cumulant-truncation scheme [66], imposing that all third- and fourth-order cumulants vanish [67]. That is, we set$${⟨\widehat{A}\widehat{B}\widehat{C}⟩}_C=0\mathrm{and}{⟨\widehat{A}\widehat{B}\widehat{C}\widehat{D}⟩}_C=0,$$ where the subscript “C” denotes the connected (cumulant) part. This approximation can capture second-order correlations while keeping the system of equations tractable.

Below, we provide the resulting coupled ordinary differential equations (ODEs). These govern the dynamics of our reduced set of expectation values under the second-order expansion:(A17)$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}⟨{\widehat{a}}_j⟩=-& \left({\gamma }_j+i{\Delta }_j\right)⟨{\widehat{a}}_j⟩-i{U}_j\left(⟨{\widehat{a}}_j^{†}⟩⟨{\widehat{a}}_j^2⟩+2⟨{\widehat{a}}_j⟩⟨{\widehat{n}}_j⟩-2⟨{\widehat{a}}_j^{†}⟩{⟨{\widehat{a}}_j⟩}^2\right)\hfill \\ \hfill & -iJ\left({\delta }_{1j}⟨{\widehat{a}}_2⟩+{\delta }_{2j}⟨{\widehat{a}}_1⟩\right)-iF\left(t\right)\hfill \end{array}$$(A18)$${\displaystyle \frac{d}{dt}}⟨{\widehat{n}}_j⟩=-2{\gamma }_j⟨{\widehat{n}}_j⟩-iJ\left({\delta }_{1j}-{\delta }_{2j}\right)\left(⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩-⟨{\widehat{a}}_2^{†}{\widehat{a}}_1⟩\right)+iF\left(t\right)\left(⟨{\widehat{a}}_j⟩-⟨{\widehat{a}}_j^{†}⟩\right)$$(A19)$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}⟨{\widehat{a}}_j^2⟩=-& 2\left({\gamma }_j+i{\Delta }_j\right)⟨{\widehat{a}}_j^2⟩-2iJ⟨{\widehat{a}}_1{\widehat{a}}_2⟩-2iF\left(t\right)⟨{\widehat{a}}_j⟩\hfill \\ \hfill & -i{U}_j\left(⟨{\widehat{a}}_j^2⟩+6⟨{\widehat{n}}_j⟩⟨{\widehat{a}}_j^2⟩-4⟨{\widehat{a}}_j^{†}⟩{⟨{\widehat{a}}_j⟩}^3\right),\hfill \end{array}$$(A20)$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩=& i\left({\Delta }_1-{\Delta }_2\right)⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩-\left({\gamma }_1+{\gamma }_2\right)⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩\hfill \\ \hfill & -iJ\left(⟨{\widehat{n}}_1⟩-⟨{\widehat{n}}_2⟩\right)+iF\left(t\right)\left(⟨{\widehat{a}}_2⟩-⟨{\widehat{a}}_1^{†}⟩\right)\hfill \\ \hfill & +i{U}_1\left(2⟨{\widehat{n}}_1⟩⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩+⟨{\widehat{a}}_1^{†2}⟩⟨{\widehat{a}}_1{\widehat{a}}_2⟩-2⟨{\widehat{a}}_1^{†2}⟩⟨{\widehat{a}}_1⟩⟨{\widehat{a}}_2⟩\right)\hfill \\ \hfill & -i{U}_2\left(2⟨{\widehat{n}}_2⟩⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩+⟨{\widehat{a}}_2^2⟩⟨{\widehat{a}}_1^{†}{\widehat{a}}_2^{†}⟩-2⟨{\widehat{a}}_1^{†}⟩⟨{\widehat{a}}_2^{†}⟩{⟨{\widehat{a}}_2⟩}^2\right),\hfill \end{array}$$(A21)$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}⟨{\widehat{a}}_1{\widehat{a}}_2⟩=& -\sum _{j=1,2}\left({\gamma }_j+i{\Delta }_j\right)⟨{\widehat{a}}_1{\widehat{a}}_2⟩-i\sum _{j=1,2}\left(J⟨{\widehat{a}}_j^2⟩+F\left(t\right)⟨{\widehat{a}}_j⟩\right)\hfill \\ \hfill & -i{U}_1\left(2⟨{\widehat{n}}_1⟩⟨{\widehat{a}}_1{\widehat{a}}_2⟩+⟨{\widehat{a}}_1^{†}{\widehat{a}}_2⟩⟨{\widehat{a}}_1^2⟩-2⟨{\widehat{a}}_1^{†}⟩⟨{\widehat{a}}_1^2⟩⟨{\widehat{a}}_2⟩\right)\hfill \\ \hfill & -i{U}_2\left(2⟨{\widehat{n}}_2⟩⟨{\widehat{a}}_1{\widehat{a}}_2⟩+⟨{\widehat{a}}_2^{†}{\widehat{a}}_1⟩⟨{\widehat{a}}_2^2⟩-2⟨{\widehat{a}}_2^{†}⟩⟨{\widehat{a}}_1⟩⟨{\widehat{a}}_2^2⟩\right),\hfill \end{array}$$ where ${\delta }_{ij}$ is the Kronecker delta. In the main text, we assume ${\Delta }_1={\Delta }_2$ and ${\gamma }_1={\gamma }_2$ for simplicity.

Figure A3. Classical mutual information for the 2nd-order cumulant description. This is to be contrasted with Figure 2 (left) to reveal how second-order description captures partial but nontrivial correlation effects correcting mean-field approximation. We set [Forumla omitted. See PDF.] to represent a dynamical regime with non-negligible correlations, motivating the use of second-order cumulants description.

Note that these second-order equations provide an approximation to go beyond a strict mean-field approximation without the full computational cost of higher-order correlation. By discarding third- and fourth-order cumulants, we retain the information that captures pairwise correlations, which often dominate many relevant dynamics, while avoiding an intractable explosion in the number of degrees of freedom.

References

1. Jaeger, H. Short Term Memory in Echo State Networks; GMD Report Fraunhofer-Gesellschaft: München Germany, 2001.

2. Lukoševičius, M.; Jaeger, H. Reservoir computing approaches to recurrent neural network training. Comput. Sci. Rev.; 2009; 3, pp. 127-149. [DOI: https://dx.doi.org/10.1016/j.cosrev.2009.03.005]

3. Nakajima, K. Reservoir computing: Theory, physical implementations, and applications. IEICE Tech. Rep.; 2018; 118, pp. 149-154.

4. Lukoševičius, M. A practical guide to applying echo state networks. Neural Networks: Tricks of the Trade; 2nd ed. Springer: Berlin/Heidelberg, Germany, 2012; pp. 659-686.

5. Fujii, K.; Nakajima, K. Harnessing disordered-ensemble quantum dynamics for machine learning. Phys. Rev. Appl.; 2017; 8, 024030. [DOI: https://dx.doi.org/10.1103/PhysRevApplied.8.024030]

6. Innocenti, L.; Lorenzo, S.; Palmisano, I.; Ferraro, A.; Paternostro, M.; Palma, G.M. Potential and limitations of quantum extreme learning machines. Commun. Phys.; 2023; 6, 118. [DOI: https://dx.doi.org/10.1038/s42005-023-01233-w]

7. Motamedi, A.; Zadeh-Haghighi, H.; Simon, C. Correlations Between Quantumness and Learning Performance in Reservoir Computing with a Single Oscillator. arXiv; 2023; arXiv: 2304.03462

8. Götting, N.; Lohof, F.; Gies, C. Exploring quantumness in quantum reservoir computing. Phys. Rev. A; 2023; 108, 052427. [DOI: https://dx.doi.org/10.1103/PhysRevA.108.052427]

9. Martínez-Peña, R.; Ortega, J.P. Quantum reservoir computing in finite dimensions. Phys. Rev. E; 2023; 107, 035306. [DOI: https://dx.doi.org/10.1103/PhysRevE.107.035306] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37072987]

10. Sannia, A.; Martínez-Peña, R.; Soriano, M.C.; Giorgi, G.L.; Zambrini, R. Dissipation as a resource for Quantum Reservoir Computing. Quantum; 2024; 8, 1291. [DOI: https://dx.doi.org/10.22331/q-2024-03-20-1291]

11. Xiong, W.; Facelli, G.; Sahebi, M.; Agnel, O.; Chotibut, T.; Thanasilp, S.; Holmes, Z. On fundamental aspects of quantum extreme learning machines. Quantum Mach. Intell.; 2025; in press

12. Mujal, P.; Martínez-Peña, R.; Giorgi, G.L.; Soriano, M.C.; Zambrini, R. Time-series quantum reservoir computing with weak and projective measurements. NPJ Quantum Inf.; 2023; 9, 16. [DOI: https://dx.doi.org/10.1038/s41534-023-00682-z]

13. Govia, L.; Ribeill, G.; Rowlands, G.; Krovi, H.; Ohki, T. Quantum reservoir computing with a single nonlinear oscillator. Phys. Rev. Res.; 2021; 3, 013077. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.013077]

14. Ivaki, M.N.; Lazarides, A.; Ala-Nissila, T. Quantum reservoir computing on random regular graphs. arXiv; 2024; arXiv: 2409.03665

15. Sornsaeng, A.; Dangniam, N.; Chotibut, T. Quantum next generation reservoir computing: An efficient quantum algorithm for forecasting quantum dynamics. Quantum Mach. Intell.; 2024; 6, 57. [DOI: https://dx.doi.org/10.1007/s42484-024-00188-7]

16. Martínez-Peña, R.; Giorgi, G.L.; Nokkala, J.; Soriano, M.C.; Zambrini, R. Dynamical phase transitions in quantum reservoir computing. Phys. Rev. Lett.; 2021; 127, 100502. [DOI: https://dx.doi.org/10.1103/PhysRevLett.127.100502]

17. Mujal, P.; Martínez-Peña, R.; Nokkala, J.; García-Beni, J.; Giorgi, G.L.; Soriano, M.C.; Zambrini, R. Opportunities in Quantum Reservoir Computing and Extreme Learning Machines. Adv. Quantum Technol.; 2021; 4, 2100027. [DOI: https://dx.doi.org/10.1002/qute.202100027]

18. Hartmann, M.; Brandão, F.; Plenio, M. Quantum many-body phenomena in coupled cavity arrays. Laser Photonics Rev.; 2008; 2, pp. 527-556. [DOI: https://dx.doi.org/10.1002/lpor.200810046]

19. Gat, I.; Tishby, N. Synergy and redundancy among brain cells of behaving monkeys. Proceedings of the Advances in Neural Information Processing Systems; Denver, CO, USA, 30 November–5 December 1998; Volume 11.

20. Brenner, N.; Strong, S.P.; Koberle, R.; Bialek, W.; Steveninck, R.R.d.R.v. Synergy in a neural code. Neural Comput.; 2000; 12, pp. 1531-1552. [DOI: https://dx.doi.org/10.1162/089976600300015259] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10935917]

21. Schneidman, E.; Bialek, W.; Berry, M.J. Synergy, redundancy, and independence in population codes. J. Neurosci.; 2003; 23, pp. 11539-11553. [DOI: https://dx.doi.org/10.1523/JNEUROSCI.23-37-11539.2003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14684857]

22. Latham, P.E.; Nirenberg, S. Synergy, redundancy, and independence in population codes, revisited. J. Neurosci.; 2005; 25, pp. 5195-5206. [DOI: https://dx.doi.org/10.1523/JNEUROSCI.5319-04.2005]

23. Bertschinger, N.; Rauh, J.; Olbrich, E.; Jost, J.; Ay, N. Quantifying unique information. Entropy; 2014; 16, pp. 2161-2183. [DOI: https://dx.doi.org/10.3390/e16042161]

24. Kamenev, A. Field Theory of Non-Equilibrium Systems; Cambridge University Press: Cambridge, UK, 2023.

25. Dudas, J.; Carles, B.; Plouet, E.; Mizrahi, F.A.; Grollier, J.; Marković, D. Quantum reservoir computing implementation on coherently coupled quantum oscillators. NPJ Quantum Inf.; 2023; 9, 64. [DOI: https://dx.doi.org/10.1038/s41534-023-00734-4]

26. Dudas, J.; Grollier, J.; Marković, D. Coherently coupled quantum oscillators for quantum reservoir computing. Proceedings of the 2022 IEEE 22nd International Conference on Nanotechnology (NANO); Palma de Mallorca, Spain, 4–8 July 2022; pp. 397-400. [DOI: https://dx.doi.org/10.1109/NANO54668.2022.9928763]

27. Kalfus, W.D.; Ribeill, G.J.; Rowlands, G.E.; Krovi, H.K.; Ohki, T.A.; Govia, L.C.G. Hilbert space as a computational resource in reservoir computing. Phys. Rev. Res.; 2022; 4, 033007. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.4.033007]

28. Khan, S.A.; Hu, F.; Angelatos, G.; Türeci, H.E. Physical reservoir computing using finitely-sampled quantum systems. arXiv; 2021; arXiv: 2110.13849

29. Milburn, G.; Holmes, C. Quantum coherence and classical chaos in a pulsed parametric oscillator with a Kerr nonlinearity. Phys. Rev. A; 1991; 44, 4704. [DOI: https://dx.doi.org/10.1103/PhysRevA.44.4704] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9906513]

30. Milburn, G.; Corney, J.; Wright, E.M.; Walls, D. Quantum dynamics of an atomic Bose-Einstein condensate in a double-well potential. Phys. Rev. A; 1997; 55, 4318. [DOI: https://dx.doi.org/10.1103/PhysRevA.55.4318]

31. Carusotto, I.; Ciuti, C. Quantum fluids of light. Rev. Mod. Phys.; 2013; 85, 299. [DOI: https://dx.doi.org/10.1103/RevModPhys.85.299]

32. Goto, H.; Kanao, T. Chaos in coupled Kerr-nonlinear parametric oscillators. Phys. Rev. Res.; 2021; 3, 043196. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.043196]

33. Breuer, H.P.; Petruccione, F. The Theory of Open Quantum Systems; Oxford University Press on Demand: Oxford, UK, 2002.

34. Bena, I. Dichotomous Markov noise: Exact results for out-of-equilibrium systems. Int. J. Mod. Phys. B; 2006; 20, pp. 2825-2888. [DOI: https://dx.doi.org/10.1142/S0217979206034881]

35. Abel, B.; Marquardt, F. Decoherence by quantum telegraph noise: A numerical evaluation. Phys. Rev. B; 2008; 78, 201302. [DOI: https://dx.doi.org/10.1103/PhysRevB.78.201302]

36. Franco, R.L.; D’Arrigo, A.; Falci, G.; Compagno, G.; Paladino, E. Entanglement dynamics in superconducting qubits affected by local bistable impurities. Phys. Scr.; 2012; 2012, 014019. [DOI: https://dx.doi.org/10.1088/0031-8949/2012/T147/014019]

37. Makkeh, A.; Theis, D.O.; Vicente, R. BROJA-2PID: A robust estimator for bivariate partial information decomposition. Entropy; 2018; 20, 271. [DOI: https://dx.doi.org/10.3390/e20040271] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33265362]

38. Soriente, M.; Chitra, R.; Zilberberg, O. Distinguishing phases using the dynamical response of driven-dissipative light-matter systems. Phys. Rev. A; 2020; 101, 023823. [DOI: https://dx.doi.org/10.1103/PhysRevA.101.023823]

39. Soriente, M.; Heugel, T.L.; Omiya, K.; Chitra, R.; Zilberberg, O. Distinctive class of dissipation-induced phase transitions and their universal characteristics. Phys. Rev. Res.; 2021; 3, 023100. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.023100]

40. Alaeian, H.; Soriente, M.; Najafi, K.; Yelin, S. Noise-resilient phase transitions and limit-cycles in coupled kerr oscillators. arXiv; 2021; arXiv: 2106.04045[DOI: https://dx.doi.org/10.1088/1367-2630/ad2414]

41. Dambre, J.; Verstraeten, D.; Schrauwen, B.; Massar, S. Information Processing Capacity of Dynamical Systems. Sci. Rep.; 2012; 2, 514. [DOI: https://dx.doi.org/10.1038/srep00514]

42. Martínez-Peña, R.; Ortega, J.P. Input-dependence in quantum reservoir computing. arXiv; 2024; arXiv: 2412.08322

43. Verstraete, F.; Wolf, M.M.; Ignacio Cirac, J. Quantum computation and quantum-state engineering driven by dissipation. Nat. Phys.; 2009; 5, pp. 633-636. [DOI: https://dx.doi.org/10.1038/nphys1342]

44. Fry, D.; Deshmukh, A.; Chen, S.Y.C.; Rastunkov, V.; Markov, V. Optimizing quantum noise-induced reservoir computing for nonlinear and chaotic time series prediction. Sci. Rep.; 2023; 13, 19326. [DOI: https://dx.doi.org/10.1038/s41598-023-45015-4]

45. Leghtas, Z.; Touzard, S.; Pop, I.M.; Kou, A.; Vlastakis, B.; Petrenko, A.; Sliwa, K.M.; Narla, A.; Shankar, S.; Hatridge, M.J. et al. Confining the state of light to a quantum manifold by engineered two-photon loss. Science; 2015; 347, pp. 853-857. [DOI: https://dx.doi.org/10.1126/science.aaa2085] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25700514]

46. Shannon, C.E. A mathematical theory of communication. Bell Syst. Tech. J.; 1948; 27, pp. 379-423. [DOI: https://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x]

47. Bell, A.J. The co-information lattice. Proceedings of the Fifth International Workshop on Independent Component Analysis and Blind Signal Separation: ICA; Granada, Spain, 22–24 September 2003; Volume 2003.

48. James, R.G.; Ellison, C.J.; Crutchfield, J.P. Anatomy of a bit: Information in a time series observation. Chaos Interdiscip. J. Nonlinear Sci.; 2011; 21, 037109. [DOI: https://dx.doi.org/10.1063/1.3637494] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21974672]

49. Harder, M.; Salge, C.; Polani, D. Bivariate measure of redundant information. Phys. Rev. E; 2013; 87, 012130. [DOI: https://dx.doi.org/10.1103/PhysRevE.87.012130] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23410306]

50. Olbrich, E.; Bertschinger, N.; Rauh, J. Information decomposition and synergy. Entropy; 2015; 17, pp. 3501-3517. [DOI: https://dx.doi.org/10.3390/e17053501]

51. Rosas, F.E.; Mediano, P.A.; Gastpar, M.; Jensen, H.J. Quantifying high-order interdependencies via multivariate extensions of the mutual information. Phys. Rev. E; 2019; 100, 032305. [DOI: https://dx.doi.org/10.1103/PhysRevE.100.032305]

52. Gutknecht, A.J.; Makkeh, A.; Wibral, M. From Babel to Boole: The Logical Organization of Information Decompositions. arXiv; 2023; arXiv: 2306.00734

53. Williams, P.L.; Beer, R.D. Nonnegative decomposition of multivariate information. arXiv; 2010; arXiv: 1004.2515

54. Bertschinger, N.; Rauh, J.; Olbrich, E.; Jost, J. Shared information—New insights and problems in decomposing information in complex systems. Proceedings of the European Conference on Complex Systems 2012; Brussels, Belgium, 3–7 September 2012; Springer: Cham, Switzerland, 2013; pp. 251-269.

55. Banerjee, P.K.; Griffith, V. Synergy, redundancy and common information. arXiv; 2015; arXiv: 1509.03706

56. Wibral, M.; Priesemann, V.; Kay, J.W.; Lizier, J.T.; Phillips, W.A. Partial information decomposition as a unified approach to the specification of neural goal functions. Brain Cogn.; 2017; 112, pp. 25-38. [DOI: https://dx.doi.org/10.1016/j.bandc.2015.09.004]

57. Makkeh, A.; Gutknecht, A.J.; Wibral, M. Introducing a differentiable measure of pointwise shared information. Phys. Rev. E; 2021; 103, 032149. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.032149] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33862718]

58. Schick-Poland, K.; Makkeh, A.; Gutknecht, A.J.; Wollstadt, P.; Sturm, A.; Wibral, M. A partial information decomposition for discrete and continuous variables. arXiv; 2021; arXiv: 2106.12393

59. Rauh, J.; Banerjee, P.K.; Olbrich, E.; Jost, J.; Bertschinger, N. On extractable shared information. Entropy; 2017; 19, 328. [DOI: https://dx.doi.org/10.3390/e19070328]

60. Griffith, V.; Koch, C. Quantifying synergistic mutual information. Guided Self-Organization: Inception; Springer: Berlin/Heidelberg, Germany, 2014; pp. 159-190.

61. McGill, W. Multivariate information transmission. Trans. IRE Prof. Group Inf. Theory; 1954; 4, pp. 93-111. [DOI: https://dx.doi.org/10.1109/TIT.1954.1057469]

62. Yeung, R.W. A new outlook on Shannon’s information measures. IEEE Trans. Inf. Theory; 1991; 37, pp. 466-474. [DOI: https://dx.doi.org/10.1109/18.79902]

63. Schneidman, E.; Still, S.; Berry, M.J.; Bialek, W. Network information and connected correlations. Phys. Rev. Lett.; 2003; 91, 238701. [DOI: https://dx.doi.org/10.1103/PhysRevLett.91.238701] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14683220]

64. Dalla Torre, E.G.; Diehl, S.; Lukin, M.D.; Sachdev, S.; Strack, P. Keldysh approach for nonequilibrium phase transitions in quantum optics: Beyond the Dicke model in optical cavities. Phys. Rev. A; 2013; 87, 023831. [DOI: https://dx.doi.org/10.1103/PhysRevA.87.023831]

65. Sieberer, L.M.; Buchhold, M.; Diehl, S. Keldysh field theory for driven open quantum systems. Rep. Prog. Phys.; 2016; 79, 096001. [DOI: https://dx.doi.org/10.1088/0034-4885/79/9/096001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27482736]

66. Kubo, R. Generalized cumulant expansion method. J. Phys. Soc. Jpn.; 1962; 17, pp. 1100-1120. [DOI: https://dx.doi.org/10.1143/JPSJ.17.1100]

67. Sánchez-Barquilla, M.; Silva, R.; Feist, J. Cumulant expansion for the treatment of light–matter interactions in arbitrary material structures. J. Chem. Phys.; 2020; 152, 034108. [DOI: https://dx.doi.org/10.1063/1.5138937] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31968946]