{'fraction_non_alphanumeric': 0.06843378156138144, 'fraction_numerical': 0.0190190577537003, 'mean_word_length': 3.688072562358277, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 39, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of a dynamical system modeling lasers and applications for optical neural networks

An analytical study of dynamical properties of a semiconductor laser with optical injection of arbitrary polarization is presented. It is shown that if the injected field is sufficiently weak, then the laser has nine equilibrium points, however, only one of them is stable. Even if the injected field is linearly polarized, six of the equilibrium points have a state of polarization that is elliptical. Dependence of the equilibrium points on the injected field is described, and it is shown that as the intensity of the injected field increases, the number of equilibrium points decreases, with only a single equilibrium point remaining for strong enough injected fields. As an application, a complex-valued optical neural network with working principle based on injection locking is proposed.


Introduction
Self-sustained oscillatory systems will synchronize with an external source of periodic perturbation, given that the frequency and the strength of the injection occur within the locking range. A laser subject to external optical injection behaves the same [16]. What sets optical oscillators apart from the electronic ones is the nature of propagating electromagnetic field that has two orthogonal polarization modes which can be observed with a pair of base polarization components (meaningful reference coordinate system), be it linear, circular, or some elliptical. In following treatment, we choose to express polarization in terms of a complex amplitude E = (E − , E + ) ∈ C 2 that multiplies carrier wave of the form e −i(kx−ωt) , where k is the wave vector, x is the spatial coordinate, ω is the angular frequency, and t is the time, such that k, x, ω, t ∈ R. Coordinates E ± of E are the right (+) and left (−) circularly polarized components, they are related to the orthogonal linear components E x and E y of the electric field by Electric field emitted by a laser is where E(t) is called a slowly varying amplitude.
In absence of laser cavity anisotropies, the temporal behavior of a semiconductor laser under external optical injection can be expressed with a spin-flip rate equations [24,19] that describe the complex-valued components E ± (t) of the slowly varying amplitude E(t) as where N (t) and n(t) are real-valued functions; N is the difference between the normalized upper and lower state populations, i.e., the normalized total carrier number in excess of its value at transparency; n is the normalized imbalance between the population inversions (in reference to the populations of the magnetic sublevels), u ± are the circularly polarized components of the electric field of an external injection u = (u − , u + ) ∈ C 2 , that is, the amplitude of the external light that goes into the laser, η is the coupling efficiency factor, α is the linewidth enhancement factor that refers to saturable dispersion (Henry factor), µ is the normalized injection current, κ is the decay rate of the cavity electric field whence (2κ) −1 is the cavity photon lifetime, γ is the decay rate of the total carrier number, and γ s is the excess in the decay rate that accounts for the mixing in the carriers with opposite spins. The rate equations (1) are derived to model and explore polarization properties of Vertical-Cavity Surface-Emitting Lasers (VCSELs). The rate equations use a normalized injection current such that the unitless injection µ of 1 refers to the laser threshold operation, and µ ≈ 3 refers to the output emission of 1 mW on a typical VCSEL. In the physical world, an array of VCSELs is produced on a semiconductor wafer, where stacks of dielectric materials form high-reflectivity Bragg mirrors on the top and bottom sides of the wafer. The mirrors confine an active region in between, comprising just a few quantum wells with a thickness of some tens of nanometers. Depending on the active region diameter, the threshold current and the maximum emission power may be tailored for specific applications.
Lasers are known to exhibit a rich dynamical behavior under external optical injection [33,14,10,18,2]. Depending on laser properties and the injected optical power and its frequency, the differential equation system may converge toward an equilibrium point (a time independent solution, also called steady state, stationary point, or critical point) with locked phase synchronization. This phenomenon is called injection locking [27]. Alternatively, the system may manifest periodic oscillations, or chaos [28,7]. In this work, we explore equilibrium points of system (1) and study their stability. While in a physical system injection locking is possible only at a stable equilibrium, understanding the unstable equilibrium points provides important insight about the phase space of the system.
In our previous work [32] we concluded that in the case of linear polarization, a stably injection-locked laser approximates normalization operation that can be used for arithmetic computations. In this paper, we widen the scope and explore the equilibrium points in greater detail. Our main results regarding the dynamics of system (1) are: (i) If the injected field u ∈ C 2 is sufficiently weak (small in magnitude), then system (1) has nine equilibrium points (Theorem 2). If u is sufficiently strong (large in magnitude), then it only has a single equilibrium point (Theorem 17).
(ii) Dependence of the equilibrium points on the injected field u is described in an asymptotic sense in the limits |u| → 0 and |u| → ∞ (Theorems 2 and 17). A method for calculating the exact values of the equilibrium points is provided for weak u in terms of an ordinary differential equation (Theorem 7).
(iii) Under the assumption that α = 0 and that the injected field u is weak, it is proved that one of the nine equilibrium points is asymptotically stable, while the remaining eight are unstable (Theorem 12).
The consequence of the aforementioned results is that under weak injection of elliptically polarized light the injection-locked laser will emit linearly polarized output such that the input state of polarization is projected to a linear state of polarization (see Figure 1a). Under strong injection of elliptically polarized light, the injection-locked laser will emit light with an elliptical state of polarization, yet, the polarization is shifted toward a linear state of polarization, as shown in Figure 1b.
In the last section of this paper, we will investigate a possibility to use lasers as nodes of an optical neural network. In general, optical technologies are commonly used for linear operations, such as Fourier transformation and matrix multiplications, which come virtually free by use of lenses, mirrors, and other common light transforming elements. In this respect, optical solutions have been proposed a b Figure 1: The state of polarization is transformed by the injection-locked laser. Schematic illustration on Poincaré sphere [26]: (a) A weak injected arbitrary state of elliptical polarization (•) is projected on equator (•) by the injection-locked laser emission. (b) Under a strong elliptical state of polarization input, the state of polarization of the injection-locked output emission is shifted toward the equator, yet, will not reach it.
for matrix multiplications in optical neural networks [25,9]. However, a neural network consisting of linear transformations only is impossible, as such a network is itself linear. As recognized by the optics community, the nonlinear functions are difficult to realize in practice, as noted in recent publication Despite these positive results, the scheme faces major challenges. [...] Then there is the question of the nonlinear operation needed to link one set of [Mach-Zehnder Interferometers] with another, which [was] simply simulated using a normal computer. [6] In this respect, we propose that a laser could provide a useful nonlinearity. More specifically, a nonlinear activation function of a node is provided by injection locking; a laser nonlinearly transforms an injected field (input) into an injectionlocked emitted field (output). As the fields are complex-valued, this also leads in a natural way to a complex-valued neural network.
Complex-valued neural networks are a less studied object than their real counterpart, nevertheless, they have attracted a considerable amount of research [11,1,12]. A desired quality of any class of neural networks is the universal approximation property, namely, that any continuous function can be approximated to any degree of accuracy by a network from that class. For real-valued neural networks, necessary and sufficient conditions for an activation function to generate a class of neural networks with the universal approximation property are known [17,13], and also quantitative bounds for the approximation exist [20,34]. Besides for the theoretical expressiveness of neural networks, the choice of an activation function affects their empirical performance, as, among others, it affects the efficacy of the training algorithms [8]. In [31] we considered universality of laser based neural networks with a complex-valued activation function.
The recent universal approximation theorem for complex-valued neural networks by F. Voigtlaender [30] characterizes those activation functions for which the associated complex-valued neural networks have the universal approximation property. In this theorem, the activation function is required to be defined globally on the complex plane. As the activation function induced by injection locking is defined only locally in a neighborhood of the origin, we extend Voigtlaender's theorem by proving a local version of the universal approximation theorem (Theorem 19 stated in the Appendix). This theorem and the results about dynamics of system (1) will prove the following: The class of complex-valued optical neural networks with nodes composed of optically injected semiconductor lasers and an activation function based on injection locking has the universal approximation property, namely, it can approximate any complex-valued continuous function to any degree of accuracy (Theorem 18).
The paper is organized as follows. In Sections 2.1 and 2.2 we assume that the injected field u is weak and consider equilibrium points of system (1) and their stability, respectively. In Section 2.3 we consider the case of a strong injected field. In Section 3 we propose a design for an optical neural network with working principle based on injection locking, provide a mathematical model for such a network, and prove that these networks have the universal approximation property. In the Appendix, we prove a local version of the universal approximation theorem for complex-valued neural networks.
2 Analysis of equilibrium points and their stability

Equilibrium points with weak injected fields
In this section, we study equilibrium points of system (1) (i.e., points (E ± , N, n) at which the right-hand side of (1) vanishes) under the assumption that the injected field u is weak and constant in time. Specifically, we consider injected fields u of the form u = λ u, where u ∈ C 2 \ {0} is fixed and λ ∈ C \ {0} is a small parameter, and we are interested in the behavior of the equilibrium points as a function of the parameter λ.
(we use everywhere (a 1 , . . . , a n ) as an alternative notation for a column vector a 1 · · · a n T ). Above | · | denotes the absolute value on C and norm on C 2 , and δ := γ s /γ > 0 is a dimensionless parameter. The parameters satisfy δ, γ, κ ∈ (0, ∞), α ∈ R, and µ > 1, and throughout this paper we take them to be fixed, so that various constants explicit or implicit (as in the little o-notation) in the equations below may depend on them. Figure 2 shows an example of a solution to system (3) with an injected field u that is piecewise constant. 1 After every abrupt change of the injected field u, the solution is seen to quickly settle at a new value (an equilibrium point of the system).
Proposition 1. For every initial value (E 0 , N 0 , n 0 ) ∈ C 2 × R × R, there exists a unique maximal solution (i.e., a solution that has no proper extension that is also a solution) to system (3) satisfying the initial value at t = 0. The solution is global in forward time, that is, its domain includes [0, ∞).
Proof. A straightforward calculation shows that the right-hand side of system (3) is locally Lipschitz, which implies that for any given initial value, there exists a unique maximal solution satisfying the value at t = 0.
Consider an arbitrary maximal solution (E, N, n) : I → C 2 × R × R, where 0 ∈ I ⊂ R, and for the sake of a contradiction assume that [0, ∞) ⊂ I. If ω ∈ R denotes the right endpoint of I, then ω / ∈ I and either lim t→ω, t∈I (see [3,Theorem 7.6]). Denote ν(t) := (N (t), n(t)) ∈ R 2 . The function Y is uniformly bounded from below, in the sense that there exists c > 0 such that for every z ∈ C 2 and y ∈ R 2 it holds that Y (z)y · y ≥ c|y| 2 .
Following theorem is the main result of this section. Its essential content is that with sufficiently weak injected fields of the form u = λ u system (3) has nine distinct equilibrium points, and that the equilibrium points depend continuously on λ ∈ C \ {0} with asymptotics given by (10). In the statement of the theorem, the requirement that u − = 0 and u + = 0 means physically that the field is not circularly polarized, while | u − | = | u + | means that the field is linearly polarized. The function y : R 2 → R 2 is defined by det Y (x) = δ + (1 + δ)|x| 2 + 4x 2 1 x 2 2 > 0 (the function Y is defined in (4b)).

Theorem 2.
Consider injected external field with amplitude λ u, where λ ∈ C and u = ( u − , u + ) ∈ C 2 satisfies u − = 0 and u + = 0. There exists a constant = ( u) > 0 and a family {E with pairwise distinct values that have the following properties: (i) If in system (3) the injected field is of the form u = λ u with 0 < |λ| < , then a triple (E, N, n) ∈ C 2 × R × R is an equilibrium point (a time-independent solution) of the system, if and only if u (λ) for some j ∈ J , and (N, n) = y(|E − |, |E + |).
(ii) The functions E (j) u have following asymptotics as λ → 0: where θ := − arg(1 + iα) and (iii) Furthermore, if | u − | = | u + | and j ∈ {0, ±x}, then for every λ with 0 < |λ| < it holds that E (j) u (λ) = ρ (j) (λ) u for some ρ (j) (λ) ∈ C.    Remark 2. At the expense of a more complicated statement, the theorem can be modified to hold also in the case u − = 0 or u + = 0. The reason why this case is special is that if a point (E − , E + , N, n) is an equilibrium point of system (3) with injected field (say) u = (0, λ u + ), then for every φ ∈ R the point (e iφ E − , E + , N, n) is an equilibrium point of the system, also. Thus, instead of distinct equilibrium points, there will be disjoint sets of equilibrium points. See also Proposition 3, Remark 3, and Theorem 7 below. Figure 3 shows values of the nine equilibrium points from Theorem 2 as the magnitude λ ∈ R of an external optical injection u = λ u varies. In the dimensionless units of system (3) the intensity of the free running laser, i.e., |E| 2 at a stable equilibrium point of (3) when u = 0, is µ − 1. In the figure u has been chosen so that at |λ| = 1 the intensity |u| 2 = | u| 2 of the external injected field is also µ − 1. For the laser parameters used in the figure, the injected field is sufficiently weak in the sense of Theorem 2, namely, in the sense that the nine equilibrium points of the theorem exist, if |λ| < 0.057, i.e., if the injected field does not exceed in magnitude 5.7 % of the emitted field of the free running laser. In practice this value would depend also on experimental setup details such as the coupling efficiency.
As a real-valued amplitude E = (E − , E + ) ∈ R 2 is linearly polarized if and only if E − = ±E + , it is seen from Figure 3 that even if the injected field is linearly polarized, only three of the nine equilibrium points have a linear state of polarization, while the remaining six equilibrium points have an elliptical state of polarization.
We prove Theorem 2 at the end of this section after developing some preliminary results. We begin by transforming the problem of finding equilibrium points of system (3) from C 2 × R × R into a problem of finding solutions from R 2 to a system of two bivariate polynomials: Proposition 3. Let X and y be the functions defined in (4a) and (7).
Let φ ± ∈ R, and define a vector E ∈ C 2 and numbers N, n ∈ R by (a) Re(E + ) Figure 3: Values on the Re(E±)-plane (circularly polarized basis) of the slowly varying amplitude E of an electric field of a laser at the stable and unstable equilibrium points as the magnitude λ of an external optical injection u = λ u varies. For small |λ| the laser has nine equilibrium points (Theorem 2). Solid lines denote paths traced by real parts of the points when u = √ µ − 1(cos θ, sin θ) and θ = π/4 (linear polarization), and λ ∈ [−1/4, 1/4] varies. In this figure, u has been chosen so that the equilibrium points are real-valued for λ ∈ R and so that the intensity of the injected field u = λ u at λ = 1 is equal to the intensity of the emitted field E of the free-running laser. As λ increases, the points move in the directions indicated by the arrows. At λ = −1/4 only one of the points exists (it is located at (i)). As λ increases, eight new points appear. First at λ ≈ −0.072 two points appear at (a) and start moving in opposite directions. At λ ≈ −0.071 one of these points has moved to (b), where it splits into three. At λ ≈ −0.057 two points appear at each (c). The circled dots denote locations of the points at λ = 0. As λ grows, eight of the points disappear (at (d) (λ ≈ 0.057), (e) (λ ≈ 0.71), and (f) (λ ≈ 0.072)). The paths were calculated from the functions h (j) r (cf. Figure 4) via (12) and (13). For −1/4 ≤ λ < 0 only the equilibrium point on the path from (i) to (ii) is stable, for 0 < λ ≤ 1/4 the same is true for the equilibrium point on the path from (iii) to (iv) (cf. Figure 5). Consequently, at λ = 0, the unique stable equilibrium point of the system jumps from (ii) to (iii). The parameters κ, µ, α, γ, and δ are those of Figure 2. The dotted and dashed paths are interpreted analogously. In these paths θ ∈ {π/6, 11π/24} (elliptical polarizations).
Then the triple (E, N, n) is an an equilibrium point of system (3) with the injected electric field (ii) Suppose a triple (E, N, n) ∈ C 2 ×R×R is an equilibrium point of system (3) with some injected electric field u ∈ C 2 . Then there exists numbers φ ± ∈ R and vectors r ∈ [0, ∞) × [0, ∞) and x ∈ R 2 such that equations (11) to (13) hold.
Remark 3. An arbitrary field u = (u − , u + ) ∈ C 2 uniquely determines the numbers r j ≥ 0 in (13). If u − = 0 and u + = 0, then also the numbers e iφ ± are uniquely determined, and therefore a solution x ∈ R 2 of (11) corresponds via (12) to a unique equilibrium point of system (3). But if (say) u − = 0 and x is a solution of (11) with x 1 = 0, then there exists a continuum of equilibrium points of system (3) corresponding to x due to the arbitrary choice of φ − ∈ R in (13).

Proof of Proposition 3. For a vector
Then for every z ∈ C 2 the matrices X(z) and J φ commute, and For proving the first part of the proposition assume that x ∈ R 2 and r ∈ [0, ∞) × [0, ∞) satisfy (11), and let E = J φ x, (N, n) = y(x), and u = (1 + iα)J φ r be as in (12) and (13). Then so the point (E, N, n) is an equilibrium point of system (3) with the injected electric field u. For proving the second part of the proposition assume a point (E, N, n) ∈ C 2 × R × R is an equilibrium point of system (3) with injected electric field u ∈ C 2 , and find vectors x ∈ R 2 and φ = ( Then from above and the definition of an equilibrium point it follows that This implies (N, n) = y(x), and consequently u = (1 + iα)J φ X(y(x))x. Now define r := X(y(x))x ∈ R 2 . Then it only remains to show that r j ≥ 0, but this follows from (14), Proof. Suppose that X(y(x))x = 0, or equivalently that where Consider first the case y 1 (x)−1 = y 2 (x) = 0. Then (16) implies that x is of the form (c, 0) for some c ∈ R. To find the possible values of c, insert the candidate vector into y 1 (x) − 1 = y 2 (x) and solve for c. This shows that x ∈ {x (+l) , x (−l) }.
Consider the last case, namely y 2 (x) = 0 and y 1 (x) = 1. Then x 2 1 = x 2 2 , and inserting x = (x 1 , ±x 1 ) into y 1 (x) = 1 and solving for x 1 shows that On the other hand, a direct calculation shows that X(y(x (j) ))x (j) = 0 for every j ∈ J .
Fix nonzero r = ( r 1 , r 2 ) ∈ R 2 and define a function Our plan is to first find all zeros of F r (s, ·) for small s, and then, assuming that the injected field u in system (3) is sufficiently weak, with Proposition 3 convert these zeros to equilibrium points of the system. The Jacobian matrix of F r with respect to x will be denoted by D x F r (x) (as the Jacobian is independent of s, it is suppressed from the notation). A straightforward calculation shows that where I 2 ∈ R 2×2 is the identity matrix, (8)).
where following estimates hold for the functions a, b : Proof. Inserting the value of x (j) from (15) into the expression (18) of D x F r and inverting yields (i).
Proof. Recall that r = 0 by assumption. By (i) of Proposition 5 and the implicit function theorem there exists neighborhoods V (j) ⊂ R of 0 ∈ R and W (j) ⊂ R 2 of x (j) and smooth functions h Regarding the other direction of (23), it is enough to show that there exists > 0 such that and that F r (s, x) = 0 implies that either (s, If a pair (s, x) ∈ R × R 2 satisfies F r (s, x) = 0 and |x| > √ 2µ, then by the Pythagorean theorem (with the notation of Proposition 5) we have where the last inequality holds because a(x) > 1/2 by (19a). This implies that |s| > √ µ/( √ 2| r|), which together with the continuity of F r shows that the set is compact. By Proposition 4 the set K and the closed set {0} × R 2 are disjoint. Let d > 0 be the distance between those sets (d = +∞ if K = ∅), and consider a pair (s, x) . Consequently, if we choose > 0 small enough so that (25) and < min{d, Finally, if r 1 = r 2 and η ∈ R, then F r (s, (η, η)) = 0, if and only if The implicit function theorem shows that in some neighborhoods of (0, 0) ∈ R×R and (0, (26) implicitly defines η = η(s), and consequently, if j ∈ {0, ±x} and s is small enough, then h (j) (s) = (η(s), η(s)).
Following theorem shows that system (3) has at least nine disjoint families of equilibrium points provided that the injected field u is weak enough. These families correspond to nine distinct solutions of F r (s, ·) = 0, where s > 0 is a fixed parameter related to the strength of the field u. These solutions can be found by solving an initial value problem for an ordinary differential equation in s. As the initial value problem is easy to solve numerically, the theorem provides a computational method for obtaining numerical values for the nine families of equilibrium points.
Theorem 7. Fix u = ( u − , u + ) ∈ C 2 (with the possibility u − = 0 or u + = 0 allowed), and consider system (3) with u = λ u. Define and choose numbers φ ± ∈ R such that Let y : R 2 → R 2 and J be as defined in (7) and (9), respectively, and define θ := − arg(1 + iα). Fix j ∈ J . Suppose I ⊂ R is an interval containing the origin and is a solution to the initial value probleṁ Then for every λ ∈ C \ {0} such that |λ| ∈ I the triple (E(λ), N (λ), n(λ)) defined by is an equilibrium point of system (3) with injected field u = λ u.
Remark 4. Initial value problem (28) is straightforward to solve numerically using the explicit expressions for x (j) and D x F r given in (15) and (18), respectively. Therefore Theorem 7 provides an easy method to trace the trajectories of the equilibrium points E (j) u starting from λ = 0 for as long as |λ| is in the domain I of existence of a solution of (28). Also, the asymptotics of E (j) u as λ → 0 immediately follow from the initial value problem (28). On the other hand, if I is a finite interval, it may be possible to continue the trajectories even beyond the interval I. In that case one can use numerical continuation techniques, such as pseudo-arclength continuation, to solve the functions h (j) r (s) from (23) and use them in (29) instead (cf. Figure 4). (27), and each one of these yields an equilibrium point when plugged into (29). If u + = 0, an analogous statement holds for φ + .
We are now ready to prove Theorem 2: . The paths were solved with BifurcationKit.jl [29]. Black lines denote complement of the domain of the initial value problem (28). A solution of (28) with initial value x (j) coincides with h (j) r (s) for as long as it does not hit the boundary of the domain (at which point the right-hand side of (28a) ceases to exist). This means that the solution of (28) starting from x (−x) follows the blue path up to the point where the path first crosses the black line, and then ends there. All other paths can be solved in full from the inital value problem (28). The solutions h (j) r correspond via (29) to the equilibrium points (depicted in Figure 3) of a laser with injected external optical field. The parameters used are those of Figures 2 and 3.
Proof of Theorem 2. We will first prove that there exists a constant > 0 and nine continuous functions E (j) u , j ∈ J , that are of the form (9), for which the points (E, N, n) with are equilibrium points of system (3) with u = λ u, and that satisfy the asymptotics (10) as λ → 0. Define and let > 0 be the constant and h Note that if | u − | = | u + |, then r 1 = r 2 , and for j ∈ {0, ±x} it follows from (35) and (24) that E (j) u (λ) = ρ(λ) u for some ρ(λ) ∈ C. Fix j ∈ J . If 0 < |λ| < , then x := h (j) r (|λ|) satisfies X(y(x))x = |λ| r, and therefore from Proposition 3 it follows that a point (E, N, n) defined by (33) is an equilibrium point of system (3) with The function s → F r (s, h (j) r (s)) vanishes identically, so differentiating it and simplifying (see (30)) gives which by Proposition 5 can be solved at s = 0 to yield The matrix [D x F r (x (j) )] −1 in (37) was calculated in Proposition 5. Substituting (37) and the value of h It only remains to prove that if a triple (E, N, n) is an equilibrium point of system (3) with injected field λ u, where 0 < |λ| < , then E = E (j) (λ) for some j ∈ J , and (N, n) = y(|E − |, |E + |). To that end, consider an arbitrary equilibrium point (E, N, n) of system (3) N n = y(x), and (38c) Equalities (38b) and (38c) imply that (N, n) = y(|E − |, |E + |). Also, positivity of the components of r together with (34) and (38d) imply that r = |λ| r. Then (38a) implies that F r (|λ|, x) = 0, so x = h (j) r (|λ|) for some j ∈ J by Proposition 6. Finally, dividing the components of (38d) by their modulus shows that λ |λ| Solving for e iφ ± and inserting these values into (38b) shows that E is equal to the right-hand side of (35).

Stability of equilibrium points with weak injected fields
In this section, we consider stability properties of the nine equilibrium points from Theorem 2. We will prove in Theorem 12 below that if α = 0 and the injected field u in system (3) is sufficiently weak, then the system has exactly one asymptotically stable (in the sense of Lyapunov) equilibrium point, while the remaining equilibrium points are unstable (for the definitions of asymptotic stability and instability of an equilibrium point, we refer the reader to [3]). By splitting the complex-valued functions E ± (t) into their real and imaginary parts, i.e., writing E ± (t) = E (re) where the function f : R 6 → R 6 is determined by system (3). A calculation shows that Df , the Jacobian matrix of f , is given by the block matrix where the superscript T denotes the transpose of a matrix, and We proved in Theorem 7 a method for calculating numerical values for the nine equilibrium points of system (3) from Theorem 2. Inserting the value of an equilibrium point into the expression (40) for Df and finding the eigenvalues of the so obtained 6×6-matrix is an easy numerical method to test the stability of the equilibrium point. Recall that if all the eigenvalues of Df at an equilibrium point have strictly negative real parts, then the equilibrium point is asymptotically stable, while if at least one of the eigenvalues has a strictly positive real part, then the equilibrium point is unstable [3]. Only if none of the eigenvalues have strictly positive real parts but at least one of them has real part equal to zero, then this test for stability is inconclusive.
In Figure 5 we have used above test to determine stability of the equilibrium points on Figure 3. As illustrated in Figure 5, of the nine equilibrium points depicted in Figure 3 that correspond to an injected field u = λ u, for each λ ∈ [−1/4, 1/4] \ {0} and u exactly one of the points is asymptotically stable, while the others are unstable.
(ii) Let θ 1 and θ 2 be the two roots of the polynomial and let θ 3 and θ 4 be the two roots of the polynomial Furthermore, let E (re) ± , E (im) ± ∈ R be any numbers such that Then for the matrix the following holds: Proof. The straightforward calculation using expression (40) for Df is omitted.
Given u ∈ C 2 such that u − = 0 and u + = 0, let E (j) u , j ∈ J , be the functions from Theorem 2. By (ii) of Proposition 3, if (E, N, n) ∈ C 2 × R × R is an equilibrium point of system (3), then (N, n) = y(|E − |, |E + |). Therefore we can define functions in a punctured neighborhood of the origin of the complex plane by requiring that the point (E We call an equilibrium point (E We can now prove instability for five of the equilibrium points from Theorem 2: Lemma 9. Assume α = 0 in system (3). Fix u ∈ C 2 with u − = 0 and u + = 0, and let > 0 and E (j) u , j ∈ J , be as in Theorem 2. Then there exists 0 < 0 ≤ such that if 0 < |λ| < 0 , then the equilibrium points corresponding to E It follows that there exists 0 < 0 ≤ 1 such that if 0 < |λ| < 0 and j ∈ {0, ±l, ±r}, then at least one of the eigenvalues of (Df )   u (λ) with 0 < |λ| < 0 and j ∈ {0, ±l, ±r} has at least one strictly positive eigenvalue. Therefore the nonlinear system (3) is unstable at such a point [3,Theorem 15.6].
Let C 6 sym denote the quotient space of C 6 by the equivalence relation that identifies vectors whose coordinates are permutations of each other, and let σ : C 6×6 → C 6 sym (49) denote the map that takes a matrix to the unordered 6-tuple of its eigenvalues (repeated according to their algebraic multiplicities). Then (C 6 sym , d) is a metric space with the optimal matching distance [5] d
Our plan is to determine stability of the remaining equilibrium points corresponding to E (j) u (λ) with j ∈ {±x, ±y} by finding all eigenvalues of (Df ) (j) u (λ). In the following lemma we will first show that for small λ = 0 the eigenvalues of H (j) u (λ) approximate those of (Df ) (j) u (λ), and after that in Lemma 11 we will determine the eigenvalues of H (j) u (λ). Combining these results will then make it possible for us to conclude stability of the equilibrium points.
Lemma 10. For every j ∈ {±x, ±y}, Here d is the optimal matching distance on C 6 sym and σ is the map (49).

Proof. A calculation shows that for every
There exists numbers r > 0 and R > 0 such that if 0 < |λ| < r and j ∈ {±x, ±y}, then (Df ) is the closed ball of radius R centered at the origin. Because the continuous map σ is uniformly continuous on the compact set B R , from (52) it follows that the limit (51) holds.
We can now prove the main result of this section.
Theorem 12. Consider system (3) under the assumption that α = 0 and that the injected field u is of the form u = λ u, where λ ∈ C \ {0} and u = ( u − , u + ) ∈ C 2 satisfies u − = 0 and u + = 0. With reference to Theorem 2, let > 0 be a constant and E (j) u (λ), j ∈ J , the functions with asymptotics (10) such that for 0 < |λ| < they determine the nine equilibrium points of system (3) with injected field u = λ u.
There exists a constant 0 < 0 ≤ such that for every 0 < |λ| < 0 the equilibrium point corresponding to E  Proof. By Lemma 9 we know that the equilibrium points corresponding to E (j) u (λ) with j ∈ {0, ±l, ±r} and λ = 0 sufficiently small are unstable. By decreasing > 0 if necessary, we can assume that this is the case for all 0 < |λ| < .
To prove the theorem, we will show that for sufficiently small λ = 0 all of the eigenvalues of (Df ) (+x) u (λ) have strictly negative real parts, and that at least one of the eigenvalues of each of (Df ) (j) u (λ) with j ∈ {−x, ±y} has a strictly positive real part. By [3,Theorem 15.6] this will imply the result.
Let θ i , i ∈ {1, 2, 3, 4}, be the roots of the polynomials (43) and (44) in Lemma 8. Because all of the coefficients in the polynomials are strictly positive, Re θ i < 0 for every i. Therefore it is possible to find a radius r > 0 such that ∪ 4 i=1 B r (θ i ) ⊂ C − := {z ∈ C : Re z < 0}, and such that this union is disjoint from B r (0). Here B r (z) ⊂ C denotes the open disk of radius r centered at z ∈ C.
Fix j ∈ {±x, ±y}. By Lemmas 10 and 11 and the definition of the optimal matching distance d, we can find 0 < 1 ≤ such that if 0 < |λ| < 1 , then (Df ) (j) u (λ) has two eigenvalues in B r (0) and four eigenvalues in u (λ)) = 0, so by (i) of Lemma 8 the two eigenvalues of (Df ) only the two eigenvalues (53) are relevant for determining the stability for small λ = 0. Consider Theorem 7 and let h be a solution to the initial value problem (28). For sufficiently small λ = 0 let E(λ), N (λ) and n(λ) be defined in terms of h by (29). Then by Theorem 2 the vector E(λ) is equal to E (k) u (λ) for some k ∈ J , and an inspection shows that k = j is the only possibility. If y 1 and y 2 are the component functions of the function y from (7), i.e., y(x) = (y 1 (x), y 2 (x)), above implies that The functions s → y k • h(s) are defined and differentiable in a neighborhood of the origin, and where r = (| u − |, | u + |) ∈ (0, ∞) × (0, ∞). Calculating the gradient and applying the value of [D x F r (x (j) )] −1 obtained in (i) of Proposition 5 to (55), we can calculate that The numbers in the parenthesis on the right-hand sides of (56) and (57) are nonzero and positive. If j = +x, then x Therefore for these λ these two eigenvalues of (Df )

Equilibrium points with strong injected fields
In this section, we consider equilibrium points of system (3) under the assumption that the injected electric field u is strong (large in magnitude). We assume that the injected field is of the form where λ ∈ C is a large parameter and u = ( u − , u + ) ∈ C 2 satisfies u − = 0 and u + = 0, and we are interested in the behavior of the equilibrium points as a function of the parameter λ. For a number 0 < η < 1 and a vector r ∈ R 2 such that r 1 > 0, r 2 > 0, and | r| = 1, let us define the compact set K(η, r) := w ∈ R 2 : w · r ≥ η|w| and 1 2 ≤ |w| ≤ 3 2 .
We will prove that given the vector r ∈ R 2 , we can choose a number η = η( r) ∈ (0, 1) and a constant L = L( r) > 0 so that for the function F r defined in (17) the following holds: If s ≥ L, then Recall that by Proposition 3 the zeros of F r (s, ·) and the equilibrium points of system (3) are in one-to-one correspondence. Once (i) and (ii) are proved, we can conclude from (i) that for s ≥ L every zero of F r (s, ·) is contained in sK(η, r), and from (ii) and the Banach fixed-point theorem that there exists exactly one such zero in sK(η, r). From this it follows that if the injected field u is strong enough, then there exists a unique equilibrium point of system (3).
Given s > 0 and r ∈ R 2 , define a mapping G s r : Obviously, for every s > 0 the set of zeros of F r (s, ·) and the set of fixed points of G s r coincide.
Below D x G s r denotes the Jacobian matrix of the map G s r defined in (62).
Proposition 16. Let r ∈ R 2 satisfy (58). There exists a constant L = L( r) > 0 such that following hold: (i) For every s ≥ L the function G s r has a unique fixed point in R 2 .
(ii) If h : [L, ∞) → R 2 denotes the function that maps s to the unique fixed point of G s r , then h is differentiable on (L, ∞).
(iii) There exists a constant C > 0 (independent of r) such that the function h from (ii) satisfies Proof. Let η = η( r) ∈ (0, 1) and L = L( r) > 0 be such that for s ≥ L inequality (64) holds for every x, x ∈ sK(η, r) and 0 ≤ ν ≤ 1. If necessary, increase L so that in addition for s ≥ L inclusion (63) holds and equality F r (s, x) = 0 implies that x ∈ sK(η, r) (cf. Lemma 13). Let s ≥ L. Then G s r maps sK(η, r) into itself, and if x, x ∈ sK(η, r), applying the fundamental theorem of calculus and estimating with (64) shows that Thus, the restriction of G s r to sK(η, r) is a contraction. By the Banach fixed-point theorem the function G s r has a unique fixed point in sK(η, r). Because G s r (x) = x if and only if F r (s, x) = 0, this fixed point is unique in R 2 , also. Part (i) is now proved.
Let s 0 > L and h be as in (ii). Consider the function (L, ∞) × R 2 (s, x) → F r (s, x) ∈ R 2 at a neighborhood of its zero (s 0 , h(s 0 )). Since it follows from inequality (64) that at the point (s, x) = (s 0 , h(s 0 )) the derivative D x F r (x) is invertible. Then by the implicit function theorem in some neighborhood (s 0 − , s 0 + ) the zero of F r (s, ·), i.e., h(s), depends differentiably on s. Because s 0 > L was arbitrary, the function s → h(s) is differentiable, and (ii) is proved.
If we write h(s) as in (66) and denote x := h(s), then in the notation of Proposition 5 we have Because 1/2 ≤ |x|/s ≤ 3/2 since x ∈ sK(η, r), we obtain from (19a) and (19b) that for some constant C > 0 depending only on µ it holds that |e(s)| ≤ C/s 2/3 , for every s ≥ L. This proves (iii).
With the previous proposition in hand, we can now prove the main theorem of this section. Note that, among others, the theorem states that unlike in the case of weak injected fields, in which case system (3) has nine equilibrium points (Theorem 2), in the case of strong injected fields, the system has a single equilibrium point.
Theorem 17. Consider u = ( u − , u + ) ∈ C 2 with u − = 0 and u + = 0. There exists a constant L = L( u) > 0 and a continuous function with the following property: If in system (3) the injected field u is of the form u = λ u with |λ| ≥ L, then a triple (E, N, n) ∈ C 2 × R × R is an equilibrium point of the system, if and only if E = E u (λ) and (N, n) = y(|E − |, |E + |) (the function y is defined in (7)). Furthermore, there exists a constant C = C( u) > 0 such that the function E u satisfies where |e(λ)| ≤ C |λ| 2/3 and θ := − arg(1 + iα). (67) Remark 6. It follows from (67) that the magnitudes of the emitted field E u (λ) and the injected field u = λ u are asymptotically related by and that as λ grows, the polarization of the emitted field E u (λ) approaches on the normalized Poincaré sphere that of u.

Optical neural networks based on injection locking
We now describe a design of an optical neural network that can be implemented with a network of lasers, and whose working principle is based on injection locking (see Figure 6a). The network consists of an input layer (Layer I), an output layer (Layer K), and one hidden layer (Layer J) in between (the working principle naturally generalizes to a network with several hidden layers): (i) In the input layer, each node (artificial neuron) is a laser. The nodes in this layer are not connected to each other, and the output of a node is the electric field emitted by the corresponding laser.
(ii) In the hidden layer, the nodes are lasers that are coupled to injected electric fields. The injected fields are composed of fixed external electric fields together with outputs of the input layer modified by some passive optical elements, e.g., polarizers or mirrors, optical isolators, and absorbing components. Due to injection locking, each laser in the hidden layer stabilizes to some equilibrium point determined by the injected field, and the output of a node is the emitted electric field.
The coupling between layers I and J is unidirectional, we note that one can use lasers of varying powers to replace the use of optical isolators.
(iii) Between the hidden layer and the output layer, the electric fields from the hidden layer are first modified by passive optical elements, and then joined to form the output of the network. The nodes in the output layer correspond to exits of optical cables or waveguides in integrated optics.
. The relation between inputs and outputs of the network is set by choosing the external electric fields that are part of the injected fields in the hidden layer, and the passive optical elements on both sides of the hidden layer. We will show that  an arbitrary continuous function can be approximated within any given accuracy by networks of this form.
The optical neural network is modeled mathematically as follows. Indexes of lasers in the input layer are denoted by I = {1, 2, . . . , i 0 }. The output of i:th laser is a linearly polarized electric field E (I) i ∈ C 2 , and all electric fields in this layer are assumed to share the same linear polarization, i.e., for all i = 1, 2, . . . , i 0 , where λ (I) i ∈ C, and u = ( u − , u + ) ∈ C 2 \{0} is fixed and satisfies | u − | = | u + |. It is also assumed that the set of all possible inputs is bounded, i.e., there exists R > 0 such that whenever (λ is an input to the network, then |(λ In the hidden layer indexes of lasers are denoted by J = {1, 2, . . . , j 0 }. The passive optical elements between the input layer and the hidden layer may induce scaling and phase shift to the electric fields, i.e., field E (I) i from the i:th laser of the input layer to the j:th laser of the hidden layer transforms to a ji E (I) i , where a ji ∈ C. The total injected field u j ∈ C 2 to the j:th laser in the hidden layer is then the sum of the modified fields and an external electric field E (ext) j , which is assumed to share the same polarization with the lasers in the input layer: By Theorems 2 and 12, if the linewidth enhancement factor α of the laser is zero (i.e., α = 0 in system (3)) and the injected field u j to the j:th laser is written as u j = λ (J) j u, then for some constant > 0 it holds that as long as 0 < |λ (J) j | < , then the j:th laser has a unique stable equilibrium point (denoted by E where the function is defined in Theorem 2. Figure 6b illustrates the function ρ corresponding to the system in Figure 2. In the output layer nodes are indexed by K = {1, 2, . . . , k 0 }, and the k:th output E (K) k ∈ C 2 of the network is a superposition of the emitted fields E (J) j of lasers in the hidden layer modified by passive optical elements represented by complex numbers c kj : whenever 0 < |λ (J) j | < for all j = 1, 2, . . . , j 0 . As the input to the network is of the form (λ i ∈ C, and the output is by (73) of the form (λ ∈ C, the network essentially computes the map (λ It follows from equations (70)-(73) that the k: where it is assumed that In (74) and (75) parameters a ji , c kj ∈ C correspond to the passive optical elements between the layers, and parameters b j ∈ C correspond to the fixed external electric fields.
Remark 7. The lasers in the input layer are not connected with each other, yet, the formulation assumes that the phase differences remain constant at the equilibrium point. As known, all oscillatory signal sources, lasers included, fluctuate in phase. This drift will inevitably invalidate the assumption of the constant phase difference between two lasers unless they share a common reference (seed) signal. Therefore, a practical implementation of a laser-based optical neural network will require a common narrow-linewidth reference signal that is used to lock enough lasers in the network. At the bare minimum, all lasers of the first layer must be injected from the same source. The phase of the injected reference light may be controlled individually for each network node, but the natural fluctuations of the reference must be experienced equally among the injected lasers. This arrangement is not unlike the clock signal of a digital computer that is used to synchronize operations between individual circuits.
Below B R ⊂ C i 0 is the closed ball of radius R centered at the origin.
Theorem 18. Fix integers i 0 > 0 and k 0 > 0 and a number R > 0, let ρ be as in (72), and consider an arbitrary continuous function f : B R → C k 0 . Let > 0. There exists an integer j 0 > 0 and numbers a ji , b j , c kj ∈ C, j = 1, 2, . . . , j 0 , i = 1, 2, . . . , i 0 , k = 1, 2, . . . , k 0 , such that following holds: (i) The inequalities (75) hold for a.e. (λ In particular ρ is not a.e. equal to a continuous function, and consequently it satisfies both (i) and (ii) of Theorem 19 stated in the Appendix (note that if ∆ m ρ ≡ 0 for some m ∈ N in the sense of distributions, then ρ is a.e. equal to a smooth function by elliptic regularity [23]). Let f : B R → C k 0 be a continuous function and fix > 0. By Theorem 19 there exists an integer j 0 > 0 and parameters a ji , b j , c kj ∈ C such that i 0 i=1 a ji λ i + b j ∈ U for every j = 1, 2, . . . , j 0 and (λ i ) i 0 i=1 ∈ B R , and such that the network N : B R → C k 0 defined componentwise by (77) satisfies Furthermore, it may be assumed that for every j either (a j1 , a j2 , . . . , a ji 0 ) = 0 or b j = 0, since otherwise the corresponding term does not affect the value of N . Observe that N is measurable, because the set

Appendix: Approximation theorem for complex-valued neural networks
In this appendix, we generalize the recent universal approximation theorem for complex-valued neural networks by F. Voigtlaender [30] to the case of activation functions defined locally in an open subset U ⊂ C, instead of globally on the whole complex plane. The gist of the proof, namely the use of Wirtinger calculus [15] to show that the functions z α z β (z is the complex conjugate of z) can be approximated by neural networks, is the same as in the proof of Voigtlaender's theorem. However, the proof is complicated by the fact that parameters for the network need to be chosen so that all inputs to the activation function stay within U .
Let B R := {z ∈ C i 0 : |z| C i 0 ≤ R}. We consider (shallow) complex-valued neural networks N : B R → C k 0 , whose k:th component function is of the form where a j · z := i a ji z i . Here the integers i 0 > 0, j 0 > 0, and k 0 > 0 are the number of inputs of the network, the width of the network, and the number of outputs of the network, respectively, and ρ : U → C, where U ⊂ C is an open set, is the activation function. The parameters a j = (a j1 , a j2 , . . . , a ji 0 ) ∈ C i 0 , j = 1, 2, . . . , j 0 , b ∈ C j 0 , and (c kj ) ∈ C k 0 ×j 0 are required to satisfy a j · z + b j ∈ U for every z ∈ B R and j = 1, 2, . . . , j 0 .
Following theorem is a local version of Voigtlaender's universal approximation theorem for complex-valued neural networks [30, Theorem 1.3]: Theorem 19. Let i 0 , k 0 , R, and ρ be as above, and suppose that (i) ρ is locally bounded and continuous almost everywhere in the nonempty open set U ⊂ C = R 2 (the measure is the two-dimensional Lebesgue measure), and (ii) ∆ m ρ does not vanish identically in U for any m = 0, 1, 2, . . . (here ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 , z = x + iy, is the Laplace operator defined in the sense of distributions).
If f : B R → C k 0 is continuous and > 0, then there exists an integer j 0 > 0 and parameters a j ∈ C i 0 , j = 1, 2, . . . , j 0 , b ∈ C j 0 , and (c kj ) ∈ C k 0 ×j 0 such that (78) holds, and that the complex-valued neural network N defined componentwise by (77) satisfies There is a slight difference in the continuity assumption for the activation function ρ between Theorem 19 and [30, Theorem 1.3]. Here we require that ρ is continuous almost everywhere, i.e., that the set D ⊂ C of its discontinuities is a null set. In [30] it is required that also the closure of D is a null set. The difference is due to how the (potentially nonsmooth) activation function is smoothly approximated; our approximation method is contained in the following two lemmas. Our approach is similar to [13,Lemma 4], in which real-valued activation functions are considered. Theorem 19 will be proved after the lemmas.
Lemma 20. For η > 0, let P(η) denote the set of countable partitions of R 2 into measurable subsets with diameter at most η, and let ψ : R 2 → C be a bounded and almost everywhere continuous function with compact support. Then where λ 2 denotes the Lebesgue measure on R 2 .
Proof. Choose a sequence of partitions ((C j (k)) ∞ j=1 ) ∞ k=1 ∈ P(1/k), and define where 1 C j (k) is the characteristic function of the set C j (k).
The functions d k are measurable, uniformly bounded by 2 ψ ∞ , and they are all supported in a fixed compact set. If x ∈ R 2 is a point of continuity of ψ, then d k (x) → 0. As a consequence, d k → 0 as k → ∞ almost everywhere in R 2 , and by the Lebesgue's dominated convergence theorem This proves the lemma as the sequence ((C j (k)) ∞ j=1 ) ∞ k=1 ∈ P(1/k) was arbitrary. Namely, if (80) did not hold, it would be possible to construct a sequence ((C j (k)) ∞ j=1 ) ∞ k=1 ∈ P(1/k) for which (81) fails.

Lemma 21.
Consider ϕ ∈ C c (R 2 ) and let ψ be as in Lemma 20. Then Proof. We can estimate The sum in A can be bounded from the above by By Lemma 20 this tends to zero as h → ∞.
The number of nonzero terms in B is bounded from the above by C/h 2 , where C > 0 is a constant independent of h. Consequently, B can be estimated from the above by C sup{|ϕ(z) − ϕ(z )| : |z − z | 2 ≤ 2h 2 }, which tends to zero as h → 0 by the uniform continuity of ϕ.
Proof of Theorem 19. It is enough to consider the case with a single output (k 0 = 1), for the general case follows from a componentwise construction of N .
For any parameters (a, b) ∈ C i 0 × U such that a · z + b ∈ U for every z ∈ B R , define a bounded function f a,b : B R → C by setting f a,b (z) := ρ(a · z + b). Then define Σ(ρ) := span{f a,b : (a, b) ∈ C i 0 × U satisfies (82)} ⊂ B(B R ).
Here B(B R ) is the complex algebra of bounded functions on B R equipped with the supremum norm, and the closure of the span is with respect to that norm. The theorem will be proved by showing that Σ(ρ) includes the subset of continuous functions of B(B R ). Let ϕ be a mollifier on R 2 and define ϕ p (s) := p 2 ϕ(ps) for p = 1, 2, . . . Fix an integer m ≥ 0 and find open sets V and W such that ∅ = V ⊂⊂ W ⊂⊂ U and that ∆ m ρ does not vanish identically in V . Let χ ∈ C c (U ) be such that χ ≡ 1 on W . The convolution (χρ) * ϕ p (s) := R 2 (χρ)(s − y)ϕ p (y) dy is then defined everywhere, and (χρ) * ϕ p | V → ρ| V as p → ∞ in the sense of distributions in V . Consequently, there exists an index p 0 such that V − supp ϕ p 0 ⊂ W and ∆ m (χρ) * ϕ p 0 does not vanish identically in V . Define ρ : C → C by ρ(s) := (χρ) * ϕ p 0 (s). Then ρ is smooth everywhere (in the sense of real differentiability), and ∆ m ρ does not vanish identically in V .
Consequently, (85) is true for every α and β. Because ∆ m ρ = (4∂∂) m ρ does not vanish identically in V , for every α and β such that |α| ≤ m and |β| ≤ m there exists b α,β ∈ V such that ∂ |α| ∂ |β| ρ(b α,β ) = 0. Then (84) and (85) with a = 0 and b = b α,β imply that z α z β ∈ Σ(ρ). Consequently, Σ(ρ) contains all functions of the form where z ∈ B R , m ∈ N and c αβ ∈ C are arbitrary. Functions of the form (87) form a self-adjoint algebra of continuous complex functions on the compact set B R , and that algebra separates points on B R and vanishes at no point of B R . By the Stone-Weierstrass theorem [22] such an algebra contains all continuous complex functions in its uniform closure, and therefore so does Σ(ρ).