1. Introduction

Complex instabilities and dynamics are common in many different laser systems. Some examples are listed. CO₂ lasers show self pulsing and deterministic chaos when the laser is subject to feedback [1]. Ammonia (NH₃) ring lasers have been shown to produce similar instabilities and period-doubling routes to chaos have been investigated in this system [2]. Semiconductor lasers can display complex, nonlinear dynamics when subject to direct current modulation [3], optical feedback [4], strong optical injection from a master laser [5], optoelectronic feedback [6] or a hybrid combination of these [7, 8]. In many instances a perturbation such as feedback or injection, which may cause complex dynamics for some parameter values, can also be used to stabilize the laser with different sets of values [9, 10].

Recently, complex laser systems have generated interest, not only for fundamental nonlinear dynamics research, but also for their potential application as chaotic sources for secure communication schemes. This is a concept that was first proposed in nonlinear electronic circuits [11]. The potential to obtain much higher frequency bandwidths in chaos using optical signals rather than electronics, has prompted development of systems that display complex instabilities, such as Erbium doped fiber ring lasers [12], solid state laser microchip lasers (Nd:YVO₄) [13], and semiconductor lasers [14]. These have been specifically researched for applications in secure communications. Schemes incorporating chaotically varying wavelength have also been proposed [15]. Systems employing semiconductor lasers with direct optical feedback and optoelectronic feedback to produce chaotic output power fluctuations have been used successfully to transmit and receive chaotically masked signals over commercial fiber-optic networks [16].

Successful development of nonlinear systems requires a thorough understanding of the dynamics generated. Dynamical maps are a good way to identify different dynamical regions within a parameter space of interest. Previous studies have used optical or radio frequency (RF) spectra to generate maps for semiconductor lasers with optical injection [17–20]. Experimental time-series from optically injected solid state lasers have also been analyzed to generate dynamical maps based on intensity [21] or bifurcation analysis [22].

Nonlinear time-series and chaos analysis techniques can be utilized to gain further insight to these systems and have been used to quantify the observed complex dynamics of real world systems from many areas of research [23]. Maps of the largest Lyapunov exponent, a nonlinear time-series measure of the sensitivity of a system to initial conditions [24, 25], have been generated from theoretically modeled optically injected semiconductor lasers for comparison with experimental results [26], for example.

Another reason chaotic and dynamically unstable lasers are useful is that the experimenter has much more control over system variables than is possible in most other nonlinear systems, such as many biological, physiological, mechanical, hydrodynamical and electronic nonlinear systems. For example, employing nonlinear time-series analysis, and having the ability to systematically vary the signal-to-noise of a photodetected chaotic semiconductor laser signal, has allowed the transition from deterministic chaos, through chaos plus noise, to a noise dominated signal, and allowed the levels of acceptable noise that would not mask the chaotic dynamics to be identified [27]. In contrast, most biological and physiological nonlinear systems do not permit independent control of signal-to-noise. Study of noise masking effects in laser-based systems can provide insight that helps interpret nonlinear data analysis in these systems.

Nonlinear laser systems are also ideal for developing and testing the nonlinear analysis concepts themselves. This has been suggested in [28], where they use nonlinear time-series techniques to analyze the chaotic output of a Nd:YAG laser with an intra-cavity second harmonic generating crystal, in order to further understand the impact noise has on the laser output. The high degree of control that one has over time-series length, signal-to-noise level and sampling rate that is experimentally achievable in laser systems make them ideal for refining robust analysis techniques. These procedures can then be implemented for use with other complex systems where short, noisy and/or under sampled data sets are unavoidable.

Despite the vast number of publications on nonlinear time-series analysis, there is still confusion and conflicting reports on whether the results of such analysis are physically meaningful. One of the main problems with implementing such analysis techniques is that the results generally require subjective interpretation by the user. For instance, calculation of the correlation dimension [29], an appropriate analysis measurand for nonlinear time-series, generates a graph of the correlation sum (explained in detail later) from which a region of constant slope needs to be identified. In most cases this identification is done by visual inspection of the graph. Consequently, discrepancies can arise due to varying criteria implemented by the user. Manual interpretation of such results from a complex system over a large parameter range, such as is required to produce a full dynamical map, is impractical. Thus, an automatic computational system able to generate dynamical maps is a tool which will significantly facilitate analysis of the complexity of nonlinear dynamical systems. An important part of any such automated computation is development of appropriate protocols to be incorporated to ensure physically meaningful dynamical measures are calculated.

In this work such an automated computational process is developed and applied to the analysis of the output power time-series from an optically injected Nd:YVO₄ solid state laser system [21]. This laser system generates nonlinear dynamics with levels of complexity that can be robustly measured using a standard chaos data analysis nonlinear time series algorithm [29]. The output power time-series data from the optically injected Nd:YVO4 solid state laser (OISSL) has been characterized in previous investigations. Maps of maximum intensity [21], and bifurcation diagrams [22], have been generated as a function of the injection strength from the master laser into the slave laser, K , and the angular frequency detuning between the free-running master laser and free-running slave laser Δω. The master/slave OISSL system is described in full detail in [21], and the theoretical model used for generating simulated output power time series for the system is described in [5, 30]. In the research reported herein it is the correlation dimension that is calculated, from both experimental and simulated data. This is a good measure of the level of complexity of the dynamics of a system and it is well suited to being extracted from experimental time-series data of only one system variable [29].

In this paper, an overview of nonlinear time series analysis techniques, and the correlation dimension in particular, is given in Section 2. The details of the experimental method, automated analysis program and theoretical model, along with main results are given in Section 3. Discussion of the results and comparisons between the experimental and simulated data is given in Section 4. The variety of dynamics identified by the automated analysis program demonstrate how this is a useful tool for investigating and understanding nonlinear systems.

2. Nonlinear time-series analysis

In general, nonlinear systems are best described and analyzed by means of a phase space trajectory or attractor, which represents all the possible states of a particular system. The dimension of the phase space required to completely describe a system needs to be equal to the number of variables describing the system, with each axis representing a single system variable. Each point on the trajectory describes a unique state of the system. Such a trajectory could be constructed from observations made on experimental nonlinear systems if it were possible to take measurements of each variable in time. In most cases this is not possible due to experimental limitations. Typically, measurements of only one system variable in time are available. In many cases even the actual number of system variables for a particular nonlinear system is unknown.

It has been shown that a ‘reconstructed’ attractor can be formed using the embedding method [31], which has all the same properties as the ‘true’ attractor. The method involves taking a single scalar measurement, in this case the intensity of the laser output power as a function of time P(T), and plotting it against a time delayed version of itself P(t+T). Since the dimension of the attractor is not known beforehand, we aim to embed the attractor sequentially in a phase space with dimension m, equal to and greater than the dimension of the attractor. This ensures that the orbit of the attractor is completely unfolded and no overlaps due to projection occur. The vectors describing the attractor are given by,

The value of time delay T and its effect on the phase space reconstruction has been previously investigated [32]. Constraints on the choice of T are that it must be a multiple of the sampling time and also ensure that data is separated as much as possible without the points becoming completely independent. While there is no rigorous way of determining the optimal value of T, a good approximation, and the method used in this work, is to use the first minimum of the mutual information function [33]. There have been many other techniques suggested, see for example [34, 35], however none have been shown to perform significantly better than others and all give values of similar magnitude [23].

The nonlinear time-series measurand used in this work, the correlation dimension, can be calculated using the Grassberger-Procaccia algorithm [29]. This involves calculation of the correlation sum defined as

where ∥y_i − y_j∥ is the distance between pairs of points in the phase space and H(x) is the Heaviside function (H = 1 if x ≥ 0 and H = 0 if x < 0). The value r is the radius of a hypersphere in an m-dimensional space, and N is the number of points in the space.

In theory the density of the points will scale as a power law with the radius

The correlation dimension is then given by

So the gradient of a graph of log C_m(r) versus log r should approach a finite value D in the limit as r → 0. For data that has finite sampling, the slope of this graph for very small r is very inaccurate so instead a ‘scaling” region of r is taken over which the slope is relatively stable. A plot of the gradient, D_m(r) = ∂ log C_m(r)/∂ log r, as a function of r makes observing this region much clearer. The value of the slope D_m(r) within the scaling region is the correlation dimension.

The graphs in Fig. 1 show the correlation sum, C_m(r), and corresponding slopes, D_m(r), as a function of radius. The data used is the x-variable time-series (10 000 points) of the chaotic Lorenz equations [36]. Looking between radius values of 3 and 10.8, the curves of D_m(r) can be seen to saturate to a fixed level as the data is embedded in higher and higher dimensions m. This ‘plateau’ is the scaling region and, if present, indicates that the data is deterministic. In the case of this Lorenz data we can infer that the minimum embedding dimension required to unfold the attractor is 8, and that the value of the correlation dimension is slightly larger than 2 (the value of the slope in the scaling region). If the time-series was from a purely stochastic process it would give rise to a slope that increases continually for higher embedding dimensions. Thus, the absence of a scaling region indicates that the data could be stochastic in nature, or severely affected by noise.

 

Fig. 1. The graph in (a) shows correlation sum as a function of radius and (b) shows the gradient of the correlation sum as a function of radius for values of the embedding dimension m between 2 and 12. Time series analyzed is the x-variable data from the chaotic Lorenz equations [36].

Download Full Size | PDF

The value of the correlation dimension characterizes the shape of a chaotic attractor [29]. A value of 1 is representative of periodic or limit-cycle dynamics and a value of 2 arises from quasi-periodic data that generates a torus in the phase space. A higher fractional value of correlation dimension indicates the phase space trajectory has a higher level of complexity, arising from a chaotic system, where the attractor is bounded within a region of phase space but does not have an exact geometrical shape. A purely stochastic signal will have no structure and will completely fill the phase space, leading to an infinitely increasing value of D_m(r).

When performing this type of analysis on experimental data care must be taken in the interpretation of the results. It has been shown that naïve implementation of the Grassberger and Procaccia algorithm [29] can lead to serious underestimation of dimension value [37]. The two most important factors to take into consideration are (i) the removal of temporally correlated points from the correlation sum calculation [38], and (ii) the requirement that a comparison of results with surrogate data sets must also be performed [39].

In order to achieve the first of these, the correlation sum (Eq. 2) needs to be slightly modified so that points in the phase space that were close in time are omitted from the calculation by the introduction of a Theiler window [38]. The modified equation is

where w is the value of the Theiler window. This omission is especially important when dealing with flow-like data since pairs of points that are close in time are generally close in phase space as well. These temporally correlated pairs introduce a bias in the correlation sum calculation. Different authors suggest different methods of determining the size of w. A mathematically rigorous approach is to use a space-time separation plot which allows the minimum value of w to be determined by looking at the number of pairs of points as a function of both time separation and spatial distance [40]. In this work, since there are a large number of data points and the dynamics are highly sampled, the size of the Theiler window was set generously at 100 points to virtually eliminate the use of temporally correlated points in the correlation sum. This is likely an over estimate, since the correlation times of all of the time series are less than 10 data points (determined as the number of points for the autocorrelation function to drop to 1/e). The effect of using a value larger than required is a small loss in statistics (about 2w/N) due to using fewer points in the correlation sum [23]. This can cause an increase in uncertainty in the value of the correlation dimension. However, in this case the loss of 4% of “point pairs” has no significant impact on the results. Omission of the space time separation plot from our automated system results in a simpler, more computationally efficient process. It was not required in this instance because there was sufficient data density captured from the OISSL system. The “space-time separation” step could be included in the automated computations for time-series data sets where the time series are short and where loss of data points needs to be minimized.

It has also been shown that the correlation dimension analysis can give spuriously low dimension values for linearly filtered random numbers [37]. Thus, a method is required to distinguish between linearly correlated noise and nonlinear structure. The use of surrogate data solves this problem and is regarded as a necessary tool for validating the results of nonlinear time-series analysis techniques [39]. The idea is that a null hypothesis is specified which then needs to be tested. In this case, the hypothesis would be that the data set being analyzed is from a linearly filtered stochastic process. If this null hypothesis can be rejected then we can have confidence that the data is truly from a deterministic nonlinear process.

The method of discriminating truly nonlinear dynamics involves performing the dynamical analysis, in this case evaluation of the correlation dimension, on the original set of data to obtain a result CD_orig. Surrogate data sets are then created from the original set, and the same dimensional analysis is performed on these. The average value of the correlation dimension obtained from the surrogates, <CD_surr>, is then compared to CD_orig to see if they are significantly different.

The surrogates data sets themselves can be formed in many ways [41]. In this work we use random phase surrogates. These are generated by taking the Fourier transform of the original data set, randomizing the phase, and performing the inverse Fourier transform. This results in the surrogates having the same power spectrum and autocorrelation function as the original data. If the results are significantly different then the null hypothesis is rejected and we can have confidence in the correlation dimension value for the original data set. The issue of how many surrogates to use has no simple solution. Statistically, the more surrogates used the better, however the fact that each measurement takes a significant amount of time means one has to balance the accuracy of the analysis with the time required to perform it.

3. Experiment and analysis

The experimental system outlined in Section 3.1 is fully described in [21, 42]. The automated computational system used to analyze this data is described in Section 3.2 and is based on the correlation dimension algorithm [29]. The model for this laser system used to generate simulated data is discussed in Section 3.3 and fully described in [5, 30].

3.1 Laser system

The experimental setup of the optically injected Nd:YVO₄ solid state laser system can be seen in Fig. 2. Two 1064 nm solid state Nd:YVO₄ lasers were operated in a master-slave configuration. The Nd:YVO₄ crystals (CASIX) are 1 mm thick with 1% Nd³⁺ atomic doping. The laser cavity was formed by HR coating on the front surface and 5% output coupling at the end facet. The laser crystals were pumped with 150 mW diode lasers (SDL) operating at 809 nm. The pump power was low enough to ensure single mode operation. Both the pump and laser polarization were along the c-axis of the crystal. Faraday isolators were used to block feedback to the pump lasers. The master laser crystal mount could be temperature controlled using a Peltier element. Interference filters removed excess pump light from the 1064 nm beams. Faraday isolators made the coupling between master and slave unidirectional and prevented unwanted feedback. Light from the master laser was injected into the slave via an acousto-optic modulator (AOM) and the relative injection power measured on the detector PD1. The beat frequency between the two lasers was measured on a 400 MHz detector (PD3). The optical input (PD2) of a Tektronix CSA 7404 oscilloscope was used to measure the intensity time-series of the slave laser.

 

Fig. 2. Experimental setup: LD, laser diode; FI, Faraday isolator; TEC, temperature control; IF, interference filter; BS, beam splitter; AOM, acousto-optic modulator; FP, Fabry-Pérot interferometer; PD, photodetectors. Figure reproduced from [21].

Download Full Size | PDF

The power of the injected beam was modulated using the AOM and the temperature of the master laser was tuned with the Peltier element to control the frequency detuning between the lasers (Δω). The maximum range over which the frequency is shifted is approximately 32 MHz. This is very narrow compared to the mode spacing so no mode hopping occurs.

For a fixed detuning, long time-series data of the laser intensity were recorded as the injection strength (K) was slowly increased from zero until injection locking was achieved. This injection sweep was repeated at a large number of fixed frequency detuning values (Δω).

3.2 Automated Computational Process

An automated computation system was implemented so that the correlation dimension could be mapped at every point in the (K, Δω) plane. The algorithm was coded using a combination of both Matlab and C++ programming languages.

The processing steps involved in the automated computation are shown in Fig. 3. The raw experimental data consists of 315 slave laser intensity time-series, each corresponding to a normalized frequency detuning value Δω between −4.11 and 3.25 units (a multiplier of the angular relaxation oscillation frequency). Each of these contains 883 000 data points of the laser intensity, sampled at 10 ns intervals, for a complete sweep of the injection strength from zero until the level at which the system reaches frequency locking of the slave to the master laser. The necessary injection strength to achieve frequency locking increases linearly with the frequency detuning between the master laser and the free running slave laser. Each of these complete time-series contains the whole range of dynamics possible for that value of frequency detuning. An example of one of these time-series is shown in Fig. 4.

The first step in the automated process was to load a complete time-series for one fixed detuning value. This was then divided up into 176 smaller subsets of 5000 points or 50 μs. Since the sweep time of the injection strength is on a much longer scale (approx. 9 ms), each subset corresponds to practically constant injection strength. The largest change in injection strength for a time-series-subset occurs at maximum detuning. Here the injection strength sweeps from K = 0 to 3.9, resulting in a maximum change of ΔK = 0.03 within the 50 μs subset. The subsets also contain at least a hundred cycles of the prevalent dynamics. This was determined to be enough for robust analysis as very similar results were obtained when analyzing the same data sets but using only the first 2000 data points.

 

Fig. 3. Structure diagram of the process of the automated computational system.

Download Full Size | PDF

Each subset was then analyzed individually. First, the time lag T required to reconstruct the phase space trajectory was calculated by locating the first minimum of the mutual information function [33]. Typically, values for this time lag were less than 15 data points. An attractor was then mapped out for embedding dimensions m = 2–10 by using Eq. (1). Each attractor was analyzed using the modified correlation sum in Eq. (5), from which a graph of the gradient as a function of radius, D_m(r), was generated.

 

Fig. 4. Experimental laser intensity time-series trace and corresponding injection strength sweep (inset) for detuning Δω = −1.0.

Download Full Size | PDF

Once D_m(r) curves for all embeddings were constructed, the scaling region, if present, was located by using a modified ‘most probable dimension value’ method [43]. This involves generating a histogram of how many times each gradient value is obtained, in all embeddings, and locating the peak value. The width of the peak in the histogram is an appropriate way to quantify confidence in the scaling region. A narrow peak indicates a significant plateau whereas a wider peak results from a plateau that is not as well defined, as seen in Fig. 5.

 

Fig. 5. Plots of D_m(r) and corresponding histograms for examples of (a) ‘good’ and (b) ‘bad’ scaling regions.

Download Full Size | PDF

If a scaling region was present then the value of D_m(r) at the histogram peak was recorded as the correlation dimension. If no plateau was detected then no value was recorded for correlation dimension. The width of this scaling region peak, measured at half the maximum, is a measure of the uncertainty in the correlation dimension. The correlation dimension values were then mapped into the corresponding positions in the (K, Δω) plane. The final dynamical map is shown in Fig. 6.

The final step in the automated system repeats the whole process shown in Fig. 3 several times, each time replacing the 50 μs time-series with a different surrogate data set. The surrogates were created by taking the Fast Fourier Transform (FFT) of the original 50 μs time-series, randomizing the phases and taking the inverse FFT. This gives a surrogate time-series that has the same spectral content as the original time-series, but with the determinism (if any) removed. A total of 10 surrogate data maps were averaged to produce a map for comparison with the original experimental map in Fig. 6. The different dynamical regions identified as I, Locked; II, Periodic; III, ‘Spiky’ Output; IV, Chaotic and V, Noisy, are discussed in Section 4.

3.3 Theoretical Model of Solid State Laser Subject to Optical Injection

A theoretical model of the specific Nd:YVO₄ solid-state laser subject to optical injection has been described in [5, 30], and is reproduced here:

Here a is the amplitude of the slowly varying field envelope and n is the population inversion density. Both are normalized to their steady state values. Ĵ = (J − J^th)/ J_th is the pump power J, normalized to the threshold value J_th. The linewidth enhancement factor, α, couples the gain variations to the refractive index and is set between 0.2 and 0.35 in this model to match experimental observations. The decay rates for the cavity, and for the upper laser level are γ_c and γ_s, respectively; γ_n and γ_p are the relaxation rates for the differential and nonlinear gain, respectively. The experimentally controllable injection parameters κ and Ω. stand for the coupling of the injected field and the angular frequency detuning between the master and slave lasers. The dimensionless normalized parameters referred to in the experimental section are related by K = κ/Ω_R for injection strength, and Δ_ω = Ω/Ω_R for angular frequency detuning, where Ω_R = 2πf_R is the angular relaxation oscillation frequency. Noise can be introduced through the Gaussian white noise source terms F_a and F_n. For the simulations in this paper no noise contribution is taken into account i.e. F_a = F_n = 0.

Table 1. Parameter values used in the simulations^*

View Table

Using the parameter values shown in Table 1 (reproduced from [21]), simulated time-series were generated, and the same automated correlation dimension analysis was performed. The correlation dimension map is shown in Fig. 7.

 

Fig. 6. Correlation dimension map in the (K, Δω) plane for experimental data. Different dynamical region are identified as I: Locked (no data), II: Periodic (CD = 1), III: ‘Spiky’ Output (CD < 1), IV: Chaotic (CD > 2), V: Noisy (CD = ∞).

Download Full Size | PDF

 

Fig. 7. Correlation dimension map in the (K, Δω) plane for simulated data without noise.

Download Full Size | PDF

3.4 Dynamical Maps Using Automated Correlation Dimension Analysis

The automated correlation dimension map method is computationally intensive but does allow maps to be generated that would be practically impossible if they were to be done manually. The actual computation time varies with data set size, maximum embedding dimension and number of time-series analyzed. For example, in this work, using a PC with modest specifications (Pentium 4 3.6 GHz CPU, 3 GB RAM), with data set size of 5000 points, maximum embedding dimension of 10, and 55 440 time-series analyzed, the generation of a complete map took approximately 315 hours to compute. Parallel processing using several PC’s meant the ‘actual’ time to generate a map was only a fraction of this. The same map to be manually constructed would require a person to view 55 440 individual graphs and visually locate plateau regions to record the dimension. Quicker computation could be effected using a more powerful computer and/or carrying out further optimization of the algorithm coding.

Previous analysis of the experimental data to explore the nature of the dynamics of the optically injected Nd:YVO₄ solid state laser system have used peak intensity amplitude [21] and bifurcation analysis [22] to generate dynamical maps. The maps from this prior work are reproduced here in Fig. 8(a) and (b) to compare with the correlation dimension map in Fig. 6. The regions identified in the bifurcation diagram in Fig. 8(b) mirror those seen in the correlation dimension map in Fig. 6.

 

Fig. 8. (a) Experimental (K, Δω) map of intensity time-series maxima. Reproduced from [21]. (b) Theoretical bifurcation diagram in the (K, Δω) plane. Showing Hopf (H) bifurcations in red, saddle-node (SN) bifurcations in blue, torus (T) bifurcations in black and period doubling (PD) bifurcations in green. Reproduced from [22].

Download Full Size | PDF

4. Discussion

The map in Fig. 6 shows the value of the correlation dimension at every point in the (K, Δω) plane, as determined by the automated process for the experimental OISSL data. The white triangular region in the middle, identified as Region I in Fig. 6, is where the laser becomes frequency locked. Here the laser output is single frequency and constant in power (with small amplitude relaxation oscillations). Several other distinct regions can be identified within this map, and they have been labeled as Regions II, III, IV and V in Fig. 6.

Figure 9 shows samples of the output power time-series, plots of gradient D_m(r) and the corresponding histograms for each of the different dynamical regions (II-V) in the (K, Δω) plane. Figure 9(a) is a typical periodic signal seen within the light blue (CD = 1) Region II of the experimental CD map. Fig. 9(b) represents the ‘spiky’ output generally seen in Region III, the dark blue (low CD) region of the CD map. Fig. 9(c) shows a ‘chaotic’ output from Region IV, identified by the red (high CD) points in the CD map. Figure 9(d) is representative of Region V, corresponding to noisy time-series that gave no scaling region in the correlation sum and therefore an infinite correlation dimension.

 

Fig. 9. Typical output power time-series, graphs of Dm(r) versus r and corresponding histograms for regions (a) II, (b) III, (c) IV and (d) V in the (K, Δω) plane.

Download Full Size | PDF

The automated computation correctly assigned a correlation dimension value of ~1, indicating a limit cycle in the phase space, to the periodic time-series seen in Region II.

In Region III the time-series are characterized by large amplitude spikes, which give rise to values for correlation dimension of 0.1 – 0.3. Within this region, the CD analysis appears to differ depending on the value of the embedding delay T used. When the program detects the first minimum of the mutual information as 1 data point, the D_m(r) graph and histogram register a narrow peak as seen in Fig. 10(a). However, it is also observed that when the delay is longer than 1 data point the D_m(r) graph obtained is similar to that seen in Fig. 10(b).

Clearly, a longer delay results in an incorrect detection of a scaling region, as in Fig. 10(b). This can be detected by looking at the width of the peak which is larger when this erroneous detection occurs. A map showing the value of the width of the scaling region peak is shown in Fig. 11. From this map it can be concluded that many instances of this are occurring within Region III as the scaling region peak width in this region is generally larger than in the other regions identified. In this case we cannot claim a correlation dimension was obtained since no true scaling region is observed. However, even though the values obtained in this region are not a true CD, the automated process does generate a low CD estimate (with a wide peak) reliably which can be identified as a consistent signature to identify the type of ‘spiky’ dynamic that occurs from the OISSL system for parameter values of Region III.

 

Fig. 10. Typical output power time-series, graphs of D_m(r) versus r and corresponding histograms (showing scaling region peak width) for data from region III with values of embedding delay (a) T = 1 data point, and (b) T = 7 data points.

Download Full Size | PDF

 

Fig. 11. Map of the width of scaling region peaks for experimental OISSL data in the (K, Δω) plane.

Download Full Size | PDF

Correct scaling region detection does occur with delay of 1 data point, as in Fig. 10(a), and in this case low values of CD are still obtained consistently. It is uncertain what such a low dimension estimate means. The output power time-series from this region shows large amplitude spikes that resemble a noisy mode locked laser. In theory an ideal mode-locked laser output would be represented by a zero dimensional point in the phase space. Correspondingly, the correlation dimension analysis may show a small, noise induced, nonzero dimension estimate. Thus, we hypothesize that correlation dimension analysis is potentially a useful tool for analyzing mode-locked lasers where dimension 0 is an ideal mode-locked laser and the magnitude of the fractional CD may be an indicator, and possibly a measure of departure from ideality. This hypothesis will need to be further investigated.

It is of interest to also look at other features of the Region III dynamics, such as the spike amplitude and frequency. This is the subject of further research.

Note that the authors recognize that correlation dimension analysis is not designed to be used with highly correlated, periodic data. However, when dealing with systems that switch between periodic and chaotic nonlinear dynamics with small parameter changes it is unavoidable if a single measurand is to be used to describe the nature of the data. It is possible the very low CD estimates seen in some regions of the map could be aberrative due to the dimensional analysis being used on data for which it is not appropriate. Despite this, the technique still gives consistent results and can identify different regions - this is significant, even if the actual value of the CD for region III may not be physically meaningful.

Regions of higher correlation dimension, Region IV in Fig. 6, generally appear for negative detunings at the transition between Regions II and III, and also between Regions III and V. These higher values of CD, in the range 2 to 2.5, indicate a chaotic output of the laser system. The complexity of this chaos is fairly low dimensional compared to other nonlinear laser systems such as a semiconductor laser subject to optical feedback which typically show correlation dimension values of between 3 and 4 [44]. This region could be of interest for chaotic secure communication schemes.

The surrogate data map was constructed by averaging 10 maps generated with surrogate time-series data sets. Inspection of each of the 10 individual surrogate maps reveals that no structure, such as that seen for the original data, is present. For the vast majority of parameter values, the automated program does not detect scaling regions in the D_m(r) graphs when using the surrogate data sets. This is as expected since the surrogate construction process should remove any determinism from the time-series, producing stochastic, infinite dimensional attractors. In each individual surrogate map there are spurious CD values recorded, but never at the same point in the (K, Δω) parameter space so their impact is averaged out in the final map which shows an infinite CD over the entire parameter range. The map of the original experimental data shows values for CD recorded at 19 369 points in the (K, Δω) range (out of 55 440 total points), compared with an average 230 values recorded from the surrogate data sets. Since this surrogate data produces a map that is significantly different to the map generated from the original data we can be confident that the correlation dimension values calculated are a result of deterministic nonlinear dynamics, rather than a linearly filtered stochastic process.

The correlation dimension map generated from the simulated data, shown in Fig. 7, closely resembles the map generated using experimental data. The model correctly predicts Region I where the output is locked and stable. The low CD, Region III, is present for negative detunings and shows the same behavior as the experimental data when the embedding delay is changed. The high CD, chaotic Region IV is also identified in the same parameter space as in the experimental map. As to discrepancies, the main difference is in the low injection strength region where the simulated data map shows no definition between Regions II and V. This arises due to the fact that this whole region is dominated by periodic oscillations – correctly identified in the simulated map as CD = 1. The reason a distinction is seen in the experimental map is due to the signal-to-noise levels. The map of time-series intensity in Fig. 8(b) shows that the two crescent shapes, corresponding to Region II, have larger amplitude than the surrounding region. In this low amplitude region, the unavoidable noise contained in the experimental time-series data results in a signal-to-noise level too low to allow the scaling region to be identified, therefore making it impossible to determine a value for the CD. Since the theoretical model does not include any noise, the small amplitude periodic oscillations of Region V can be successfully analyzed and assigned a value for CD.

On comparison with the bifurcation diagram from [22], reproduced in Fig. 8(a), it can be seen that all the identified regions I-V are separated by bifurcations lines. A period-doubling bifurcation surrounds Region II, whilst a torus bifurcation separates Region III from V. Hopf and saddle-node bifurcation lines occur at the boundaries of Regions I and III for negative detuning, and between Regions I and V for positive detunings. The chaotic dynamics seen in Region IV primarily occur at the higher injection end of the period-doubling bifurcation and along the torus bifurcation, both for negative detuning.

The peak amplitude map from previous investigations [21], also allows identification of different regions, but does not provide any insight to the nature of the dynamics. In contrast, the correlation dimension map (Fig. 6) gives a quantitative measure of the level of complexity of the dynamics at each point in the map. When used in combination with the amplitude map in Fig. 8(a) and the bifurcation diagram in Fig. 8(b), a much clearer picture of the different dynamical regions of the OISSL system emerges.

5. Conclusion

An automated computation process was developed to generate a dynamic map of the correlation dimension of the output power time series from an optically injected Nd:YVO₄ solid state laser with varying injection strength and frequency detuning. The map is in excellent agreement with prior investigations into the amplitude and bifurcation analysis of this laser system [21, 22] and gives further insight in to the level of complexity of the different nonlinear dynamics produced by this system. Maps generated for both experimental and simulated laser outputs show similar structure with discrepancies only arising due to signal-to-noise limitations in the experimental data. The automated process could easily be implemented to analyze other complex systems, from which experimental time-series can be measured, for large ranges of relevant parameter values.