*Article* **Orbital Instability of Chaotic Laser Diode with Optical Injection and Electronically Applied Chaotic Signal**

#### **Satoshi Ebisawa 1,2,\* and Shinichi Komatsu <sup>2</sup>**


Received: 21 February 2020; Accepted: 25 March 2020; Published: 30 March 2020

**Abstract:** We numerically studied the chaotic dynamics of a laser diode (LD) system with optical injection, where a chaotic signal, which is generated by an LD with optical feedback, is applied to the drive current of the master LD. To quantify the orbital instability of the slave LD, the Lyapunov exponent was calculated as a function of the optical injection ratio between the master and slave LDs and the optical feedback ratio of the applied signal. We found that the Lyapunov exponent was increased and the orbital instability was enhanced by applying a chaotic signal when the inherent system without the applied signal was in a "window". Next, we investigated the orbital instability of the slave LD in terms of statistical and dynamical quantities of the applied chaotic signal. The maximal value of the Lyapunov exponent for a certain range of the injection ratio was calculated and we showed that a chaotic pulsation is suitable for enhancing the orbital instability of the LD system. We then investigated chaos synchronization between the LDs. It is concluded that the orbital instability of an LD with optical injection can be enhanced by applying chaotic pulsation without chaos synchronization.

**Keywords:** laser chaos; semiconductor laser; chaotic laser diode; optical injection

#### **1. Introduction**

Since the chaotic oscillation of a laser diode (LD) [1–7] has a high frequency and broad bandwidth, its potential applications, such as physical random bit generation [8–11], reservoir computing [12–14], decision-making [15], and chaotic communication [16–19], have been widely studied. In these applications, more chaotic oscillation can contribute to increasing the performance, for example, randomness in random bit generation, a high bit rate in security in chaotic communication. Various methods of generating chaotic oscillation with a broad bandwidth and large chaotic property have been studied, for example, using a master–slave LD system with frequency detuning [20,21], an external feedback system with dual feedback [22], or an external feedback system with random feedback [23].

We previously proposed a method using a master–slave LD system with a random signal applied to the drive current of the LDs [24]. In the optical injection system, which consists of master and slave LDs, various dynamics of the slave LD appear. For a small optical injection ratio, the slave laser oscillates stably, periodically or quasi-periodically. Then, the dynamics develop into a chaotic state with increasing optical injection ratio, and periodic oscillation is observed between chaotic states, which is called a "window". In a window, the chaotic dynamics are concealed. We have shown numerically that the chaotic dynamics are revealed by applying a pseudorandom signal to the drive current of the master LD, and the chaotic property, that is, the orbital instability of the slave LD, is enhanced by increasing the standard deviation of the applied random signal. The orbital instability of a chaotic system can be controlled by applying a statistical random signal to a deterministic chaotic system.

In this work, a deterministic chaotic signal is adopted as the applied signal. We numerically investigate an optical injection system with unidirectional coupling from a master LD to a slave LD by applying a signal, which is generated by the chaos source LD with external optical feedback, to the drive current of the master LD. First, we compare the system with an applied chaotic signal and the system with an applied random signal having the same mean and standard deviation as the chaotic signal. It is shown that, in the window, the chaotic dynamics of the slave LD are revealed by the applied chaotic signal as well as the applied pseudorandom signal. Moreover, the applied chaotic signal more greatly enhances the orbital instability of the slave LD than the applied pseudorandom signal. Next, to explore the factor causing the enhanced orbital instability of the slave LD, we estimate the orbital instability of the slave LD in terms of statistical and dynamical quantities of the applied chaotic signal. Then, we discuss the suitable conditions of the applied chaotic signal for enhancing the orbital instability of the LD system, and the chaos synchronization between the applied signal and LD system.

#### **2. Chaotic Laser System and Lyapunov Exponent**

We consider the optical injection system consisting of two laser diodes (LDs), that is, a master LD (LD1) and a slave LD (LD2) in Figure 1a, which are driven by a DC source. The optical coupling from LD2 to LD1 is restricted by an optical isolator (ISO) and the coupling ratio is controlled by a variable attenuator (VA). An external signal is electronically applied to the drive current of LD1, which is generated by an external applied signal source, and the amplification of the applied signal is controlled by a variable electric attenuator and an amplifier. In the following sections, we consider three kinds of applied signals, that is, a chaotic signal, a pseudorandom signal and a DC. The applied pseudorandom signal, and applied DC are generated by a signal generator, and the chaotic signal is generated by an external chaos source LD (LD0) with optical feedback whose ratio is controlled by VA (Figure 1b). The chaotic signal is detected and converted into electric signal by a photo detector (PD). Since actual electric circuits have a frequency response and a cutoff frequency, impacts of the frequency band of the applied signal are needed to consider like Refs. [24,25]. In this work, we ignore the frequency response of the electric circuit to focus on the impacts of chaotic signal. The dynamics of LD0, LD1 and LD2 are described by the following rate equations [26,27]:

$$\frac{\mathrm{d}A\_1(t)}{\mathrm{d}t} = \frac{1}{2} G\_N u\_1(t) A\_1(t),\tag{1}$$

$$\frac{d\Phi\_1(t)}{dt} = \frac{1}{2} aG\_N n\_1(t),\tag{2}$$

$$\frac{\mathrm{d}n\_1(t)}{\mathrm{d}t} = [1 + \mathrm{g} \cdot \mathrm{C}(t)](p - 1)I\_{\mathrm{th}} - \gamma n\_1(t) - [\Gamma + \mathrm{G}\_N n\_1(t)]A\_1^2(t),\tag{3}$$

$$\frac{\mathrm{d}A\_{2}(t)}{\mathrm{d}t} = \frac{1}{2} G\_{N} n\_{2}(t) A\_{2}(t) + \kappa\_{\mathrm{in}\dot{\eta}} A\_{1}(t - \tau\_{\mathrm{in}\dot{\eta}}) \cos[\omega \tau\_{\mathrm{in}\dot{\eta}} + \phi\_{2}(t) - \phi\_{1}(t - \tau\_{\mathrm{in}\dot{\eta}})],\tag{4}$$

$$\frac{d\phi\_2(t)}{dt} = \frac{1}{2}aG\_N n\_2(t) - \kappa\_{\text{in}\mathfrak{j}}\frac{A\_1(t-\tau\_{\text{in}\mathfrak{j}})}{A\_2(t)}\sin[\omega\tau\_{\text{in}\mathfrak{j}} + \phi\_2(t) - \phi\_1(t-\tau\_{\text{in}\mathfrak{j}})],\tag{5}$$

$$\frac{d n\_2(t)}{dt} = (p - 1)l\_{\text{th2}} - \gamma n\_2(t) - \left[\Gamma + G\_N n\_2(t)\right] A\_2^2(t),\tag{6}$$

*Photonics* **2020**, *7*, 25

$$\frac{\mathrm{d}A\_{0}(t)}{\mathrm{d}t} = \frac{1}{2} G\_{\mathrm{N}} n\_{0}(t) A\_{0}(t) + \kappa\_{\mathrm{fb}} A\_{0}(t - \tau\_{\mathrm{fb}}) \cos[\omega \tau\_{\mathrm{fb}} + \phi\_{0}(t) - \phi\_{0}(t - \tau\_{\mathrm{fb}})],\tag{7}$$

$$\frac{d\Phi\_0(t)}{dt} = \frac{1}{2}aG\_N n\_0(t) - \kappa\_{\rm fb} \frac{A\_0(t-\tau\_{\rm fb})}{A\_0(t)} \sin[\omega \tau\_{\rm fb} + \phi \imath(t) - \phi \imath(t-\tau\_{\rm fb})],\tag{8}$$

$$\frac{d n\_0(t)}{dt} = (p-1)l\_{\rm th} - \gamma n\_0(t) - \left[\Gamma + G\_N n\_0(t)\right] A\_0^2(t),\tag{9}$$

where *A*(*t*), *φ*(*t*), and *n*(*t*) are the amplitude, the phase of the laser field, and the carrier number above the value for the solitary LD, respectively. The subscripts 1, 2, and 0 denote LD1, LD2, and LD0, respectively. *GN* is the differential optical gain, *α* is the linewidth enhancement factor, *γ* is the carrier decay rate, and Γ is the cavity decay rate. The angular frequency of the solitary LD is described as *ω* = 2*πc*/*λ*, where *c* is the velocity of light and *λ* is the wavelength. The drive current of the system without an applied signal is expressed as *p J*th.

Equations (1)–(3) describe the dynamics of LD1. The external signal is applied to ensure that the drive current of LD1 is above the threshold [1 + *g* · *C*(*t*)](*p* − 1)*J*th, where *A*<sup>0</sup> is the amplitude of LD0, *<sup>g</sup>* is the amplification coefficient, and *<sup>C</sup>*(*t*) = *<sup>a</sup>* · *<sup>A</sup>*0(*t*)<sup>2</sup> represents the applied signal for LD1, which is normalized by the parameter *a*. Equations (4)–(6) describe the dynamics of LD2, which has the optical injection from LD1. The second terms on the right side of Equations (4) and (5) describe the optical injection from LD1 to LD2. *A*1(*t* − *τ*inj) and *φ*1(*t* − *τ*inj) are the amplitude and phase of the laser field injected into LD2 from LD1, respectively. Equations (7)–(9) describe the dynamics of LD0, which has the optical feedback. The second terms on the right side of Equations (7) and (9) describe the optical feedback for LD0. *A*0(*t* − *τ*fb) and *φ*0(*t* − *τ*fb) are the amplitude and phase of the laser field fed back from the external cavity to LD0, respectively. *τ*inj is the injection time from LD1 to LD2, and *τ*fb is the round-trip time of the external cavity for LD0. The injection and feedback coefficients are expressed as *<sup>κ</sup>*inj = (<sup>1</sup> − *<sup>r</sup>*<sup>2</sup> <sup>0</sup>)*r*inj/*r*0*τ*in and *<sup>κ</sup>*fb = (<sup>1</sup> − *<sup>r</sup>*<sup>2</sup> <sup>0</sup>)*r*fb/*r*0*τ*in, respectively, where *r*inj is the injection ratio of the output injected into LD2 to the output of LD1, *r*fb is the feedback ratio of the output fed back from the external cavity to LD0, and *τ*in is the round-trip time in the inner cavity. In our simulation using the Runge–Kutta method, where the step size is 1ps, the following values are assigned to the parameters, which are taken from Ref. [26]: *GN* = 2.142 × <sup>10</sup>4[s−1], *<sup>α</sup>* = 5.0, *<sup>λ</sup>* = <sup>635</sup>[nm], *<sup>c</sup>* = 3.0 × 108[m/s], *<sup>γ</sup>* = 0.909[ns−1], <sup>Γ</sup> = 0.357[ps−1], *<sup>r</sup>*<sup>0</sup> = 0.556, *<sup>τ</sup>*in = 8.0[ps−1], *<sup>N</sup>*sol = 1.708 × 108, *<sup>τ</sup>*fb = 5.0[ns] and *<sup>τ</sup>*inj = 5.0[ns]. The initial values *<sup>A</sup>*(0) and *<sup>n</sup>*(0) are the convergent values of the solitary LD, and *φ*(0) = 0 is utilized. Then, the pseudorandom signal is generated by the Mersenne Twister random number generator [28] and Box–Muller transform [29].

**Figure 1.** Schematic diagram of the optical injection LD system with an applied chaotic signal which consists of (**a**) master-slave LD system and (**b**) applied signal source.

*Photonics* **2020**, *7*, 25

In this study, the maximal Lyapunov exponent is estimated to quantify the orbital instability of the chaotic LD. We describe how to estimate the maximal Lyapunov exponent by linear stability analysis [30–32]. When we estimate the Lyapunov exponent of LD2, the small variations *δA*2(*t*), *δφ*2(*t*) and *δn*2(*t*) of the dynamic variables of Equations (4)–(6) from the reference orbit, respectively written as *As*2(*t*), (*ωs*2(*t*) − *ω*)*t* and *ns*2(*t*), are considered. Since LD2 is the optical injection system, *δA*2(*t*), *δφ*2(*t*), and *δn*2(*t*) for LD2 satisfy

$$
\begin{pmatrix}
\frac{\mathrm{d}\delta\_{A2}(t)}{\mathrm{d}t} \\
\frac{\mathrm{d}\delta\_{\Phi2}(t)}{\mathrm{d}t} \\
\frac{\mathrm{d}\delta\_{n2}(t)}{\mathrm{d}t}
\end{pmatrix} = J\_{\mathrm{inj}} \begin{pmatrix}
\delta\_{A2}(t) \\
\delta\_{\Phi2}(t) \\
\delta\_{n2}(t)
\end{pmatrix}.
\tag{10}
$$

Here, *J*inj is the Jacobian matrix of order 3 × 3, and is given in the Appendix A. In this work, the time delay terms, *A*1(*t* − *τ*inj) and *φ*1(*t* − *τ*inj), are dealt with as external parameters since these parameters are not the dynamic variables of LD2 but those of LD1; in other words, the dynamics of LD2 is approximated using only three variables of LD2. On the other hand, when we estimate the Lyapunov exponent of LD0, since LD0 is the optical feedback system, the small variations *δA*0(*t*), *δφ*0(*t*), *δn*0(*t*), *A*0(*t* − *τ*fb) and *φ*0(*t* − *τ*fb) of Equations (7)–(9) from the reference orbit, respectively written as *As*0(*t*), (*ωs*0(*t*) − *ω*)*t*, *ns*0(*t*), *As*0(*t* − *τ*fb) and (*ωs*0(*t* − *τ*fb) − *ω*)(*t* − *τ*fb), are considered. Then, *δA*0(*t*), *δφ*0(*t*), *δn*0(*t*), *δA*0(*t* − *τ*fb) and *δφ*0(*t* − *τ*fb) satisfy

$$\begin{pmatrix} \frac{\mathbf{d}\delta\_{A0}(t)}{\mathbf{d}t} \\\\ \frac{\mathbf{d}\delta\_{\phi0}(t)}{\mathbf{d}t} \\\\ \frac{\mathbf{d}\delta\_{n0}(t)}{\mathbf{d}t} \end{pmatrix} = J\_{\mathbf{fb}} \begin{pmatrix} \delta\_{A0}(t) \\\\ \delta\_{\phi0}(t) \\\\ \delta\_{n0}(t) \\\\ \delta\_{A0}(t - \tau\_{\mathbf{fb}}) \\\\ \delta\_{\phi0}(t - \tau\_{\mathbf{fb}}) \end{pmatrix},\tag{11}$$

where *J*fb is the Jacobian matrix of order 3 × 5, and is given in the Appendix A. These equations are calculated numerically using the Runge–Kutta method, where the step size is 1 ps, and the norm *Dj* = \* ∑*t* - *δ*2 *Ai*(*t*) + *<sup>δ</sup>*<sup>2</sup> *φi* (*t*) + *δ*<sup>2</sup> *ni*(*t*) (*i* = 0, 2 *j* = 1, 2, 3, ···) is calculated by the method in Refs. [33,34]. The subscript *j* indicates the time section [(*j* − 1)*τ*, *jτ*) and the term in the square root is the summation in the range of [(*j* − 1)*τ*, *jτ*), where *τ* indicates the injection time *τ*inj for LD2 or the feedback time *τ*fb for LD0. Since the norm between the chaotic orbit and the reference orbit is gradually large and the local approximation can not be used, we initialize and replace the small variation *δAi*(*jτ*), *δφi*(*jτ*) and *δni*(*jτ*) with *δAi*(*jτ*)/*Dj*, *δφi*(*jτ*)/*Dj* and *δni*(*jτ*)/*Dj*, respectively, at intervals of *τ*. The rate of increase in the norm is considered and the Lyapunov exponent is represented by

$$
\lambda\_{\rm LSA} = \frac{1}{N\pi} \sum\_{j=1}^{N} \ln \frac{D\_{j+1}}{D\_j}. \tag{12}
$$

We use the discrete optical outputs *Ai*(*t*), *φi*(*t*) and *ni*(*t*) (*i* = 0, 1, 2), which are sampled at intervals of 10 ps over 5 μs. Here, to show the robustness of *λ*LSA against initial conditions, we investigate *λ*LSA plotted against length of time series for the calculation of *λ*LSA . We consider the LD used in the figure of Section 3.2, which has the parameters *r*inj = 0.06 and *r*fb = 0.10 as a typical example, and show the mean of *λ*LSA in Figure 2. The error bars represent the standard deviations. The number of this population is a thousand and the initial value of *A*2(*t*) are given randomly in the range of [0.9 × *A*2(0), 1.1 × *A*2(0)]. As the length increases, *λ*LSA converges and the standard deviation

sufficiently becomes small; for example, the standard deviation for 5 μs is 0.025. In this work, we adopt 5 μs as the length of time series, and the results shown in the following figures are not means but values which are obtained by a single calculation. Then, in the following numerical simulations, some *λ*LSA diverge when *A*<sup>2</sup> → 0 or *A*<sup>0</sup> → 0 and are not shown in the following figures.

**Figure 2.** Mean of Lyapunov exponent plotted against length of time series for numerical simulation. The error bars represent standard deviations.

#### **3. Orbital Instability of Chaotic Laser Diode with Chaotic Applied Signal**

#### *3.1. Mean and Standard Deviation of Applied Signal*

In this section, we investigate the orbital instability of LD2 by applying a chaotic signal to the drive current of LD1. The chaotic signal is normalized to a value of [0, 1] using the parameter *a* in Equation (3). The Lyapunov exponent *λ*LSA is plotted against the optical injection ratio *r*inj in Figure 3. In Figure 3a, the black circles and red squares indicate *λ*LSA of LD2 without any applied signal and with an applied chaotic signal for *g* = 5.0 and *r*fb = 0.05, respectively. The mean and standard deviation of the applied chaotic signal are 0.113 and 0.080, respectively. The gray plots are the extrema of the intensity of LD2 without any applied signal, which shows the bifurcation diagram. When LD2 has no applied signal, for small *r*inj, LD2 oscillates stably or periodically and the corresponding *λ*LSA is nonpositive. With increasing *r*inj, the intensity of LD2 has a large number of extrema and the dynamics are chaotic, with positive *λ*LSA. Then, a large window is observed between the chaotic states around *r*inj ∼ 0.10 and small windows are observed for some other *r*inj, where LD2 oscillates periodically and the corresponding *λ*LSA is nonpositive. However, chaotic dynamics appear upon applying a chaotic signal to the drive current of LD1, and *λ*LSA > 0, as shown by the red squares in Figure 3a. This phenomenon is similar to that observed when by applying a pseudorandom signal in Ref. [24]. In addition, it seems that the red squares in Figure 3a shift slightly away from the black circles in the negative direction of *r*inj.

Next, we consider the applications of a pseudorandom signal with a mean of 0.113 and standard deviation of 0.080, which are the same values as those of the applied chaotic signal, and DC with *C*(*t*) = 0.113. In Figure 3a, the blue diamonds and purple triangles indicate *λ*LSA of LD2 with the pseudorandom signal and with the applied DC, respectively. Since the mean of the drive current increases by the applied signal for both plots, the chaotic dynamics of the injection system are enhanced and the plots are shifted away from the black circles in the negative direction of *r*inj. When the applied signal is a DC but not a pseudorandom signal, windows are observed.

In Figure 4, we confirm that the distribution of *λ*LSA is shifted away from the intrinsic distribution in the negative direction of *r*inj by applying a signal to the drive current of LD1, which is a DC and pseudorandom signal with a standard deviation of 0.10 in Figure 4a,b, respectively. The black circles, red squares, blue diamonds, and purple triangles indicate *g* · *C*(*t*) = 0, 0.10, 0.50 and 1.00, respectively. The gray plots show the bifurcation diagram for the intrinsic system without the applied signal. With increasing applied signal, the shift of the plots increases. Thus, the orbital instability is sensitive to the mean of the applied signal for the DC or pseudorandom signal.

However, the orbital instability is not always sensitive to the mean of the applied signal for a chaotic signal. In Figure 3b, we show *λ*LSA for the system with the chaotic applied signal of *r*fb = 0.20, where the mean and standard deviation of the signal are 0.108 and 0.102, respectively. The black circles, red squares, blue diamonds, and purple triangles indicate the injection system without any applied signal, with the chaotic signal, with the pseudorandom signal with a mean of 0.108, and standard deviation of 0.102, and with DC with a mean of 0.108, respectively. The plots for the system with the applied pseudorandom signal and the applied DC are shifted away from the black circles in the negative direction of *r*inj. Since the mean of the applied signal is larger than that in Figure 3a, the shift of the plots is larger. On the other hand, when the chaotic signal is applied, the plots are shifted away from the black circles in the positive direction of *r*inj. Therefore, it is considered that the factor contributing to the enhanced orbital instability is not the mean or standard deviation but another factor.

**Figure 3.** Bifurcation diagram and Lyapunov exponent plotted against injection ratio when (**a**) *r*fb = 0.05 and (**b**) *r*fb = 0.20. The black circles, red squares, blue diamonds, and purple triangles indicate the system without the applied signal, with the applied chaotic signal, with the pseudorandom signal, and with applied DC, respectively. The gray plots are the local maximal values of the intensity of LD2 without the applied signal.

**Figure 4.** Bifurcation diagram and Lyapunov exponent plotted against injection ratio with (**a**) applied DC and (**b**) applied pseudorandom signal. The black circles, red squares, blue diamonds, and purple triangles indicate *g* · *C*(*t*) = 0, 0.10, 0.50 and 1.00, respectively. The gray plots are the local maximal values of the intensity of LD2 without the applied signal.

#### *3.2. Optical Feedback Ratio and Optical Injection Ratio*

Here, to consider the effect of the applied chaotic signal, we show the orbital instability of the system with the applied chaotic signal as a function of the optical feedback ratio *r*fb of LD0 and the optical injection ratio *r*inj from LD1 to LD2. Figure 5 shows *λ*LSA plotted against *r*fb and *r*inj. According to the previous discussion, the mean of the applied signal may contribute to the orbital instability of LD2. Thus, all the applied signals in this subsection are normalized by the parameter *a* in Equation (3), and the mean of the applied signal is fixed as *g* · *C*(*t*) = 0.5 and 5.0, shown in Figure 5a,b, respectively. When the amplitude of the applied chaotic signal is small, *r*fb makes a small contribution to *λ*LSA for small *r*inj (Figure 5a). With increasing *r*inj, when *r*fb is large, *λ*LSA increases gradually. For larger *r*fb, *λ*LSA has a peak around *r*inj ∼ 0.14. On the other hand, when the amplitude of the applied chaotic signal is large, *r*fb makes a larger contribution to *λ*LSA (Figure 5b) than that in Figure 5a. The window around *r*inj ∼ 0.10, which is observed in the inherent system without the applied signal, is not observed, the peak around *r*inj ∼ 0.05 is shifted in the positive direction of *r*inj and the peak around *r*inj ∼ 0.14 gradually becomes large.

Figure 5c shows the maximal value of *λ*LSA in the range of 0 < *r*inj ≤ 0.20 and the corresponding optical injection ratio *r* inj plotted against *r*fb for *g* · *C*(*t*) = 0.5. When *r*fb ≤ 0.04, *λ*LSA has the maximal value at *r* inj ∼ 0.058. On the other hand, when *r*fb ≥ 0.04, *λ*LSA depends on *r*fb, and 0.10 ≤ *r* inj ≤ 0.15, where the window is observed in the inherent system without the applied signal. Similarly, we shows the maximal value of *λ*LSA in the range of 0 < *r*inj ≤ 0.20 and the corresponding optical injection ratio *r* inj plotted against *r*fb for *g* · *C*(*t*) = 5.0. When *r*fb ≤ 0.03, *λ*LSA has the maximal value at *r* inj ∼ 0.037. Then, the maximal value of *λ*LSA depends on *r*fb. Since the maximal value of *λ*LSA gradually increases and saturates with increasing *r*fb, it is controlled by *r*fb of the applied chaotic signal. In the next section, we discuss some quantities of the applied signal to study the conditions of the applied signal that cause large orbital instability of LD2.

**Figure 5.** Lyapunov exponent of the system with the applied chaotic signal as a function of the feedback ratio of LD0 and the injection ratio from LD1 to LD2 when (**a**) *g* · *C*(*t*) = 0.5 and (**b**) *g* · *C*(*t*) = 5.0, and maximal value of Lyapunov exponent in the range of 0 < *r*inj ≤ 0.20 and corresponding optical injection ratio *r* inj as a function of *r*fb when (**c**) *g* · *C*(*t*) = 0.5 and (**d**) *g* · *C*(*t*) = 5.0.

#### **4. Orbital Instability against Statistical and Dynamical Quantities of Applied Signal**

First, we study the applied chaotic signal, which is generated by LD0 with optical feedback. Figure 6 shows the extrema of the intensity of LD0 and the Lyapunov exponent *λ*LSA plotted against the optical feedback ratio *r*fb. Different symbols are used for different ranges of *r*fb, that is, purple down-pointing triangles, blue up-pointing triangles, green diamonds, orange squares, and red circles indicate *λ*LSA for 0 < *r*fb ≤ 0.040, 0.040 < *r*fb ≤ 0.080, 0.080 < *r*fb ≤ 0.120, 0.120 < *r*fb ≤ 0.160 and 0.160 < *r*fb ≤ 0.200, respectively. In the range of 0 < *r*fb ≤ 0.040, the fluctuation of the intensity is small and *λ*LSA is small. With increasing *r*fb, *λ*LSA inceases for 0.040 < *r*fb ≤ 0.120 and gradually decreases for 0.120 < *r*fb.

**Figure 6.** Bifurcation diagram and Lyapunov exponent of LD0. Purple down-pointing triangles, blue up-pointing triangles, green diamonds, orange squares, and red circles indicate *λ*LSA for 0 < *r*fb ≤ 0.040, 0.040 < *r*fb ≤ 0.080, 0.080 < *r*fb ≤ 0.120, 0.120 < *r*fb ≤ 0.160 and 0.160 < *r*fb ≤ 0.200, respectively. The gray plots are the local maximal values of the intensity of LD0.

Next, we study the orbital instability of LD2 with the applied chaotic signal for statistical and dynamical quantities of the applied chaotic signal to show the characteristics of the applied chaotic signal that can control the orbital instability of LD2. The maximal values of *λ*LSA of LD2 in the range of 0 < *r*inj ≤ 0.20 for certain *r*fb of LD0 are calculated and plotted against the standard deviation, skewness, kurtosis, Lyapunov exponent, bandwidth, and mode of the histogram of LD0 in Figure 7. The symbols correspond to those in Figure 6. In the range where the orbital instability of the applied signal is small (0 < *r*fb ≤ 0.040), the maximal value of *λ*LSA does not vary with *r*fb (purple down-pointing triangles in Figure 7). We then consider the range of 0.040 < *r*fb ≤ 0.200 where the orbital instability of the applied chaotic signal is sufficiently large. The correlation between the maximal value of *λ*LSA and the standard deviation of the applied signal is low in Figure 7a. On the other hand, in Figure 7b–d, the maximal value of *λ*LSA is nonlinear with the skewness, kurtosis, and Lyapunov exponent of the applied signal in the range of 0.040 < *r*fb ≤ 0.200, respectively. The plots for 0.040 < *r*fb ≤ 0.160 and 0.160 < *r*fb ≤ 0.200 have different gradients: thus, it is difficult to identify *λ*LSA from the statistical quantities. However, in Figure 7e–f, the maximal value of *λ*LSA is linear to the bandwidth and mode of the histogram of the applied signal in the range of 0.040 < *r*fb ≤ 0.200. Therefore, we can identify *λ*LSA from these quantities. Since the large bandwidth and small mode of the histogram of the applied signal contribute to the large Lyapunov exponent, the orbital instability of LD2 can be enhanced by applying a chaotic pulsation having a broad bandwidth.

**Figure 7.** Maximal value of Lyapunov exponent in the range of 0 < *r*inj ≤ 0.20 plotted against (**a**) standard deviation, (**b**) skewness, (**c**) kurtosis, (**d**) Lyapunov exponent, (**e**) bandwidth, and (**f**) mode of histogram of LD0.

Next, we discuss the shift of the maximal value of *λ*LSA upon applying the chaotic signal in Figure 5. We consider the optical injection ratio *r* inj0, which is *r* inj for the inherent system without the applied signal, and introduce the difference Δ*r* = *r* inj − *r* inj0. Figure 8 shows Δ*r* plotted against statistical and dynamical quantities of the applied chaotic signal, that is, the standard deviation, skewness, kurtosis, Lyapunov exponent, bandwidth, and mode of the histogram of LD0 as in Figure 7. In the range where the orbital instability of the applied signal is small (0 < *r*fb ≤ 0.040), Δ*r* is small for most plots. However, in the range of 0 < *r*fb ≤ 0.010, the orbital instability of LD2 is reduced since LD0 oscillates periodically or quasi-periodically. Since an additional optical injection is needed to obtain similar orbital instability, Δ*r* becomes large. In the range of 0.040 < *r*fb ≤ 0.160, the orbital instability of LD2 is enhanced in the range of 0.09 *r*inj - 0.13, where a window can be observed in the inherent system, and the plots are concentrated around Δ*r* ∼ 0.07.

The correlation between Δ*r* and the standard deviation, skewness and kurtosis of the applied signal is low in Figure 8a–c, respectively. In Figure 8e,f, Δ*r* is nonlinear to the bandwidth and mode of the histogram of the applied signal in the range of 0.040 < *r*fb ≤ 0.200, respectively. The plots for 0.040 < *r*fb ≤ 0.160 and 0.160 < *r*fb ≤ 0.200 have different gradients. On the other hand, Δ*r* seems to depend on the Lyapunov exponent of the applied signal in the range of 0.040 < *r*fb ≤ 0.200 (Figure 8d).

**Figure 8.** Optical injection ratio where the Lyapunov exponent is maximum in the range of 0 < *r*inj ≤ 0.20 plotted against (**a**) standard deviation, (**b**) skewness, (**c**) kurtosis, (**d**) Lyapunov exponent, (**e**) bandwidth, and (**f**) mode of histogram of LD0.

Finally, we discuss the chaos synchronization between LDs. For example, we assume the application of the present system to chaotic secure communication that is a digital scheme by using a difference of the orbital instability of chaotic LD [18]. The scheme is hardware-dependent, where the key to communication is based on the parameter of the LD system. LD1 and LD2 act as the transmitter and receiver LDs, respectively. Then, LD0 is the driver used to control the dynamics of LD1 and the message is modulated by LD0 and applied to LD1. The orbital instability of LD2 is controlled by LD0 through LD1 and corresponds to the digit. The proper receiver quantifies from only the optical intensity of LD2 at a certain interval, for example, using the method in Ref. [32], and compares the quantified orbital instability with the predetermined threshold to decide the digit. Since the dynamics of LD2 are decided by the parameters of three LDs, it is difficult for eavesdroppers to decode the digit with only the transmitting signal. However, if LD2 synchronizes with the other LDs, the eavesdropper can estimate the digit from the transmitting signal. Thus, we investigate the chaos synchronization between LD2 and the other LDs. In Figure 9, we calculate the correlation coefficient between the LDs plotted against *r*fb and *r*inj to quantify the chaos synchronization. Figure 9a,b show the correlation

coefficients between LD0 and LD2 and between LD1 and LD2, respectively. The correlation coefficient is expressed as

$$\rho\_{l2} = \max\_{\Delta t} \frac{\left\langle \left( I\_i(t) - \langle I\_i \rangle \right) \left( I\_2(t + \Delta t) - \langle I\_2 \rangle \right) \right\rangle}{\sqrt{\left\langle \left( I\_i(t) - \langle I\_i \rangle \right)^2 \right\rangle \left\langle \left( I\_2(t + \Delta t) - \langle I\_2 \rangle \right)^2 \right\rangle}},\tag{13}$$

where *Ii* and *I*<sup>2</sup> indicate the optical outputs of LD*i* (*i* = 0, 1) and LD2, respectively, and · indicates the ensemble average. The roundtrip time of the external cavity *τ*fb and the trip time of the injection light from LD1 to LD2 *τ*inj, have the same value, and the correlation coefficient is calculated in the range −10*τ*fb ≤ Δ*t* ≤ 10*τ*fb. Figure 9c,d show the maximal value of *ρ*<sup>02</sup> and *ρ*<sup>12</sup> in the range of 0 < *r*inj ≤ 0.20 and the corresponding optical injection ratio *r* inj plotted against *r*fb. In Figure 9a,c, the correlation coefficient *ρ*<sup>02</sup> between LD0 and LD2 is small in the entire range and the maximum is 0.11, showing that LD2 does not synchronize with LD0. On the other hand, as shown in Figure 9b,d, the correlation coefficient *ρ*<sup>12</sup> between LD1 and LD2 is larger than that in Figure 9a. For all *r*fb, the correlation coefficient is small in the range of *r*inj ≤ 0.05. With increasing *r*inj, the correlation coefficient becomes larger since the orbital instability of LD2 is enhanced (*r*inj ∼ 0.05). With further increase of *r*inj, the correlation coefficient becomes small again in the range of *r*inj ≥ 0.10, where a window can be observed in the inherent system. Since the maximal correlation coefficient between LD1 and LD2 is 0.40, chaos synchronization between LD1 and LD2 is not achieved. Therefore, it is concluded that the orbital instability of LD2 can be controlled by varying the parameters of LD0, which generates the applied chaotic signal, without chaos synchronization between the LDs.

**Figure 9.** Correlation coefficient as a function of the feedback ratio of LD0 and the injection ratio from LD1 to LD2: (**a**) between LD0 and LD2 and (**b**) between LD1 and LD2, and maximal value of correlation coefficient in the range of 0 < *r*inj ≤ 0.20 and corresponding optical injection ratio *r* inj as a function of *r*fb: (**c**) between LD0 and LD2 and (**d**) between LD1 and LD2.

#### **5. Conclusions**

We numerically studied the orbital instability of a chaotic laser diode (LD) system with optical injection, which consists of the master LD (LD1) and slave LD (LD2). The drive current of LD1 is modulated by the chaotic applied signal, which is generated by LD0 with optical feedback. First, we showed that chaotic behavior in the window is actualized by applying the chaotic signal as well as a pseudorandom signal. The optical injection ratio required to oscillate LD2 chaotically is decreased by applying the pseudorandom signal or DC but increased by applying the chaotic signal.

Next, we investigated the maximal value of the Lyapunov exponent of LD2 in the range of 0 < *r*inj ≤ 0.20 as a function of the optical feedback ratio of LD0 and the optical injection ratio from LD1 to LD2. When the amplitude of the applied chaotic signal is sufficiently large and the inherent system without the applied chaotic signal is in the window, the Lyapunov exponent of LD2 can be controlled by varying the optical feedback ratio of LD0.

Then, we discussed the effect of statistical and dynamical quantities of the applied chaotic signal on the orbital instability of LD2. The bandwidth and mode of the histogram of the applied chaotic signal are linear to the maximal value of the Lyapunov exponent of LD2. It was shown that the orbital instability of LD2 can be enhanced efficiently by applying a chaotic pulsation having a broad bandwidth.

Finally, we investigate the chaos synchronization between LDs. The LDs do not synchronize with each other. It was shown that the orbital instability of the chaotic LD can be controlled without chaos synchronization. Since it is difficult to estimate the dynamics of LD0 from the optical intensity of LD1, the characteristics is useful to the application of chaotic LD like a secure communication.

**Author Contributions:** Conceptualization, S.E.; investigation, S.E.; supervision, S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Uchida Energy Science Promotion Foundation (Grant No. 1-1-38).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Jacobian Matrix**

The Jacobian matrix on the right side of Equation (10) is defined as

$$J\_{\rm inj} = \begin{pmatrix} \frac{\partial f\_{A2}}{\partial A\_2} & \frac{\partial f\_{A2}}{\partial \phi\_2} & \frac{\partial f\_{A2}}{\partial n\_2} \\ & \frac{\partial f\_{\phi 2}}{\partial A\_2} & \frac{\partial f\_{\phi 2}}{\partial \phi\_2} & \frac{\partial f\_{\phi 2}}{\partial n\_2} \\ & \frac{\partial f\_{n2}}{\partial A\_2} & \frac{\partial f\_{n2}}{\partial \phi\_2} & \frac{\partial f\_{n2}}{\partial n\_2} \end{pmatrix},\tag{A1}$$

where *fA*2, *fφ*2, and *fn*<sup>2</sup> indicate the function for the right side of Equations (4)–(5), respectively. Then, the matrix used in our work described by

$$\mathbf{J}\_{\rm inj} = \begin{pmatrix} \frac{1}{2} G\_N n\_2(t) & -\kappa\_{\rm inj} A\_1(t - \tau\_{\rm inj}) \mathcal{S}\_{\rm linj}(t) & \frac{1}{2} G\_N A\_2(t) \\\\ \kappa\_{\rm inj} \frac{A\_1(t - \tau\_{\rm inj})}{A\_2^2(t)} \mathcal{S}\_{\rm linj}(t) & -\kappa\_{\rm inj} \frac{A\_1(t - \tau\_{\rm inj})}{A\_2(t)} \mathcal{C}\_{\rm linj}(t) & \frac{1}{2} a G\_N \\\\ -2[\Gamma + G\_N n\_2(t)] A\_2(t) & 0 & -\gamma - G\_N A\_2^2(t) \end{pmatrix} . \tag{A2}$$

Here, Sinj(*t*) = sin[*ωτ*inj + *φ*2(*t*) − *φ*1(*t* − *τ*inj)] and Cinj(*t*) = cos[*ωτ*inj + *φ*2(*t*) − *φ*1(*t* − *τ*inj)]. Similarly, since the matrix on the right side of Equation (11) is defined as

$$J\_{\mathbf{fb}} = \begin{pmatrix} \frac{\partial f\_{A0}}{\partial A\_0} & \frac{\partial f\_{A0}}{\partial \phi\_0} & \frac{\partial f\_{A0}}{\partial n\_0} & \frac{\partial f\_{A0}}{\partial A\_{\mathbf{fb}}} & \frac{\partial f\_{A0}}{\partial \phi\_{\mathbf{fb}}}\\ \frac{\partial f\_{\phi0}}{\partial A\_0} & \frac{\partial f\_{\phi0}}{\partial \phi\_0} & \frac{\partial f\_{\phi0}}{\partial n\_0} & \frac{\partial f\_{\phi0}}{\partial A\_{\mathbf{fb}}} & \frac{\partial f\_{\phi0}}{\partial \phi\_{\mathbf{fb}}}\\ \frac{\partial f\_{n0}}{\partial A\_0} & \frac{\partial f\_{n0}}{\partial \phi\_0} & \frac{\partial f\_{n0}}{\partial n\_0} & \frac{\partial f\_{n0}}{\partial A\_{\mathbf{fb}}} & \frac{\partial f\_{n0}}{\partial \phi\_{\mathbf{fb}}} \end{pmatrix},\tag{A3}$$

where *fA*0, *fφ*<sup>0</sup> and *fn*<sup>0</sup> indicate the function for the right side of Equations (7)–(9), respectively, the matrix used in our work described by

$$\begin{pmatrix} \begin{aligned} &=\\ & \begin{pmatrix} \frac{1}{2}\mathsf{G}\_{\mathrm{N}}u\_{0}(t) & -\mathsf{x}\_{\mathrm{fb}}A\_{0}(t-\mathsf{z}\_{\mathrm{fb}})\mathsf{S}\_{\mathrm{lb}}(t) & \frac{1}{2}\mathsf{G}\_{\mathrm{N}}A\_{0}(t) & \mathsf{x}\_{\mathrm{fb}}\mathsf{C}\_{\mathrm{lb}}(t) & \mathsf{x}\_{\mathrm{fb}}A\_{0}(t-\mathsf{z}\_{\mathrm{fb}})\mathsf{S}\_{\mathrm{lb}}(t) \\\\ & \mathsf{x}\_{\mathrm{fb}}\frac{A\_{0}(t-\mathsf{z}\_{\mathrm{fb}})}{A\_{0}^{2}(t)}\mathsf{S}\_{\mathrm{lb}}(t) & -\mathsf{x}\_{\mathrm{fb}}\frac{A\_{0}(t-\mathsf{z}\_{\mathrm{fb}})}{A\_{0}(t)}\mathsf{C}\_{\mathrm{lb}}(t) & \frac{1}{2}\mathsf{a}\mathsf{G}\_{\mathrm{N}} & -\frac{\mathsf{x}\_{\mathrm{fb}}}{A\_{0}(t)}\mathsf{S}\_{\mathrm{lb}}(t) & \mathsf{x}\_{\mathrm{fb}}\frac{A\_{0}(t-\mathsf{z}\_{\mathrm{fb}})}{A\_{0}(t)}\mathsf{C}\_{\mathrm{lb}}(t) \\\\ & -2[\Gamma + \mathsf{G}\_{\mathrm{N}}u\_{0}(t)]A\_{0}(t) & 0 & -\gamma - \mathsf{G}\_{\mathrm{N}}A\_{0}^{2}(t) & 0 & 0 \end{aligned} \end{cases} \tag{A4}$$

Here, Sfb(*t*) = sin[*ωτ*fb + *φ*0(*t*) − *φ*0(*t* − *τ*fb)] and Cfb(*t*) = cos[*ωτ*fb + *φ*0(*t*) − *φ*0(*t* − *τ*fb)].

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
