Chaos-preserving reduction of the spin-flip model for VCSELs: failure of the adiabatic elimination of the spin-population difference

TL;DR

This paper investigates the spin-flip model for VCSELs, demonstrating that adiabatic elimination of the spin-population difference fails to preserve key chaotic dynamics, emphasizing the variable's essential role.

Contribution

The study shows that the spin-population difference cannot be adiabatically eliminated without losing critical dynamical features like polarization chaos in the model.

Findings

Adiabatic elimination fails to reproduce polarization chaos.

The spin-population difference is essential for complex dynamics.

Minimal embedding dimension analysis supports model reduction limitations.

Abstract

When studying the dynamics of Vertical-Cavity Surface-Emitting Lasers, and their polarization properties, the spin-flip model appears to be the simplest model qualitatively reproducing all dynamical features that have been observed experimentally. Nonetheless, because of the fast time-scale of the spin-relaxation processes, the specific role and the importance of the spin-population difference - which is one of the specific feature of the spin-flip model - has been continuously questioned. In fact, the debate regarding the possible adiabatic elimination of the spin-population remains fully open. In this paper, our goal is to bring new light into this issue by demonstrating that this variable is essential to preserve the most complex dynamical features predicted by the spin-flip model, such as polarization chaos, and, therefore, needs to be conserved. To do so, we first perform a detailed analysis, focusing on the chaotic dynamics, to determine the minimal embedding dimension for the spin-flip model. As the latter confirms that a reduction of the model could be envisaged, we then consider the adiabatic elimination of the spin-population difference to highlight and explain its failure to reproduce essential dynamical features obtained in the original model.

Figures

Equations

$$\frac{d{R}_{+}}{dt}=$$

$$\frac{d{R}_{-}}{dt}=$$

$$\frac{d\Phi }{dt}=$$

$$+\left(\frac{R_{+}}{R_{-}}+\frac{R_{+}}{R_{-}}\right){\gamma }_{a}sin(\Phi )$$

$$\frac{dN}{dt}=$$

$$\frac{dn}{dt}=$$

$$\mathbf{P}=\mathbf{XW}$$

$${\mathbf{X}}_{rec}=\mathbf{P}(:,1:r)\mathbf{W}(:,1:r{)}^T={\mathbf{P}}_{trunc}{\mathbf{W}}_{trunc}^T$$

$$\mathbf{P}(:,k)=\mathbf{XW}(:,k)$$

$$n=\frac{N(R_{-}^2-R_{+}^2)}{{\gamma }_s/\gamma +R_{+}^2+R_{-}^2}$$

$${\lambda }_1=\left(-2{\gamma }_a-\chi \right)\pm \sqrt{{\chi }^2-4\alpha {\gamma }_p\chi -4{\gamma }_p^2}$$

$$\chi =\frac{\gamma \kappa }{{\gamma }_s}\left(\mu -1+\frac{{\gamma }_a}{\kappa }\right)$$

$$\mu >1-\frac{{\gamma }_a}{\kappa }-2\frac{{\gamma }_a{\gamma }_s}{\kappa \gamma }$$

$$\sqrt{{\chi }^2-4\alpha {\gamma }_p\chi -4{\gamma }_p^2}<2{\gamma }_a+\chi$$

$$0<({\gamma }_a^2+{\gamma }_p^2)+({\gamma }_a+\alpha {\gamma }_p)\chi$$

$$\chi$$

$$\text{i.e.}\mu$$

$$N=1-\frac{{\gamma }_a}{\kappa },\varphi =\pi ,R=\sqrt{\frac{\mu -N}{2N}}$$

$$\left[\begin{array}{c}\dot{\delta R_{+}}\\ \dot{\delta R_{-}}\\ \dot{\delta \varphi }\\ \dot{\delta N}\end{array}\right]=\left[\begin{array}{cccc}\kappa (N-1)-2\kappa \gamma N{R}^2/{\gamma }_s& {\gamma }_a+2\kappa \gamma N{R}^2/{\gamma }_s& {\gamma }_{p}R& \kappa \gamma R\\ {\gamma }_a+2\kappa \gamma N{R}^2/{\gamma }_s& \kappa (N-1)-2\kappa \gamma N{R}^2/{\gamma }_s& -{\gamma }_{p}R& \kappa \gamma R\\ -4\kappa \gamma \alpha NR/{\gamma }_s-2{\gamma }_p/R& 4\kappa \gamma \alpha NR/{\gamma }_s+2{\gamma }_p/R& -2{\gamma }_a& 0\\ -2\gamma NR& -2\gamma NR& 0& -\gamma (1+2R^2)\end{array}\right]\left[\begin{array}{c}\delta R_{+}\\ \delta R_{-}\\ \delta \varphi \\ \delta N\end{array}\right]$$

$$\chi =\frac{\gamma \kappa }{{\gamma }_s}\left(\mu -1+\frac{{\gamma }_a}{\kappa }\right)$$

$${\lambda }_1=\left(-2{\gamma }_a-\chi \right)\pm \sqrt{{\chi }^2-4\alpha {\gamma }_p\chi -4{\gamma }_p^2}$$

$${\lambda }_2=\frac{\gamma \kappa \mu }{2({\gamma }_a-\kappa )}\pm \sqrt{\frac{{\gamma }^2{\kappa }^2{\mu }^2}{4({\gamma }_a-\kappa {)}^2}-2\gamma \kappa \mu +2\gamma \kappa (\frac{{\gamma }_a}{\kappa }-1)}$$

Full text

In this paper, our goal is to bring new light into this issue by demonstrating that this variable is essential to preserve the most complex dynamical features predicted by the spin-flip model, such as polarization chaos, and, therefore, needs to be conserved. To do so, we first perform a detailed analysis, focusing on the chaotic dynamics, to determine the minimal embedding dimension for the spin-flip model. As the latter confirms that a reduction of the model could be envisaged, we then consider the adiabatic elimination of the spin-population difference to highlight and explain its failure to reproduce essential dynamical features obtained in the original model.

I Introduction

Due to their intrinsic advantages over standard edge-emitting semiconductor lasers, Vertical-Cavity Surface-Emitting Lasers (VCSELs) are surely making their way towards a widespread use not only in telecom but also for illumination - as in the iPhone X - or gesture recognition and other smart sensing applications. Their main drawback has been identified early on as polarization instabilities [1]: the most striking example being polarization switching (PS) events where the linearly polarized laser light suddenly rotates by $90^{\circ}$ [2, 3, 4, 5, 6]. But a wide range of behaviour have been reported including elliptical polarization [7], noise-induced hopping [8, 9] and even chaotic dynamics with the so-called polarization chaos [10, 11]. Although these instabilities have been first studied to be avoided [12, 13], it was also remarked that polarization dynamics could be embraced instead: thus leading to the tuning of polarization oscillations for high-speed spin-VCSELs [14, 15], polarization switching as optical memory mechanism [16, 17, 18] or random bit generation using polarization chaos [19].

The polarization features of VCSELs are accurately described using the so-called spin-flip model (SFM) [20, 21]. The main features of the SFM are that it takes into account: 1) the phase and amplitude anisotropies of the gain medium; 2) two distinct pools of carrier depending on their spin, each of them connected to circularly polarized light of a certain handedness. The spin relaxation processes connecting the two carrier pools are then modeled by a simple coefficient commonly called the spin-flip rate. The SFM is, so far, the simplest model qualitatively predicting all the dynamical polarization features experimentally observed in VCSELs. Indeed, other models have been proposed for a given behaviour [22, 8], but none - apart from the SFM - reproduced the wide variety of features observed experimentally. On the other hand, the nonlinearities of the SFM hinder the physical interpretation. To get a better understanding and gain further physical insight, different model reductions - based on various approximations - have been proposed [23, 24, 25]. Experimentally, a rather large value for the spin-flip rate has typically been reported [26, 27, 28]: the spin relaxation rate typically seems to be much faster than the other time-scales of the system, in particular the carrier and field decay rates. As a result, adiabatically removing the spin-population difference, i.e. only considering one carrier pool, was often considered to be a reasonable and realistic approximation, especially when focusing on the system steady-states[23, 24, 25]. On the contrary, a rather low value of the spin-flip rate seems to be a crucial ingredient to obtain polarization chaos in solitary VCSELs [29]. Obviously, this observation suggest that, at least for these peculiar experimental cases [10, 11], the spin-population difference could actually play a significant role in the laser dynamics. Hence, we can naturally wonder where is the limit of validity for the approximation allowing the adiabatic elimination of the spin-population difference.

In this work, we first aim at confirming that the SFM dimension, i.e. the number of independent variables, could indeed be reduced below the current state of the art. We determine the minimal embedding dimension of the system, i.e. the minimal number of variables required to accurately represent its behaviour, using numerical processing of simulated data based on Principal Component Analysis (PCA) [30] along with the ”False Neighbors” technique [31, 32]. On top of this important insight, which indicates that a reduction could indeed be considered, we also identify from our analysis data the spin-population difference as a potential candidate for such reduction. Yet, we then show that an adiabatic elimination of the corresponding variable performs, in fact, very poorly at preserving system behavior such as polarization chaos, but also the stability of steady-state solutions. We therefore highlight that the spin-flip model could indeed be further reduced without qualitative loss, but the model-order reduction approach exploiting the adibatic elimination of the spin-population difference is not a suitable solution to accomplish this goal.

II Spin-flip model and parameters

In its original version, the SFM was a 6-dimension model which has first been introduced in [20]. When investigating its dynamical properties [21], it was quickly observed, however, that only the phase difference between the two linear polarization modes played a significant dynamical role. Using phase-amplitude decomposition, it was then proposed to use this phase-difference to get a 5-dimensional model [25]. This version of the SFM is the one we exploit here and reads as follows:

[TABLE]

with $R_{\pm}$ the amplitude of the right and left circular polarizations, $\Phi$ the phase difference between them, $N$ the total carrier population and $n$ the carrier population difference between the two carrier reservoirs for each circular polarization [20]. Because the carriers in these two reservoirs have different spins, the carrier population difference is often referred to as the spin-population difference. The parameters are as follows: $\kappa$ and $\gamma$ are the field and carrier decay rates respectively, $\gamma_{s}$ is the spin flip relaxation rate, $\alpha$ is the linewidth enhancement factor, $\gamma_{p}$ and $\gamma_{a}$ are the phase and amplitude anisotropies respectively. For simplicity, no misalignment between amplitude and phase anisotropies is considered [3, 33, 34].

Unless stated otherwise, we use the same parameter values as in previous work [10, 35]: $\kappa=600\,ns^{-1}$, $\gamma_{s}=100\,ns^{-1}$, $\gamma_{a}=-0.7\,ns^{-1}$, $\gamma=1\,ns^{-1}$ and $\alpha=3$. We also consider two distinct cases of chaotic dynamics as discussed in [35] which are obtained for $\gamma_{p}=4\,ns^{-1}$ (case 1) and $\gamma_{p}=25\,ns^{-1}$ (case 2). The injection current is normalized by the current at threshold, i.e. $\mu=1$ corresponds to the laser threshold when no anisotropies are considered [21]. The normalized injection current $\mu$ is varied from 1 to 10; the higher limit corresponds, in this theoretical framework, to 10 times the laser threshold without anisotropies which is already largely above the current range commonly considered [11]. Although more complex dependencies appear to be needed for an accurate estimation of the laser threshold, see e.g. eq. 11 in [27], we believe that this higher limit provides a reasonable first order approximation for comparison with experimental observations. Finally, the linearly polarized (LP) steady-state stable at threshold for the parameters considered here will be identified as X-LP (which is also the low-frequency eigenmode) while the orthogonal LP steady-state will be identified as Y-LP (the high-frequency eigenmode). This notation is identical to the one used in previous work [10, 35].

III Optimal embedding of polarization chaos

The first part of this work is dedicated to the determination of the minimal number of variables required to accurately describe the polarization chaos dynamics. As such, this will also tell us whether our model could potentially be further reduced and simplified without any loss of generality. A common approach to do so is to use the so-called Principal Component Analysis (PCA) technique[30]. This technique performs a linear mapping of the data to a lower-dimensional space in such a way that the variance of the data in the low-dimensional representation is maximized. It is a popular technique for dimensionality reduction. The original variables are transformed into new variables, called principal components (PCs), by an orthonormal linear transformation and are ordered by decreasing variance. Using the 5-dimension simulated data generated by the SFM model, we generate a set of 5 linearly independent PCs describing the system trajectory in the variable space. Here, it is however important to note that different pre-processing steps might applied to the data before applying PCA. As any pre-processing step on the data can influence the outcome of the PCA, the user needs to carefully choose them. In this paper, we will use the common centering pre-processing step, i.e. the removal of the mean value of the different variables, and will compare results of the PCA with and without this centering.

The PCA technique provides a decomposition of the original SFM data matrix $\mathbf{X}$ of size $\#time\ samples$ x $\#states$ (i.e. $\#time\ samples$ x $\#variables$) as

[TABLE]

where the matrix $\mathbf{W}$ has dimension $\#states$ x $\#states$. The symbol $\#$ is equivalent to ”number of”. The number of columns of $\mathbf{P}$ and $\mathbf{W}$ can be truncated and a reconstructed $\mathbf{X}_{rec}$ can be computed as

[TABLE]

where $\mathbf{P}(:,1:r),\mathbf{W}(:,1:r)$ denote the first $r$ columns of the corresponding matrices. It is important to highlight that each principal component $\mathbf{P}(:,k),k=1,...,5$ can be expressed as a linear combination as

[TABLE]

where $\mathbf{P}(:,k),\mathbf{W}(:,k)$ denote the $k-th$ column of the corresponding matrices. Here, we obtain a set of five linear combinations of variables that are linearly independent of each other in order of decreasing variance. Considering both case 1 ($\gamma_{p}=4\,ns^{-1}$) and case 2 ($\gamma_{p}=25\,ns^{-1}$) for increasing injection currents, we perform PCA independently for each current step. We thus obtain Fig. 1 that shows the variance associated with each PC after a PCA decomposition with and without centering of the initial data. The regions for which polarization chaos is observed are highlighted by the light gray background. The results are very consistent between case 1 and case 2, but we also see an influence of the centering. In all panels, we clearly see that at least one PC exhibits a variance significantly smaller than the others, more than 50 dB weaker than the variance of the dominant PCs. This therefore suggests that there would indeed be some potential to reduce the number of variables of the model without loss of generality. When the centering pre-processing is used, two PCs - instead of only one - show smaller variance than the dominant PCs. Interpreting this result as a possible reduction of two orders, i.e. possible reduction to a 3D model, might be an overstatement, but these results strongly suggest that at least one variable could be removed.

When dealing with chaotic dynamics, it is obvious that any small changes of the system should be considered with the highest level of caution as these could have a dramatic impact on the system behaviour. To verify the suitability of an embedding for a chaotic systems, a common test is the so-called ”False Neighbors” test [31, 32]. If the embedding dimension is too small, the chaotic attractor will necessarily be folded, i.e. some parts of the chaotic trajectory will artificially become close and thus create the so-called ”False Neighbors”. To detect such folding, we use the following approach. After applying PCA, we reconstruct the 5-D states trajectories using only a subset of the calculated PCs using (7). We then compare the distance between the data point and its closest neighbor in the initial and reconstructed 5-dimension state trajectories space. From these data, we derive two figures of merit:

• •

the Mean Distance Growth (MDG): the average growth for all selected closest neighbors for all data points considered. Obviously, a large MDG $>>1$ indicates a probable folding while a MDG very close to 1 will tend to confirm that the embedding dimension is large enough to contain the chaotic trajectory.

• •

the amount of “False Neighbors” (FN): the number of closest neighbors - for all data points - for which the distance growth exceed a given threshold $R_{tol}$. The amount of FN of course depends on the threshold value, but, as discussed in [31], the approach is quite robust against the choice of threshold. Here, we show results for $R_{tol}=20$, but verified that the same conclusions hold for threshold values from 2 to 30.

Here, we consider 10000 data point per time-series with a time-step of 25 ps. We take the single closest neighbor for each data point [31] and impose a minimal time separation of 1000 ps between them to avoid correlation in time. The evolution of these two figures of merit for increasing currents and for different embedding dimensions are shown in Figs. 2-3 without and with centering, respectively. Considering both data sets, it is clear that considering only the first 2 PCs is largely insufficient as the amount of false neighbors is quite large: for both cases, inside the chaotic region, the closest neighbor after embedding is most of the time a ”false” neighbor. Adding the third PC induces a significant improvements, and only a few false neighbors are detected - between 1 to 4 % without centering (Fig. 2) and less than 0.5 % with centering (Fig. 3) - and only inside the chaotic region. But it is only when considering the first 4 PCs that the amount of detected False Neighbors effectively reach and remain at 0 for all conditions. From a system analysis viewpoint, the imperfect embedding using only 3 PCs could be considered for a first approximation even though it would miss some details of the dynamics. On the other hand, from a nonlinear dynamics viewpoint, a 4 dimension embedding using the first 4 PCs seems to be sufficient to contain the whole chaotic attractor without inducing any detectable folding. Indeed, in this case absolutely no false neighbors is detected and low value of the mean distance growth is systematically observed. In essence, these results confirm that the 5th PC, i.e. the PC with the smallest variance, does not provide essential information over the system behaviour or chaotic features. Therefore, we can conclude that these results confirm that the spin-flip model could be reduced from 5 to 4 variables while preserving all complex features including chaos.

Considering this perspective, we can further analyze the different outputs of the PCA to identify the variables that could effectively be dismissed. Such information can be provided by the $\mathbf{W}$ matrix which describes the linear combination of the original variables leading to the PC description. To put it differently, the $\mathbf{W}$ matrix shows the contribution of each variable to the different PCs. Most variables are actually contributing to different PCs in a rather complex way. In Fig. 4, we show the contribution of the spin-population difference $n$ to each of the PCs. Thus, it is rather striking that the spin-population difference is, in fact, almost only contributing to the 5th PC, with the exception of Fig. 4(b.1) that we discuss below. It is however important to understand that this does not mean that the 5th PC is almost perfectly equal to the $n$ variable, since the other SFM variables can also have an impact on the 5th PC in the linear combination formula (8). The result shown in Fig. 4 suggests that the spin-population difference $n$ could be a good candidate to be considered for the model reduction as we will do in the next section using an adiabatic elimination. The failure of this adiabatic elimination approach, which we will discuss in what follows, should however not lead to the conclusion that the SFM cannot be reduced to 4 variables despite the outcome of the PCA and ”False Neighbors” analyses, but only to the fact that another approach needs to be considered to obtain a suitable model-order reduction. Finally, in Fig. 4(b.1), i.e. for PCA with centering of case 2: in this case, the spin-population difference also appears to contribute to the 4th PC shown in purple. Yet, this result needs to be put in perspective with the variance of this 4th PC as shown in Fig. 1(b.1): there, we see that the 4th and 5th PC both have very low variance, more than 50 dB below the three other PCs. Hence it does not invalidate the proposed interpretation as the spin-population difference $n$ never contributes to one of the dominant PC. Nevertheless, this might eventually be seen as a hint of the failure of the adiabatic elimination reduction, but this interpretation would be far too speculative at this stage.

IV Adiabatic elimination of the carrier population difference

As discussed above, the dynamical analysis that we performed suggests that the spin-population difference $n$ could be a good candidate for the reduction of the spin-flip model. Since a similar conclusion can be reached from physical considerations [23, 24], this approach obviously need to be further investigated.

Based on the rate equations describing the spin-flip model, and similarly to what has been done in previous works [23, 24], we can eliminate the carrier population difference adiabatically. Starting with (5) and considering that the spin-flip rate $\gamma_{s}$ is large enough for $n$ to reach a steady-state much faster than the other variables, we directly obtain that:

[TABLE]

The differences with the expression found in [24] are only due to the particular approximation made therein and that we do not apply here. We can then easily confirm that this expression of the carrier population difference is an accurate approximation by comparing it to the simulated values. By doing so, we obtain that the error is typically two-orders of magnitude smaller than the value taken by the carrier population difference variable. As such, although such small error will have an obvious impact on the simulated chaotic time-series, we would expect at this point a negligible qualitative impact on the laser behaviour. On the contrary, when computing bifurcation diagrams - as shown in Fig. 5 - using this reduced model, we can only observe rather dramatic qualitative changes. In particular, in case 2, the initially wide chaotic region is shifted and shrunk into a narrow range of injection current. Thus, for currents above 2.7, the laser only exhibits a stationary behaviour corresponding to Y-LP, i.e. the linearly polarized steady-state orthogonal to the steady-state stable at threshold X-LP. This is particularly puzzling considering that, in case 2, the Y-LP steady-state is supposed to be unstable up to extremely high values of current ($\mu>50$) [35].

To clarify this issue, we have re-investigated the stability of both the X-LP and Y-LP steady-state for the full SFM and for the reduced model. To compute it, we linearize the rate equations around the steady-states and numerically extract the eigenvalues of the linear system. By adiabatically eliminating the carrier population difference, we obviously remove one equation. But using the expression of $n$ given in (9) also creates several additional dependencies that are nonexistent in the full model. While no impact of the adiabatic elimination is observed on the stability of the X-LP state (not shown), we observe that the stability of the Y-LP state is not well-preserved. As displayed in Fig. 6, the lower boundary of the stability region with respect to the injection current appears to be independent of the birefringence for the reduced model while a clear dependence can be observed in the full-model. As could be expected, this discrepancy is reduced as the value of the spin-flip rate is increased, but this agreement is still imperfect and restricted to low values of the birefringence $\gamma_{p}$.

We now investigate analytically the stability of the Y-LP steady-state. As described in the appendix, the linearized system has two pairs of eigenvalues $\lambda_{1}$ and $\lambda_{2}$ of the form $A\pm\sqrt{B}$. Because $\lambda_{2}$ always corresponds to a pair of complex conjugated eigenvalues with a negative real part for the typical range of parameter values that we consider here, we will focus on $\lambda_{1}$ which effectively limits the stability of Y-LP. The latter can be expressed as follows:

[TABLE]

with the variable $\chi$ being defined as:

[TABLE]

In practice, $\chi$ corresponds to the injection current, normalized by the Y-LP steady-state threshold $\mu_{Y-LP,th}=1-\gamma_{a}/\kappa$ - meaning that $\chi$ is necessarily positive - and with a peculiar scaling. Inside the square root term in (10), we have a convex $2^{nd}$ order polynomial equation in $\chi$. The roots of the equation are $2\gamma_{p}(\alpha-\sqrt{\alpha^{2}+1})$ and $2\gamma_{p}(\alpha+\sqrt{\alpha^{2}+1})$, which are respectively negative and positive for $\gamma_{p}>0$. In this case, the second term represents a crucial threshold with respect to $\chi$: above it, we have real eigenvalues, while, below it, we have a pair of complex conjugated eigenvalues. Thus, it is important to remark that, with the set of parameters considered in this work, this threshold is increasing along with $\gamma_{p}$ and is already above $\mu=5$ for $\gamma_{p}=2\,ns^{-1}$. As a result, for most cases, the latter will apply, and the stability is then defined by the first term of ( 10), i.e. the real part of the eigenvalue pair. From its expression, we immediately get that the real part is negative if:

[TABLE]

This equation is almost the same as eq. 55 of [25], except for the 2nd term $\gamma_{a}/\kappa$ which is however expected to have a negligible impact as we typically have $\gamma_{a}<<\kappa$. This is clearly the stability limit for the Y-LP state visible in Fig. 6 which correspond to a well-known Hopf bifurcation [21, 25].

When we have two real eigenvalues, since the square root term is necessarily positive, the largest eigenvalue will necessarily be the one for which the two terms are added. In this case, it can only be negative if the following condition is fulfilled:

[TABLE]

To have a real eigenvalue, the left hand side must be positive. This implies that $2\gamma_{a}+\chi>0$ - which corresponds to (12) - is a necessary condition for the previous inequality to be verified. If both sides are positive, the previous equation is equivalent to:

[TABLE]

If $(\gamma_{a}+\alpha\gamma_{p})>0$, this inequality is obviously always verified. For $\gamma_{a}<0$ and $\gamma_{p}>0$, this comes down to a threshold on the birefringence: $\alpha\gamma_{p}>|\gamma_{a}|$, which, with our typical set of parameters, corresponds to values of the birefringence above approximately $0.23$. On the other hand, when $\gamma_{a}+\alpha\gamma_{p}<0$ the Y-LP steady-state appears to be only stable when:

[TABLE]

Apart from the extra $\gamma_{a}/\kappa$ term which is expected to be negligible, this expression is similar to eq. 53 found in [25] but for the X-LP steady state and with a change of sign. This difference is easily explained by the fact that X-LP and Y-LP are interchangeable by simultaneously changing the sign of $\gamma_{a}$ and $\gamma_{p}$. This leads to the same change of sign, as shown in fig. 1 (a) and (b) of [25], where the same curve is obtained for $\mu_{xs}$ and $\mu_{ys}$. Using the same approach as for Fig. 6, we could also confirm that an excellent agreement is obtained between the reduced and complete model for $\gamma_{p}<\gamma_{a}/\alpha$.

To conclude, for the case of $\gamma_{a}<0$ and $\gamma_{p}>0$, the stability of the Y-LP steady-state is largely impacted by the adiabatic elimination of the spin-population inversion. While for sufficiently small values of the birefringence $\gamma_{p}$, an acceptable agreement between the reduced and complete model can be obtained, significant discrepancies arise when $\gamma_{p}$ is increased. In the reduced model as soon as $\gamma_{p}$ is larger than $\gamma_{a}/\alpha$, the Y-LP state is always stable when the condition of eq. 12 is met. This means that, intrinsically, the stability of Y-LP becomes independent of $\gamma_{p}$. Although this reduction could be a sufficient but limited approximation for very large spin-flip rates $\gamma_{s}>1000\,ns^{-1}$ and small birefringence values $\gamma_{p}\approx 5\,ns^{-1}$, we showed in our previous work [29] that these are, by far, conditions not suitable to generate polarization chaos dynamics.

V Discussion

Looking at all in-depth investigations of VCSELs dynamics, it is undeniable that the spin-flip rate is an essential parameter controlling the dynamical behaviour of VCSELs. Yet, the spin-flip rate only plays a role of the laser dynamics through its influence on the spin population difference. Often, it is assumed - also considering the typically large values of the spin-flip rate reported experimentally - that the adiabatic elimination of the spin-population difference could be a suitable approximation and would retain all essential non-linearities in the SFM model. Yet, in this work, we highlight that the impact of this adiabatic elimination on the polarization chaos and stability of the system steady-states is far too large to be neglected, even for large values of the spin-flip rate and, especially, for birefringence values in the order of tens of $ns^{-1}$ and above [14, 15]. However, our analysis using the PCA and ”False Neighbors” techniques confirm that there is potential for model order reduction, but it needs to be performed using methods that, unlike the adiabatic elimination of the spin-population, would preserve all key dynamical features of the model. However, is is an aspect that we will leave for future work. Model order reduction is a vast field [36, 37, 38, 39]. Depending on the complexity of the system representation (e.g., linear time-invariant, time-variant, non linear, etc.) to be reduced, the model order reduction approach has an increasing complexity especially if system properties need to be preserved in the reduced model. Model order reduction schemes very often use specific mathematical transformation to transform the original state vector space into a reduced space, while keeping accuracy and properties with respect to the original model. The property-preserving model order reduction is much more complex than just an accuracy-preserving model order reduction. Already in the case of linear time-invariant systems, model order reduction techniques can provide accurate reduced order models of the transfer function behavior in a frequency bandwidth of interest, but they can fail in preserving properties such as stability and passivity. For nonlinear systems, it is obviously even more complex. Therefore, a dedicated investigation of such techniques that can perform model order reduction of the original SFM equations are beyond the scope of this contribution and deserves a dedicated effort in a separated work.

Finally, we would like to discuss why the work presented here leads to conclusions that significantly differ from previous reports, in particular [23, 24]. We believe that this can be mostly attributed to the different focus of our analysis: while the authors of [23, 24] were working to get a better physical understanding of the stability boundaries of the linearly polarized steady-states, we look here at the modelling accuracy of the dynamical features. Similarly, while strong approximations are taken in [23, 24], we tried to keep a very general viewpoint, thus did not use similar approximation. In the end, we are convinced that these two aspects are of course complementary. However, as highlighted in this work, when considering dynamical behaviours, the spin-population difference appears as an essential piece of the puzzle that must be taken into account to its full extent to ensure qualitative relevance of the modelling.

Appendix A Stability of the Y-LP steady-state in the adiabatically reduced spin-flip model

The Y-LP steady-state, is defined as follows:

[TABLE]

For simplicity, we will approximate the expression of the adiabatically eliminated population difference by considering that the denominator can be directly expressed as $\gamma_{s}/\gamma$. Since this term is typically of the order of 100 to 1000 while the amplitude terms are expected to be of the order of 1, the approximation is small for the practical cases considered here. After linearization of the equations, we obtain the following system whose eigenvalues will determine the stability of the steady-state.

{strip}

[TABLE]

Using a symbolic computation software, we obtain two pairs of eigenvalues of the form $\lambda_{1}=A1\pm\sqrt{B1}$ and $\lambda_{2}=A2\pm\sqrt{B2}$. Using the following change of variable

[TABLE]

the first eigenvalue pair can be efficiently expressed as follows:

[TABLE]

The second pair of eigenvalue can be expressed as:

[TABLE]

As mentioned in the core of the article, $\lambda_{1}$ is the interesting pair of eigenvalue and is therefore further described therein. The impact of $\lambda_{2}$ is far more limited as will be described here.

Since we typically have $|\gamma_{a}|<<\kappa$ with $\kappa$ being positive by design, the first term will always be negative for the considered parameter range. If the square root is imaginary - which is typically the case for the standard parameter range considered here -, we obtain two complex conjugate eigenvalues with negative real part. On the other hand, if it is real, its absolute value is smaller than the absolute value of the first term, and we therefore have two negative real eigenvalues. Either way, the system will not be destabilize by this second pair of eigenvalue.