Neuromorphic Functions of Light in Parity‐Time‐Symmetric Systems

Abstract As an elementary processor of neural networks, a neuron performs exotic dynamic functions, such as bifurcation, repetitive firing, and oscillation quenching. To achieve ultrafast neuromorphic signal processing, the realization of photonic equivalents to neuronal dynamic functions has attracted considerable attention. However, despite the nonconservative nature of neurons due to energy exchange between intra‐ and extra‐cellular regions through ion channels, the critical role of non‐Hermitian physics in the photonic analogy of a neuron has been neglected. Here, a neuromorphic non‐Hermitian photonic system ruled by parity‐time symmetry is presented. For a photonic platform that induces the competition between saturable gain and loss channels, dynamical phases are classified with respect to parity‐time symmetry and stability. In each phase, unique oscillation quenching functions and nonreciprocal oscillations of light fields are revealed as photonic equivalents of neuronal dynamic functions. The proposed photonic system for neuronal functionalities will become a fundamental building block for light‐based neural signal processing.

Interacting nonlinear channels allow complex dynamical functions in open systems.For example, the competition between sodium and potassium channels with opposite gradients of ionic concentrations constructs a biological building block 1 for neural dynamics, including the Hopf bifurcation 14 with the 'all-or-none' characteristic and the repetitive firing 3 that enables the synchronization of neural networks 15 .Nonlinear resistance of memristors 2 has also been intensively applied to the electrical realization of chaotic systems 16 for encryption and random number generation.In these dynamical phenomena, the state-dependent nonlinear operation of each channel (e.g.potential-dependent gate dynamics 1 ) and the interaction between channels (e.g.Kirchhoff's law for entire ionic currents 1 ) determine the detailed form of the resulting dynamical function (e.g.threshold and repetitive firing in Hodgkin-Huxley neurons 1 ).
Nonlinear channels for electromagnetic waves also lead to critical phenomena in nonequilibrium regimes, especially from saturable dynamics of amplification or absorption.In classical laser systems, the nonlinear response of the gain saturation is one of the key parameters determining the laser throughput 17 , and the saturable absorber embedded in a laser cavity allows the mode locking and Qswitching for pulsed emission 7,8 .Recently, inspired by the concept of parity-time (PT) symmetry 18 in non-Hermitian quantum mechanics 19 , the investigation on saturated optical channels has been revisited in terms of PT-symmetric phase transition 20,21 .The alteration of the PT-symmetric phase from the gain saturation has been explored to derive the interesting phenomena of non-reciprocal transparency of light 9,10 and robust wireless power transfer in radio-frequency circuits 11 .The evolution around the PTsymmetric exceptional point, also known as dynamical encircling, is now receiving growing attention for non-reciprocal mode conversion 13,22 with its topological robustness to nonlinear channels 12 .Notably, in contrast to competing ionic channels in Hodgkin-Huxley neurons 1 , only a single nonlinear channel (saturable gain or loss) has been considered in the field of PT symmetry, missing possible important forms of dynamical functions from the interaction between channels.
In this paper, we investigate the competition between saturable gain and loss channels in PTsymmetric systems, and explore the dynamical phase diagram defined by equilibrium and oscillation parameters.We reveal that three different phases exist: fast-and slow-response equilibrium phases, and a nonequilibrium phase with weak oscillations.At the phase boundary between fast-response equilibrium and oscillatory nonequilibrium phases, the existence of the non-reciprocal strong oscillation state is demonstrated, which is a missing dynamical function in single nonlinear channel systems [9][10][11][12][13]22 . Utiizing the notion of the limit cycle in the real-valued parameter space 3,14 , the chaotic nature of this phase boundary is analyzed for the applications of the intensity-based bifurcation and fully tunable time delay of repetitive resonator excitations.Due to their neuromorphic behaviors of threshold and repetitive firing, merged with PT-symmetric chaotic properties, the proposed PT-symmetric systems with competing saturable channels may become a building block for light-based artificial neural networks and encrypting systems.
The schematic of the dynamical PT-symmetric system is illustrated in Fig. 1a.The gain and loss resonators with field amplitudes ψG,L are coupled each other with the coupling coefficient κ, while gain and loss channels possess intensity-dependent nonlinear saturations each with the gain 4 and loss 5 coefficients of γG,L = γG,L(|ψG,L|) = γG0,L0/(1+|ψG,L/ψGs,Ls| 2 ), where γG0,L0 are the linear coefficients, and ψGs,Ls are the saturation intensities.While the gain saturation generally occurs in pumped gain media (e.g.erbium-doped film 4,9 ), the saturated loss can be obtained by utilizing organic dyes 23 , carbon nanotube 24 , graphene 25 , and artificial realizations 26 .The amplification and absorption inside the resonators then decrease with increased field intensity (Fig. 1b), dynamically affecting the channels of the PT-symmetric system.From the temporal coupled mode theory 6 , the governing equation of the PTsymmetric system is then written as where ω0 is the resonant frequency of each isolated resonator.
The nonlinear temporal Eq. ( 1) can be solved numerically by applying the 6 th -order Runge-Kutta method 27 (time step = 2π/(200•ω0)) with the initial condition of ψG,L(t=0).To analyze dynamical behaviors of the system from the solutions ψG,L(t), we define the characteristic parameters α and β, each measuring the equilibrium and oscillation condition of the system during the time range T1 ≤ t ≤ T2.In detail, for the field intensity of the entire system I(t) = |ψG(t)| 2 + |ψL(t)| 2 , the 1 st order polynomial fitting of log[I(t)] ~ α•t + δ presents the averaged equilibrium condition during T1 ≤ t ≤ T2: α = 0 for the equilibrium phase, α > 0 for amplifying, and α < 0 for attenuating nonequilibrium phases with exponential evolutions.As well, inspired by the activity oscillation amplitude in neural networks 28 , the temporal intensity oscillation inside the PT-symmetric system during T1 ≤ t ≤ T2 is then quantified by β, which measures the temporal deviation from the average signal activity e α•t+δ as where a high value of  implies a high amplitude of the oscillations of the average activity inside the system 28 .
Figures 1c and 1d each represent the calculated dynamical phase diagram of α and β, as functions of the relative strength (γL0/γG0) and saturation level (ψLs/ψGs) of the channels, with the initial condition of ψG(t=0) = 1 and ψL(t=0) = 0. From the dynamical characteristic parameters α and β, three distinct phases of nonequilibrium (N), fast-response-(Ef) and slow-response-(Es) equilibriums are then observed with different equilibrium and oscillation conditions, also accompanying phase boundaries and a triple point.When the loss channel approaches the saturation level faster than the gain channel (ψLs < ψGs), the nonequilibrium phase N with the field amplification is obtained (Fig. 1e), while faster saturation of the gain channel (ψGs < ψLs) leads to the equilibrium phases Ef and Es.The equilibrium phases can then be divided into two distinct phases of Ef with fast temporal response (Fig. 1f) and Es with slow temporal response (Fig. 1g) according to the relative strength of the gain and loss channels (γL0/γG0).It is also noted that a single nonlinear channel system [9][10][11][12][13] (ψLs→∞) corresponds to a single 'line' in those two-dimensional phase diagrams.
In contrast to continuous phase transition in linear PT-symmetric systems 20,21 , the transition between different dynamical phases occurs discontinuously across phase boundaries.For example, considering that the quasi-static phase of PT symmetry (see Supplementary Note S1) is unbroken in N and Es regimes, and is broken in the Ef regime, exotic 'phase boundaries' in terms of the oscillation are observed between unbroken and broken PT-symmetric phases of N-Ef or Ef-Es (red filled symbols in Fig. 1c,d).As shown in those boundaries in Fig. 1d, strong dynamical oscillation (with large β) inside the system can be achieved at the boundary between unbroken and broken PT-symmetric phases, having distinct eigenmodal profiles and real-complex eigenvalues 20 .The field-amplitude-dependent variation of the gain and loss coefficients, for different saturation levels (ψ Gs ,ψ Ls ).
Dynamical phase diagrams defined by (c) the equilibrium parameter α and (d) the oscillation parameter β: N for the nonequilibrium phase, and E f (or E s ) for the equilibrium phase of the fast (or slow) temporal response.Yellow filled circles in (c,d) present the cases in (e-g), for the temporally varying field intensities inside the gain (ψ G ) and loss (ψ L ) resonators.Red filled circles in (c,d) present the case of the oscillating phase boundaries N-E f and E f -E s , while a red empty circle denotes the non-oscillating phase boundary N-E s .The initial condition is ψ G (t = 0) = 1 and ψ L (t = 0) = 0 for (c-g).The time range for calculating α and and ψ Gs = 2 for all cases.
Among two phase boundaries N-Ef and Ef-Es with strong oscillations (or large β), we now focus on the phase boundary N-Ef, which offers the exotic and unique dynamical function of non-reciprocal strong oscillations, absent in the phase boundary Ef-Es, the triple point N-Ef-Es, and the single channel system [9][10][11][12][13]22 (see Supplementary Note S2-S4). We irst analyze the behavior of the phase boundary N-Ef, in terms of the initial excitation channel.Figure 2a-d represents the enlarged α-β phase diagrams around the phase boundary N-Ef for the initial excitation conditions of ψG(t=0) = 1, ψL(t=0) = 0 (Fig. 2a,b) and ψG(t=0) = 0, ψL(t=0) = 1 (Fig. 2c,d).Due to the discontinuity at the boundary N-Ef, the switching of 'phase' with strong dependence on the initial excitation condition is observed in Fig. 2a-d, for both α (from N to Ef phase) and β (from non-oscillating to oscillating state) parameters (red dashed lines).In detail, for the gain-resonator excitation (Fig. 2a,b), the overall system remains at the N phase with a low value of β, leading to the weakly-oscillating amplifying response (Fig. 2e).In contrast, for the lossresonator excitation (Fig. 2c,d), the system stays at the equilibrium phase Ef with a high value of β.This large β with α = 0 corresponds to the repetitive excitation of both resonators (Fig. 2f) with the oscillation across the exceptional point (see Supplementary Note S5 for oscillating PT-symmetric transition).
This 'chaotic' phenomenon with a highly sensitive response to the initial excitation condition originates from the non-reciprocity of the dynamical PT-symmetric system, furthermore, resembling a repetitive firing response of Hodgkin-Huxley neurons 1,3 : the analogy of sodium (incoming current) and potassium (outgoing current) gated channels, with gain (amplified light) and loss (dissipative light) saturated channels.To further look into the chaotic, oscillatory behavior at the phase boundary N-Ef, we now develop the limit cycle analysis which has been utilized for analyzing repetitive firing in Hodgkin-Huxley neurons 3,14 .For the field amplitude ψG,L, we employ the representation based on real-valued parameters 29 I(t), θ(t), ϕ(t), and ω(t), as With some elaborations (see Supplementary Note S6), Eq. ( 1) becomes the closed equation of the form of The evolutions of optical states are shown in Fig. 3a,b, each for the gain and loss resonator excitation.It is noted that the directional rotation in the θ-ϕ parameter space occurs for both cases, which corresponds to the alternating transition between amplitude (θ) and phase (φ) differences of the field in each resonator, in accordance with the oscillating PT-symmetric phase transitions (Supplementary Note S5).However, depending on the initial excitation channel (red and blue symbols in the I(t=0) = 1, θ-ϕ green plane of Fig. 3), the bifurcation of the system between the nonequilibrium state (Fig. 3a) and the equilibrium state with the limit cycle operation (Fig. 3b) is evident, representing the 'open' and 'closed' cycles for each case in the θ-I plane (bottom figures in Fig. 3).This non-reciprocity originates from the distinct initial variation of gain-loss distributions determined by the excited channel, leading to the chaotic state between weakly-oscillating nonequilibrium and strongly-oscillating equilibrium phases (Fig. 2e,3a versus 2f,3b).Black dashed lines show the evolution of the excitation peak, controlled by the excited field intensity.Only the loss resonator is excited at t = 0.All other parameters are the same as those in Fig. 1.
In summary, we have shown that the competing phenomenon between saturable gain and loss channels is described in the form of the dynamical phase diagram.Each phase exhibits different equilibrium and oscillation features, in relation to the quasi-static variation of PT-symmetric gauges and phases.At the abrupt phase boundary between these dynamical phases, the chaotic state with nonreciprocity exists with the realization of non-reciprocal repetitive firing in resonators.From the chaotic feature, the convergence speed to the limit cycle can also be manipulated, allowing the fully tunable time delay of the resonator excitation.Considering the practical realization from the controlled saturation with engineered energy bands 30 or artificial realizations 26 , the fully tunable repetitive firing of electromagnetic waves may become a proper building block for neuromorphic wave circuits, underpinning the critical role of the competition between nonlinear dynamics for high-level functions or network synchronizations.
For the comprehensible understanding of the dynamical PT-symmetric system, in this note, we assume the adiabatic evolution of the system, which enables the quasi-static analysis.For Eq. ( 1) in the main text, the instantaneous eigenvalues ω e (t) of the system then becomes (S2) The temporal function γ avg (t) determines the gauge of the PT-symmetric system 1,2 , as γ avg > 0 for the active gauge, γ avg < 0 for the passive gauge, and γ avg = 0 for the exact PT symmetry.Meanwhile, D(t) defines the phase of PT symmetry 3 : D(t) > 0 for unbroken phase, D(t) < 0 for broken phase, and D(t) = 0 for exceptional point.
Figure S1 shows the temporal evolution of gauges γ avg (t) (Fig. S1a-c) and PT-symmetric phases D(t) (Fig. S1d-f) in the regimes of N, E f , and E s phases of the dynamical PT-symmetric system, with the initial condition of ψ G (t=0) = 1 and ψ L (t=0) = 0. Firstly, N and E s phases which have the strong oscillations of γ avg (t) and D(t) maintain the unbroken PT-symmetric phase (Fig. S1d,f), due to the power exchange from two distinct eigenvalues with different eigenmodes in the unbroken PT symmetry 1,3 .In contrast, the phase of E f corresponds to the broken PT-symmetric phase with the identical real part of eigenvalues 1,3 (Fig. S1e), leading to the rapid convergence to the steady state.
From the gauge analysis, it is also shown that the nonequilibrium property of the N phase originates from the spiking of the gauge γ avg (t) of the system (enlarged plot of Fig. S1a).However, as shown in the similarity of the gauge and PT-symmetric phase between N and E s phases (Fig. S1a,d versus Fig. S1c,f), it is noted that the dynamical phase diagrams of α and β provide deeper understanding of the temporal dynamics with the clear distinction of the nonequilibrium (N) and equilibrium (E s ) phases.

Note S2. Quasi-reciprocal phase boundary E f -E s
The phase boundary N-E f in between the nonequilibrium (N) and equilibrium (E f ) phases is highly sensitive to the initial condition, leading to the chaotic property and the following non-reciprocity, intensity-based bifurcation and time delay (see Fig. 2-4 in the main text).In contrast, at the phase boundary E f -E s in between equilibrium (E f and E s ) phases, the overall energy inside the system is converged, with the significantly suppressed chaotic nature.
Figure S2a-d represents the enlarged phase diagrams of α and β around the phase boundary E f -E s , for the initial condition of ψ G (t=0) = 1, ψ L (t=0) = 0 (Fig. S2a,b) and ψ G (t=0) = 0, ψ L (t=0) = 1 (Fig. S2c,d).Because of the equilibrium condition in E f and E s phases, the initial condition does not affect much the nonlinear gain and loss channels, leading to quasi-reciprocal responses (Fig. S2a,b versus Fig. S2c,d).However, around the phase boundary E f -E s which is near the exceptional point of linear PT-symmetric systems (D(t) = 0 when γ G0 = γ L0 , for κ = γ G0 and ψ Gs , Ls → ∞, Eq. (S2)), the equilibrium condition of α varies drastically (Fig. S2a,c).Figure S3 shows the temporal evolution at each point near the phase boundary E f -E s .While preserving quasi-reciprocal responses, the convergences to non-oscillating constant states (Fig. S3a,d) and oscillating states (Fig. S3b,c) are observed with relatively slower response than that at the phase boundary N-E f (Fig. 2e,f in the main text).

Note S4. Comparison with a single-nonlinear-channel system
In the two-dimensional phase diagrams from competing saturable channels (Fig. S5a), a single nonlinear channel system [4][5][6][7][8] corresponds to the single 'line'.We compare the regimes including the phase boundaries of N-E f (Fig. S5b) and E f -E s (Fig. S5c) with the single channel system (Fig. S5d), which has the gain saturation and linear loss (ψ Ls →∞, red dashed lines in Fig. S5a).In Fig. S5b-d

Note S5. Quasi-static analysis of the phase boundary N-E f
The results of the quasi-static analysis for the phase boundary N-E f are shown in Fig. S6, in terms of gauges γ avg (t) (Fig. S6a,b) and PT-symmetric phases D(t) (Fig. S6c,d).When the gain resonator is excited (ψ G (t=0) = 1 and ψ L (t=0) = 0, Fig. S6a,c), there exists the oscillation between active and passive gauges (Fig. S6a) which results in the nonequilibrium condition as shown in the main text.It is also noted that because the PT-symmetric phase is at the unbroken regime (Fig. S6c, D(t) > 0), the evolution of optical energy has the form of the directional coupling 3 (Fig. 2e in the main text).
In contrast, for the case of the loss-resonator excitation, while the passive gauge is maintained (Fig. S6b), the oscillation of PT-symmetric phase occurs around the exceptional point (Fig. S6d, D

Note S6. Derivation of the limit cycle equation
From ψ G (t) = [I(t)]

Figure 1 .
Figure 1.Dynamical phase diagram of the PT-symmetric system with competing saturable channels.(a) A schematic for the system with saturable channels, from the nonlinear coefficient functions γ G , L = γ G , L (|ψ G , L |).(b)

Figure 2 .
Figure 2. Non-reciprocal emergence of repetitive resonator firing of the PT-symmetric system.The enlarged phase diagrams around the phase boundary N-E f in the cases of (a,b) gain-resonator excitation and (c,d) lossresonator excitation, for (a,c) the equilibrium parameter α and (b,d) the oscillation parameter β. (e,f) The temporally varying field intensities in the gain (ψ G ) and loss (ψ L ) resonators for each case of the yellow filled circles in Fig. 2a,b and Fig. 2c,d.All other parameters are the same as those in Fig. 1.

Figure 3 .
Figure 3. Non-reciprocity of the dynamical evolutions of optical states in the real-valued parameter space.Temporal variation of optical states in the I-θ-ϕ space for the (a) gain-resonator excitation and (b) loss-resonator

Figure S1 .
Figure S1.Quasi-static analysis of N, E f , and E s phases.(a-c) Temporal gauges γ avg (t) and (d-f) PT symmetry phase parameter D(t) for (a,d) N, (b,e) E f , and (c,f) E s phases.The dynamical gain and loss coefficients are calculated from the 6 th -order Runge-Kutta analysis, for the calculation of Eq. (S2).All other parameters are the same as in Fig. 1 in the main text.

Figure S2 .
Figure S2.Quasi-reciprocal response at the phase boundary E f -E s .The enlarged phase diagrams around the phase boundary E f -E s in the cases of (a,b) gain-resonator excitation and (c,d) loss-resonator excitation, for (a,c) the equilibrium parameter α and (b,d) the oscillation parameter β.Each yellow filled circle denotes the distinguished case in the equilibrium phase diagram α, with the temporal evolutions in Fig. S3.All other parameters are the same as those in Fig. 1.

Figure S3 . 1 .
Figure S3.Temporal dynamics near the phase boundary E f -E s .(a-d) The temporally varying field intensities in the gain (ψ G ) and loss (ψ L ) resonators, for each case of the yellow filled circles in Fig. S2, with different values of the relative strength of the loss channel γ L0 /γ G0 : (left) gain-resonator excitation and (right) loss-resonator excitation.All other parameters are the same as those in Fig. 1.

Figure S4 .
Figure S4.Phase diagrams near the triple point.The enlarged phase diagrams around the triple point N-E f -E s in the cases of (a,b) gain-resonator excitation and (c,d) loss-resonator excitation, for (a,c) the equilibrium parameter α and (b,d) the oscillation parameter β.All other parameters are the same as in Fig. 1 in the main text.
, the uniqueness of the non-reciprocity (difference between red and blue lines) and strong oscillation (large β) is evident only at the phase boundary N-E f , when compared with the quasi-reciprocal and weakly oscillating responses at the phase boundary E f -E s and in the single channel case.

Figure S5 .
Figure S5.Comparison with a single nonlinear channel system.(a) Phase diagrams of α and β, representing the plotted regimes in (b-d).Black dashed lines denote the regimes including the phase boundaries of N-E f and E f -E s .Red dashed line denotes the single nonlinear channel system.(b-d) The evolutions of α and β, as a function of the relative strength (γ L0 /γ G0 ) of the channels for the regimes across (b) N-E f , (c) E f -E s , and (d) for the single channel case.Red lines represent the cases of the gain resonator excitation, and blue lines represent the cases of the loss resonator excitation.All other parameters are the same as in Fig. 1 in the main text.
(t) ~ 0) which directly determines the form of the repetitive resonator excitation (Fig. 2f in the main text).For example, because of the evolution near the exceptional point, two resonators are almost synchronized from the coalescence of eigenmodes 3 .

Figure S6 .
Figure S6.Quasi-static analysis of the phase boundary N-E f .(a,b) Temporal gauges γ avg (t) and (c,d) PT symmetry phase parameter D(t) for the (a,c) gain-resonator and (b,d) loss-resonator excitations.The saturated gain and loss coefficients are calculated from the 6 th -order Runge-Kutta analysis, for the calculation of Eq. (S2).All other parameters are the same as in Fig. 1 in the main text.
the first and second equations in Eq. (S5) derives the first and second equations of Eq. (3) in the main text, while the combination of the third and fourth equations of Eq. (S5) derives the third equation of Eq. (3) in the main text with dω/dt = ω 0 + κ•cos(ϕ)/sin(2θ).