## 3 Results  
The model equations were integrated with the following (typical) parameters that fit the experimental conditions in [19, 20, 21]: \(k=300\) ns\({}^{-1}\), \(\gamma_{N}=1\) ns\({}^{-1}\), \(\alpha=4\), \(\epsilon=0.01\), \(\eta=30\) ns\({}^{-1}\), \(\tau=5\) ns, \(\mu_{dc}=0.99\) and \(D=10^{-5}\) ns\({}^{-2}\). The Euler-Maruyama method with an integration step of \(dt=1\) ps was used to integrate the equations. The initial conditions were chosen with the laser off and after disregarding a transient time of 50 us, the equations were integrated for 10 us. To detect the spike times the intensity time series, \(|E(t)|^{2}\) was band-pass filtered [27] to simulate the finite bandwidth of the experimental detection system, and was normalized to zero mean and unit variance. Then, a spike was detected when the intensity decreased below the threshold \(Th=-1.1\). In Figure 1, we analyze the statistics of the spike sequences, as a function of the modulation amplitude and frequency. Figure 1 displays in color code the number of spikes in the intensity time series (panel a), the mean inter-spike-interval (ISI) normalized to the modulation period (panel b), and the normalized standard deviation of the ISI distribution (panel c, \(\sigma_{IS1}/\langle ISI\rangle\), in logarithmic scale). In panel (a), we see that the spikes become more frequent with increasing \(a_{mod}\) and \(f_{mod}\), i.e., the spikes become faster with \(a_{mod}\) and \(f_{mod}\). In panel (b) we can differentiate five regions: in blue, \(\langle ISI\rangle/T_{mod}\approx 1\), in cyan, \(\langle ISI\rangle/T_{mod}\approx 2\), in green, \(\langle ISI\rangle/T_{mod}\approx 3\), in yellow \(\langle ISI\rangle/T_{mod}\approx 4\), and in red \(\langle ISI\rangle/T_{mod}\geq 5\). We do not find 5:1 locking, and we interpret that this is likely due to the fact that the 5:1 locking region would be located at high modulation frequencies, where the spikes are poorly defined, as the intensity displays quite chaotic fluctuations. To determine where, within these regions, the spikes are regularly locked to the sinusoidal modulation, we need to inspect panel (c), where very low values of \(\sigma_{ISI}/\langle ISI\rangle\) (\(<\)10\({}^{-2}\), blue regions) reveal a narrow ISI distribution (i.e., a regular spike timing). In panel (c), we see two well-defined blue regions (light blue and dark blue). In these regions, in panel (b) we see that \(\langle ISI\rangle/T_{mod}=2\) or 3. Comparing panels (b) and (c), we see that for the blue regions in panel (b) \(\left\{\left\langle ISI\right\rangle/T_{mod}\approx 1\right\}\), \(\sigma_{ISI}/\langle ISI\rangle>10^{-2}\) in panel (c); therefore, the timing of the spikes when \(\langle ISI\rangle\approx T_{mod}\) is significantly more irregular than the timing of 2:1 and 3:1 locked spikes. Examples of the intensity time series in the different locking regions are presented in Figure 2. The shape of the oscillations after each spike is similar for 2:1 and 3:1 locking, but it is very different for 1:1 locking, as in this case the spike is of very small amplitude. For 2:1 and 3:1 locking (panels b and c), in our previous works, we have reported very similar experimental time series (see Figure 1 in [20]). In our previous works, the modulation was of a small amplitude and we did not observe 1:1 locking of the form shown in Figure 2a. However, a very similar experimental time series of 1:1 locked spikes was reported in Ref. [8] (Figure 3f, which is for a modulation amplitude of 13.8% of the DC level of the laser current). Figure 1: **(a)** Number of spikes in the intensity time series, in logarithmic color code, as a function of the modulation amplitude and frequency. **(b)** Mean inter–spike–interval, normalized to the period of the modulation. Here, the color scale is saturated such that red represents \(\left\langle ISI\right\rangle/T_{mod}\geq 5\). **(c)** Standard deviation of the ISI distribution normalized to the mean ISI. Note the logarithmic color scale.