Oscillator networks with N‐shaped nonlinearities: Electrical modeling and wave digital emulation

This paper proposes contributions to the efficient wave digital (WD) modeling of large oscillator networks which are emerging as energy‐efficient alternatives to traditional computers. The WD concept enables in‐operando parameter tuning, real‐time testing, and the associated algorithms are highly parallelizable. We present a general electrical model of N‐shaped nonlinearities that are commonly found in nonlinear oscillators. Our model offers the flexibility to design the current–voltage characteristic based on specific requirements. We show how this model can be used to derive efficient and explicit WD algorithms for nonlinear oscillators. Furthermore, we propose the use of lossless transmission lines between the oscillators and the coupling network to obtain an ideal circuit for an oscillator network that can function as an Ising machine and be efficiently and exactly evaluated in the WD domain. The proposed algorithms are compared against the classical method involving iterative techniques, and their capabilities are evaluated through the emulation of a single FitzHugh‐Nagumo oscillator as well as an Ising machine involving transmission lines. In the latter case, we show that, for large networks, the proposed methods decrease the runtime by up to 75% compared to using iterative techniques.


| INTRODUCTION
Computational paradigms based on oscillator networks have been gaining popularity as a promising energy-efficient supplementation to the classical von Neumann computer. 1,2These networks are usually large and consist of many oscillators performing computations in parallel. 3One prominent example are Oscillator-based Ising Machines (OIMs).Ising machines are systems that are governed by the Ising Hamiltonian and evolve towards its ground state.They are of interest because all NP-problems can be mapped to them. 4,5,3Here, NP stands for Non-deterministic Polynomial-time, referring to a set of computational problems for which there is no known algorithm to compute the solution in polynomial-time.Polynomial-time means that the number of steps that a Turing machine requires to solve a problem can be expressed as a polynomial of the input size. 6This makes NP-problems difficult to solve for conventional digital computers.As OIMs are inherently massively parallel, they may offer an advantage for finding solutions for such problems. 4We discuss them in more detail in Section 4.2 and also refer to Reference 4 for a great overview.
Emulating such networks, however, requires great computational effort.A suitable tool are wave digital (WD) algorithms 7 which are robust, massively parallel, energetically consistent and allow for in-operando parameter tuning and real-time testing.They involve mapping a circuit port-wisely onto a so-called WD model, preserving the original topology.This enables efficient modeling and emulation of circuits.Especially repetitive circuit structures like oscillator networks can be efficiently programmed and evaluated with vector-valued components.Emulating circuits with nonlinear components using the WD concept usually leads to implicit relationships.These so-called delay-free loops must be resolved by iterative methods. 8,9However, closer physical modeling using transmission lines to impose finite signal propagation speed can help avoid iterations and thus computational load.
At the heart of OIMs lay self-sustained oscillators.To generate sustained oscillations, many oscillators rely on locally active components.These destabilize the equilibrium point of the system while semi-passivity guarantees bounded trajectories, leading to oscillations. 10,11For example, the relaxation oscillator, 12 the Hindmarsh-Rose oscillator, 13 the FitzHugh-Nagumo oscillator (FNO), 14 or the classical differential LC oscillator 15 the locally active component takes the form of an N-shaped nonlinearity.The name of these one-ports stems from their i, u ð Þ curve which assumes the shape of the letter "N".Its outer strokes correspond to dissipative regions with positive resistance that are separated by the locally active region with negative differential resistance.The first negative-resistance solid state device was the tunnel diode, 16 used in the circuit of the FNO 14 which is the oscillator considered in this paper.An algorithm for the systematic generation of passive devices with negative differential resistance, consisting of linear resistors and two transistors, was proposed in Reference 16.In Reference 17, operational amplifiers are used to generate a piecewise linear N-shaped nonlinearity for a circuit implementation of the FNO.Based on the nonlinear function synthesis with transconductance-amplifiers from References 18,19, generates the piecewise-linear nonlinearity in a CMOS implementation of the FNO.A variety of memristor-based nonlinear oscillators is discussed in Reference 20, where (locally) active components consist of a memristor and a negative resistance, implemented with a standard op-amp circuit.In Reference 21, an N-shaped nonlinearity is modeled using hyperbolic tangent memristors.An N-shaped locally active memristor is investigated in Reference 22 to construct oscillators.
In this work, we present a general electrical model for N-shaped nonlinearities that aims to be efficient and easily parametrizable to simulate systems with a broad variety of N-shaped nonlinearities, regardless of their implementation.Concretely, based on technical requirements, the i,u ð Þ curve of our model can be actively designed and we can derive from it an explicit WD model, 7 similar to References 23,24.We then make use of the proposed electrical model to derive efficient and explicit WD algorithms for nonlinear oscillators that we couple dissipatively to obtain oscillator networks.To deepen the modeling and take finite propagation speeds into account, we introduce transmission lines between the oscillators and the coupling network.We obtain an oscillator network that can be parametrized to be functionally equivalent to a network modeled with iterative techniques while reducing the computational load for evaluation.Most notably, the resulting WD model of the oscillator network is explicit, that is, no delay-free loops occur.
The remainder of this paper is divided as follows.In Section 2, we present our general electrical model for N-shaped nonlinearities.In Section 3, we make use of this model to derive an explicit WD model.This model is then used in Section 4 for the modeling of oscillator networks, OIMs specifically, employing transmission lines.In Section 5, we present and compare emulation results and, finally, in Section 6 we summarize our main contributions.

| MODELING N-SHAPED NONLINEARITIES
In this section, first, we briefly recapitulate the Lambert W function, which is necessary for the second part in which we propose a general electrical model for N-shaped nonlinearities.

| Lambert W function
5][26] It is defined implicitly as Þ¼x and with it the solution of the equation can be written as 25 : x It cannot be written analytically but be computed iteratively as is readily available in many toolboxes.

| General electrical model
While N-shaped nonlinearities are commonly realized by interconnections of transistors, 11 we propose to model their i, u ð Þ curve by the simple circuit depicted in Figure 1A.Inspired by, 23,24 we combine an active clipper diode model with a negative resistance, where the latter element can, for example, be implemented by the negative impedance converter in Figure 1B.The i, u ð Þ curve resulting from this interconnection is depicted in Figure 1C.The four sources j 0 , e 0 , e 1 and e 2 as well as the resistors R 1 , R 2 and G 0 are design parameters to shift and shape the nonlinearity.Writing out all current and voltage relations gives, where we have modeled the diode characteristics by the Shockley equation.The diodes are denoted by D ν where ν ¼ 1, 2 such that i ν , u ν and R ν are the current through, voltage across and resistance in series with the diode, respectively.Furthermore, by convention, for AE and ∓ , ν ¼ 1 selects the upper and ν ¼ 2 the lower sign throughout the paper.The N-shape in the i,u ð Þ curve arises from the exponentially growing diode currents with a negative sign for D 1 , a positive sign for D 2 and the linearly decreasing current i 3 through the negative differential resistance.The superposition of these currents leads to the curve depicted in Figure 1C.Next, we discuss the influence of the design parameters.From (2) can be seen that e 0 shifts the i, u ð Þ curve horizontally and j 0 shifts it vertically.Further analysis is impeded by the diode currents being defined implicitly: Inserting (2d) into (2b) reveals a dependence of i ν on itself in the exponential of the Shockley equation.Using the Lambert W function, this implicit relation can be resolved to F I G U R E 1 (A) General equivalent circuit of an N-shaped nonlinearity with G 0 > 0. (B) Practical implementation of a negative conductance using a negative impedance converter.(C) i,u ð Þ curve corresponding to the circuit in (A) where all sources are zero and R ν !0.
The extrema are marked as i m,ν ,u m,ν ð Þ .
Plugging (2b) to (2d) and ( 3) into (2a) yields an explicit formulation, To determine the location of the extrema, we make two assumptions: First, the reverse currents of the diodes are constant such that their derivatives vanish.Second, the two diodes are never simultaneously in forward mode as that does not result in an N-shaped curve.Solving di=du ¼ 0 for u, the Lambert W function conveniently vanishes in the expression for the locations of the extrema: Here, ν corresponds to the conductive diode D ν .This result shows that the sources e ν can individually shift the extrema along the u-axis.The conductance G 0 only appears in i 3 .It allows adjusting the negative differential resistance and thus also the currents i m,ν at the extrema.Lastly, the resistors R ν appear in the exponential of (2b) and decrease the slope of the i,u ð Þ characteristic of the diodes.They thus control the slope of the left and right strokes of the N-shape.The effects of all parameters are visualized in Figure 2. Together, they allow fine-grained control when designing a nonlinearity to specific requirements.Moreover, as some of the parameters are sources, the curve can be dynamically translated in the i,u ð Þ plane and the locations of the maxima can be placed freely at any time.

| EXPLICIT WAVE DIGITAL MODEL
Given the general electrical model, we now derive an explicit wave digital model of the N-shaped nonlinearities.To that end, we first briefly introduce the wave digital concept.

| Wave digital concept
When emulating electrical circuits, it is advantageous to have a modular digital representation.However, attempting to obtain such a representation by discretizing circuit components and connecting them to a signal flow diagram (SFD) leads, even for linear elements, to delay-free loops.Consequently, evaluating the SFD is not possible.While monolithic models, like differential equations, sacrifice modularity to solve this issue, the Wave Digital (WD) concept was originally introduced to emulate linear filter networks in a fully modular fashion. 7Classical WD models of linear circuits Influence of the parameters e 0 , j 0 , R ν , G 0 , and e ν on the shape of the nonlinearity.The black curve is identical to Figure 1C and results from all sources being zero and R ν !0. The green dashed curves stem from varying the stated parameters.keep advantageous properties like passivity and losslessness from the reference circuit and are furthermore very robust to finite word lengths of coefficients. 27,28One drawback is that traditional WD models are restricted to a tree-like structure and can only accommodate a single one-port nonlinearity at the root. 8However, more recent advances expand the WD concept towards a more general circuit emulation technique.One approach to allow for multiple nonlinearities and more general topologies is using a fixed-point iteration 8 : Here, delay-free loops are broken by introducing fractional delays.At every timestep of the simulation, the states of dynamic elements are held constant while iterating the fractional delay, converging to a fixed point, and solving the implicit relations.Due to its straightforward implementation, this is the technique to which we compare our explicit WD model of an oscillator network.There are, however, other approaches: An alternative to fixed-point iterations is the use of the Newton-Raphson method in the WD domain, introduced in Reference 29 to handle circuits with multiple nonlinearities.Here quadric convergence speeds can be attained.Circuits with reactive elements and multiple nonlinearities can also be efficiently emulated using the Scattering Iterative Method (SIM) when discretizing the circuit elements with Linear Multi-Step Discretization Methods (LMSDMs). 30Vectorization of memoryless linear two-ports can enable connecting them to topological junctions without delay-free loops. 31Today, the WD concept is a powerful tool for circuit emulation that we consider because of the mentioned favorable passivity, losslessness and robustness properties, the modularity of the resulting models as well as efficient computability.Furthermore, using WD algorithms, software emulations of large oscillator networks are easily implementable.
Deriving a WD model involves, first, decomposing a circuit into a set of one-and multiports, including Kirchhoff connections.Second, their dynamic equations are discretized.In this work, we use the trapezoidal rule which preserves passivity properties.However, other discretization methods are applicable to, for example, gain accuracy or computation speed. 30Third, a bijective variable transform, maps from the i, u ð Þ domain to the wave digital domain.Here, the new variables a and b correspond to voltage wave quantities and denote the incident and reflected wave of a port, respectively.The constant R is the port resistance, which, in theory, can have an arbitrary positive value.However, an appropriate choice of it can simplify the resulting WD structures 7 and indeed be necessary for the realizability of the wave flow diagram which is, fourth, obtained by reconnecting the WD elements according to their original topology.In the remainder of the paper, we will introduce further WD components and related concepts as they are needed.For a more in-depth introduction to the topic and a list of translated components, see Reference 7.

| Explicit wave digital N-shaped nonlinearity
The goal is now to derive a scattering function b ¼ S N a ð Þ for the voltage-controlled N-shaped nonlinearity that is equivalent to (4).This scattering relation will describe the entire one-port from Figure 1 monolithically, without needing to translate each of its components to the WD domain.This scattering function should be explicit in the WD domain since when evaluating a WD model, voltages and currents are unknown.Furthermore, for an explicit scattering relation, no delay-free loop arises when connecting it to a reflection-free port.Preliminarily, we check under which conditions such a function exists. 32 If g is strictly monotonic, g À1 will exist.From the definition of strict monotony, boundaries for R can be derived for which a scattering function exists: Only one of the conditions can hold at a time since a function cannot be strictly monotonically increasing and decreasing at once.For (4), the evaluation of the boundaries yields À inf u di=du ¼ G 0 and Àsup u di=du ¼ À∞.As R is a positive constant, RG 0 < 1 must hold for the WD function to exist.We note that in this work, we seek an exact mapping of the nonlinearity to the WD domain.Alternatively, Bernardini and Sarti 33 start with a piecewise-linear approximation of f in the i,u ð Þ domain to construct an explicit WD function.When modeling a Shockley diode, they report significant speedups over an exact representation using the Lambert W function.To derive the explicit WD function, we start by bringing the implicit i, u ð Þ curve (2) into the WD domain using the variable transform (6).Both a and b appear in the exponentials of the diode currents, rendering the function implicit.This can be resolved when making an assumption: At least one diode D ν is always non-conductive with i ν ≈ 0. To decide which diode that is, we use the middle point between the extrema of the i, u ð Þ curve as a decision boundary: With this distinction, we can write the still implicit a, b ð Þ curve: where, Like the diode current before, this can be written explicitly with help of the Lambert W function by bringing (10a) into the form of (1a) such that, where the appearing constants are given by: With (1b) and x ¼ b, the explicit, nonlinear scattering relation can be written as, with different parameters for u larger or smaller than u 0 , depending on which diode D ν is active.However, as discussed, the voltage u is unknown during the evaluation of the WD model.Since due to (8), RG 0 < 1 must hold, we know that d du a > 0 which means that we can transform the decision boundary u 0 into a boundary a 0 to make the same decision.
Given that for u ¼ u 0 both diodes are non-conductive, we can write i 0 ¼ i N u 0 ð Þ¼ÀG 0 e 2 À e 1 ½ =2 þ j 0 .With (6), we get, where the sign indicates the decision: For a < a 0 we thus set ν ¼ 1 and otherwise ν ¼ 2.

| OSCILLATOR NETWORKS
Given the explicit WD model of the N-shaped nonlinearity, in this section, we introduce the Fitzhugh-Namugo oscillator (FNO) which we emulate to demonstrate the correctness of the WD model.Furthermore, we consider a network of dissipatively coupled oscillators with transmission lines capable of solving combinatorial optimization problems, functioning as an Ising machine.

| FitzHugh-Nagumo oscillator
The circuit implementing the FNO is depicted in Figure 3A, marked by the gray dashed box.It is based on the electrical implementation of FitzHugh's simplification of the Hodgkin-Huxley model of the squid giant axon 34 by Nagumo et al. 14 The active N-shaped nonlinearity i N u ð Þ enables self-sustaining oscillations and is implemented as presented in Section 2 and shown in Figure 1A.The electrical parameters of the oscillator circuit are given in Table 1.Like the nonlinearity before, the FNO circuit will be translated into the WD domain to allow for efficient emulation.The resulting WD model is depicted in Figure 4, indicated by the gray dashed box.As per, 7 it results from port-wisely decomposing the circuit in Figure 3A, discretizing the elements, translating them into the WD domain via (6), and interconnecting them with the same topology using WD adaptors that implement the Kirchhoff laws.The capacitor translates to a delay with the port resistance R C ¼ T= 2C ½ , where T denotes the sampling period of the WD algorithm.The inductor translates to a delay element with sign inversion and its port resistance is R L ¼ 2L=T.The resistor R e translates to a zero source and a sink, while, as per Section 3.2, the N-shaped nonlinearity, i N u ð Þ, translates to S N from (12).Lastly, the series (parallel) connection translates to a series (parallel) adaptor, marked by a line with a dot (two parallel lines), implementing the corresponding Kirchhoff laws.Both adaptor types can have one reflection-free port, indicated by a ⊥ symbol.At that port, the adaptor's reflected wave is independent of the incident wave.This enables connecting components whose reflected waves instantaneously depend on their incident waves without utilizing computationally costly iterative techniques. 8We use it to connect the discussed N-shaped nonlinearity, obtaining an explicit WD model, that is, a model without delay-free loops.To demonstrate the correctness of the WD model, we emulate a single FNO and compare the resulting output voltage to an LTspice simulation, where the diodes and negative resistance are replaced by models of real components.Here, we use the op-amp (TL071) and the diode model (1N4148).When fitted correctly, we see in T A B L E 1 Emulation parameters used in Figures 3 and 6.

Electrical parameters
Figure 3B that the WD model produces the same output as the LTspice model.When emulating n FNOs, we make use of a vector-valued WD model as depicted in Figure 4 where all waves are stacked in vectors such that efficient matrix computations can be used to evaluate the wave flow diagram in parallel.

| Oscillator-based Ising machine
Having established the single oscillator, we now consider a set of coupled FNOs, forming an oscillator network, as depicted in Figure 3C.More specifically, we consider an oscillator-based Ising machine.The foundation for the Ising machine and the source of its name is the Ising model. 35Originally, it was developed for understanding ferromagnetism where it describes the energy of magnetically coupled spins with nearest-neighbor interactions.Mathematically, the Ising model describes a system, where binary random variables s μ À1, þ1 f g are the nodes in a weighted graph.The edge weights describe their interaction strengths. 36The Ising Hamiltonian describes the total energy of this system and can be written as, where J ¼ J T is an adjacency matrix containing the weights J μν between spins s μ and s ν .Finding the ground state of the Ising Hamiltonian is called the Ising problem.Importantly, all NP-problems can be mapped to it. 37,5,3A well-known problem of the sort is the max-cut problem.Solving it requires to partition the nodes of a graph into two sets such that the sum of the cut (weighted) edges is maximized.Mathematically, it can be expressed as where W is a weighted adjacency matrix.For J ¼ ÀW , the Ising and the max-cut problem are solved by the same argument.Due to this direct mapping, max-cut problems are a prime example problem for Ising machines.Ising machines are systems that are governed by the Ising Hamiltonian such that over time they move towards the ground state of the Ising Hamiltonian.Coupled oscillators with phase-encoded binary variables are Ising machines: When the oscillators are coupled according to some graph, the max-cut on that graph can after some time be read from the relative phases of the oscillators. 38The spins can be phase encoded by injecting a subharmonic injection locking signal (SHIL) with double the oscillator's natural frequency Ω 0 , causing the phases to asymptotically lock to a relative value of 0 or π. 39 We implement this by a current source added in parallel to the capacitor in Figure 3A, which, when choosing e 0 ¼ 0V, can be subsumed under the source j 0 of the N-shaped nonlinearity circuit in Figure 1A.The SHIL current can be written as, The coupling topology is given by the adjacency matrix J and is equivalently described by the incidence matrix N obtained by assigning an arbitrary direction to the undirected edges.Every edge is electrically implemented by an elementary coupling circuit as depicted in the light green dashed box in Figure 3C.As we will only consider max-cut problems on graphs with positive weights, the Ising machine only needs couplings implementing negative weights.Negatively coupling oscillators can be implemented by inverting the output voltage of one oscillator and resistively coupling it to the non-inverted output of another. 40As, for the sake of simplicity, we will only consider a unit weight problem, all coupling conductances are chosen uniformly to be G μν ¼ R À1 c where a fitting value is determined experimentally and given in Table 1.For an arbitrarily directed edge, let the voltage of the FNO at the start of the edge be denoted u s and the voltage of the FNO at the end be u e .Nodal analysis yields that the evoked current at the elementary coupling circuit implementing the edge e ϰ ¼ u s,ϰ , u e,ϰ ð Þcan be written as . Note that the current only depends on the summation of the voltages and not on them individually, making the coupling dissipative.For vector-valued notation, the incidence matrix can be used to write the vectors u s and u e containing the u s,ϰ and u e,ϰ of all edges: With N j j denoting the elementwise absolute of N, define Since N s þ N e ¼ N j j, the vector of elementary coupling currents can be written as, where G is a diagonal matrix containing the coupling conductances G ϰ .The considered problems graphs are undirected and we imposed an arbitrary direction in which the elementary coupling currents flow.Because in reality the graphs are bidirectional, the elementary coupling currents also contribute to the coupling current at the oscillator where an edge starts, such that i ¼ N e þ N s ½ j.The coupling can then in total be written as, where G c is a conductivity matrix that appears in the green dashed box in Figure 3 to draw the vector-valued circuit of the Ising machine.Now, we derive the WD description of the coupling circuit by applying the bijective transformation ( 6) to (19).We obtain the scattering relation, with the port resistance matrix In the WD model of the OIM in Figure 4, this scattering relation is marked by the dashed green box.As discussed above, the reflection-free port of the parallel adaptor of the FNO is occupied with the N-shaped nonlinearity.Connecting the coupling network to another port of the parallel adaptor results in a delay-free topological loop where incident and reflected waves depend on each other instantaneously.This loop can be resolved with a fixed-point iteration, which, however, impairs simulation time. 8Instead, we propose a lossless transmission line with minimal delay T between the output port of the FNO and the coupling structure, as marked by the dashed blue boxes in Figures 3A,C and 4.This does two things: First it deepens the modeling as it takes the finite propagation speed of information into account.][43] In the WD domain, a transmission line equals a delay in the path of both the incident and reflected wave as depicted in Figure 4, marked by the blue dashed box.That is equal to a 4-port circulator with two capacitors with capacitance C, dictating the port resistance R T ¼ T= 2C ½ . 7This means that when employing the transmission line, the port resistance R T can no longer be arbitrarily chosen as it reflects a component in the electrical circuit.This is illustrated in Figure 5: For values of R T that are too low or too high, the signal shapes and amplitudes degenerate.
We aim thus to find a value for R T that will keep the signal shapes and amplitudes close to an oscillator network that employs iterative techniques to break the delay-free loop while stressing that there cannot be strict equality between the two.Recall that a transmission line with matched termination transforms the load impedance to the input.Perfectly matched termination would thus yield a model of an oscillator network without a transmission line and without delay-free loops.We are, however, restricted to R T ¼ 1R T .Consequently, we aim to find a heuristic for R T that approximates matched termination as best as possible, minimizing the power reflected by the transmission line.To achieve that, we choose the wave resistance R T as the solution of where Á k k 2 denotes the spectral norm.If G c is diagonalizable it can be written as, where P is a matrix containing the eigenvectors of G c and , where λ G c is a vector containing the eigenvalues λ G c ,μ of G c .With this, the above problem can now equivalently be rewritten as, holds such that we find a compact heuristic when instead solving which yields where n ¼ 1 T 1 is the number of nodes in the graph.

| SOLVING AN ISING PROBLEM
Now, we will emulate the oscillator network depicted in Figure 3A,C by evaluating its wave digital model from Figure 4.The coupling topology is given by the adjacency matrix J ¼ ÀW , where, for the sake of simplicity, the entries of W are 0 or 1.This configuration corresponds to an Ising machine which solves the max-cut problem on W .As there is only one weight value, the coupling resistance for every edge can be chosen to be G ϰ ¼ 1=R c .The problem that we consider features n ¼ 30 oscillators.The value for R T is determined by the heuristic ( 25) and all parameters are given in Table 1.The results of the emulation are given in the top row and bottom left plot of Figure 6.It can be seen that after initial disarray, the oscillator voltages synchronize in two distinct groups with a phase difference of π.This is illustrated by the estimate of the spin variables s μ as well.The bottom left plot shows the cut size of W that corresponds to the spins encoded in the oscillators.It can be seen that when the oscillator voltages have converged towards synchronization, the oscillator network has found the solution to the max-cut problem.It functions thus as an Ising machine.As has already been alluded to, the comparison to Figure 5 shows that this desired behavior is impaired when deviating strongly from the found value for R T .For values that are too small, the signal shapes degenerate compared to the desired shape from Figure 3B.Values that are too large impede synchronization.
Moving on, we compare simulation times for oscillator networks of increasing size when employing our explicit model versus a model that utilizes the traditional iterative approach 8 with 5 iterations per time step.To implement the latter model, the transmission line in Figure 4 is replaced by a fractional delay.At every timestep of the simulation, a solution to the implicit relations is then approximated by iterating 5 times over the fractional delay, keeping the state of reactive elements constant.The simulated time is 100 ms for all networks and the topology of all networks is given by complete graphs.Note that we do not compare to LTspice as it is not well suited for scaling to large networks.All simulations were carried out on a lab computer (in MATLAB) with an Intel Core i7 processor running at 2.9 GHz, using 32 GB of RAM and Windows 10.The results are depicted in the bottom right plot of Figure 6.As is expected, the simulation time T sim increases with the number of oscillators n.Additionally, we can see that using the proposed WD algorithm can significantly decrease the simulation time in large networks.In fact, we see that we can reduce the simulation time by 75% in the case of large networks.We want to stress that this speedup results from the combination of both the explicit WD implementation of the nonlinearity and the introduction of the transmission lines.The former enables employing the nonlinearity without delay-free loops when coupling it to the reflection-free port of the parallel adaptor.The latter allows connecting the coupling network to a reflective port without an arising delay-free loop.We stress that while fulfilling the same function, the machine modeled with the fixed-point iteration and the one with the transmission lines are not expected to produce identical behavior.Rather the former is a more idealized version of the second.

| CONCLUSION
In this work, firstly, we have presented a general electrical model for N-shaped nonlinearities.Using the Lambert W function, we derived an explicit scattering relation in the wave digital (WD) domain with which such a nonlinearity can be emulated in real-time.To test the validity of the proposed model, we emulated a single FitzHugh-Nagumo oscillator and showed the results to be identical to LTspice.Secondly, we introduced an Ising machine with a transmission line separating the oscillators and the coupling circuit to deepen the physical modeling and prevent delay-free loops in the WD domain.We demonstrated that the resulting machine is capable of solving computationally difficult max-cut problems.We also showed that the evaluation of its WD model takes 75% less time in large oscillator networks than evaluating a traditional model using iterative techniques to break delay-free loops, while still modeling an ideal circuit exactly.These results enable efficient simulation of large oscillator networks.Future research will deal with expanding the coupling elements to be memristive, enabling self-organized topology formation and moving towards neuroinspired computing hardware.
(A) The sources e 0 and j 0 translate the curve horizontally and vertically, respectively.(B) R ν controls the slope of the diode currents that form the left and right stroke of the N-shape.(C) G 0 controls the negative differential resistance and thus the slope between the extrema.(D) e ν allows for horizontal translation of the extreme points.

F
I G U R E 3 (A) Vector-valued circuit representation of coupled FNOs.The gray dashed box marks the individual oscillators.The transmission lines model the spatial distance between the oscillators and the coupling circuit.The admittance matrix G c encodes the coupling shown in (C).(B) Output of an uncoupled FNO using WD emulation and LTspice.(C) Oscillator-based Ising machine consisting of coupled FNOs with transmission lines.The light green dashed box marks an elementary coupling circuit where the voltages u 1 and u 2 evoke the current j 12 .

F I G U R E 5
Oscillator voltages for 0:05 RT , RT and 20 RT , increasing from top to bottom.The middle plot is identical to the top left plot from Figure6.As can be seen, too large and small values for R T impede the functionality of the oscillator network as the synchronization in the middle plot is the expected behavior.

F
I G U R E 6 Top: Output voltages u μ over time (left) and the corresponding spin value s μ (right).Bottom left: Cut size converges to the max-cut over time.Bottom right: Comparison of simulation time T sim using the classical iterative and proposed method for different numbers of oscillators n.