Time crystal dynamics in a weakly modulated stochastic time delayed system

Time crystal oscillations in interacting, periodically driven many-particle systems are highly regular oscillations that persist for long periods of time, are robust to perturbations, and whose frequency differs from the frequency of the driving signal. Making use of underlying similarities of spatially-extended systems and time-delayed systems (TDSs), we present an experimental demonstration of time-crystal-like behavior in a stochastic, weakly modulated TDS. We consider a semiconductor laser near threshold with delayed feedback, whose output intensity shows abrupt spikes at irregular times. When the laser current is driven with a small-amplitude periodic signal we show that the interaction of delayed feedback and modulation can generate long-range regularity in the timing of the spikes, which lock to the modulation and, despite the presence of noise, remain in phase over thousands of modulation cycles. With pulsed modulation we find harmonic and subharmonic locking, while with sinusoidal modulation, we find only subharmonic locking, which is a characteristic feature of time-crystal behavior.

In many-particle systems, when the translation symmetry in space (in time) is spontaneously broken, the result is a space (time) crystal. Time-crystal states are characterized by highly regular oscillations that are stable over very long times, are robust under perturbations ("rigidity") and break time-translation symmetry 1,2 . Time-crystal states were originally introduced for many-particle systems in equilibrium 3 , and after the possibility of such states in equilibrium systems was ruled out 4,5 , they were investigated in non-equilibrium systems under a periodic drive 6 . In this case, the periodic forcing defines the discrete time translation symmetry and this symmetry is broken in the time crystal phase, where ' discrete time-crystals' states occur [7][8][9] . In these states the system's variable displays sub-harmonic oscillations that are rigidly locked to the driving signal. In recent years, time-crystal behavior has been studied in a wide range of systems [10][11][12][13][14][15][16][17][18][19] .
A key requirement for observing discrete time-crystal dynamics is that the system's oscillations have longterm regularity that results from the collective synchronization of many interacting degrees of freedom. This excludes period-doubling oscillations in low-dimensional dynamical systems, and it also excludes oscillations in mode-locked lasers, which arise due to interactions of many modes but which lack long-term regularity because, due to noise, they do not remain in phase for long times 1 .
Here we address the following questions: Can a feedback loop counteract the effect of noise? Can a delayed feedback loop generate long-term order? Can we find discrete time-crystal-like behavior in periodically driven stochastic systems with feedback loops?
Time-delayed systems (TDSs) with feedback loops, governed by equations of the form du(t)/dt = f (u(t), t) + Ku(t − τ ) , have an infinite phase space because the initial condition is the function u(t) defined in [-τ,0] 20,21 . In this type of TDS, when the feedback delay time, τ , is longer than the internal characteristic time-scale of system, analogies have been found with the dynamics of one-dimensional spatially extended systems (1D SESs) governed by equations of the form ∂u(x, t)/∂t = f (u, x, t) + D∂ 2 u/∂x 2 with x(t) in [0, L]. Specifically, in the TDS, τ plays the role of the size, L, in a 1D SES. When τ is long enough, using the so-called space-time representation 22,23 , complex spatio-temporal phenomena has been found, such as pattern formation and propagation of defects and localized structures, analogous to that occurring in 1D SESs [24][25][26][27][28][29][30][31][32][33] . These analogies make periodically driven stochastic TDSs promising test benches for finding non-trivial time-crystal-like oscillations.
A semiconductor laser with optical feedback is a well-known TDS, where the delay is due to the finite propagation time of the feedback light. Optical feedback from a distant reflector introduces a set of new modes (the www.nature.com/scientificreports/ so-called external cavity modes) and generates a complex intensity dynamics 34 that, when viewed using the space-time representation, reveals spatio-temporal patterns and defects typical of 1D SES 35 . Near threshold and for appropriated feedback parameters the laser output intensity displays spikes that occur at irregular times. A typical intensity time series is shown in Fig. 1a and a video of the spiking dynamics can be found in 36 . In this regime the laser dynamics is strongly influenced by noise: simulations of the well-known Lang-Kobayashi model 37 indicate that when a noise term is not included in the model, the spikes are periodic or are transient [38][39][40] . When the laser current is modulated with a small-amplitude periodic signal, the spikes tend to occur at intervals of time that are integer multiples of the period of the modulation [41][42][43][44][45][46][47][48][49] .
We have recently shown 50 that a small-amplitude pulsed modulation (less than 2.5% of the dc value of the laser current) generates regular 1:1 locked spikes. A typical intensity time series is shown in Fig. 1b (details of the experimental setup and parameters can be found in 50 ). In Fig. 1b we see that the spikes are periodic but in between the spikes the intensity shows irregular fluctuations.

Results
To determine whether the laser with delayed feedback and weak modulation can display non-trivial discrete time-crystal-like dynamics, we begin by analyzing how regular the spike timing is, and how the spike timing regularity is affected by the waveform and by the parameters of the driving signal. The experimental dataset (described in 50 ) consists of time series of the laser output intensity recorded over an interval of 5 ms. We analyze the effects of sinusoidal and pulsed modulation, varying the modulation frequency, f mod , and the dc value of the laser current, I dc , keeping constant the modulation amplitude (0.631 mA). The intensity time traces contain between 9000 spikes (for low I dc and f mod ) and 120000 spikes (for high I dc and f).
In Fig. 2 we compare the distribution of time intervals between consecutive spikes (inter-spike-intervals, ISIs) generated by pulsed modulation (left) and by sinusoidal modulation (right), as a function of the modulation frequency. The vertical axis displays the ISI normalized to the modulation period, T mod , and the distribution of ISI values in shown in log color code. We see that for particular modulation frequencies the ISI distribution is narrow, and the ISIs are multiples of T mod . For the pulsed waveform, Fig. 2a, there are three "plateaus" where ISI/T mod =1 or 2 or 3 (a spike is emitted every one, two or three modulation cycles), but for the sinusoidal waveform, Fig. 2b, there are only two plateaus where ISI/T mod =2 or 3. www.nature.com/scientificreports/ Therefore, pulsed modulation generates harmonic and subharmonic locking, while sinusoidal modulation, only subharmonic locking. In Fig. 2b we note that for f mod ≈ 3 MHz, �ISI� /T mod ≈ 1 , i.e., the sinusoidal waveform generates, on average, a spike per modulation cycle; however, there is no locking because the ISI distribution is broad and there is no plateau, i.e., there is no interval of frequencies in which �ISI� /T mod ≈ 1 . These observations are consistent with previous experiments using small-amplitude sinusoidal modulation 47,48 , where we found subharmonic locking but not 1:1 locking.
The lack of 1:1 locking generated by small-amplitude sinusoidal modulation contrasts with the dynamics of low-dimensional (non-delayed) dynamical systems, where small-amplitude sinusoidal modulation can generate harmonic locking 51,52 . It also differs from the stochastic resonance phenomenon 53 that has been observed in many noisy oscillators, including in a similar laser system 45 , whose characteristic feature is a peak in the ISI distribution at ISI/T mod =1, whose strength is maximum for appropriated parameters.
With sinusoidal current modulation harmonic locking has been observed experimentally; however, under large-amplitude modulation ( ∼14% and 20% of I dc in 42 and 49 respectively). With large-amplitude modulation the feedback-induced spikes are a perturbation of a sinusoidal-like oscillation. To shed light on this situation we have carried out simulations of the Lang-Kobayashi (LK) model 37 and found a good agreement with the observations: large-amplitude sinusoidal current modulation produced harmonic locking, while small-amplitude modulation produced sub-harmonically locked spikes 54 .
To precisely quantify the regularity of the spike timing we calculate the Fano factor 55 , F, which is a well-known measure of the variability of sequences of events 56 , described in Sec. Methods. Figure 3 displays the Fano factor (in log color code) of sequences of spikes recorded for different I dc and f mod . In panel (a) the modulation is pulsed, while in (b), it is sinusoidal. In both panels we see parameter regions where F < 10 −2 (dark-blue) and regions where F > 1 (yellow). We remark that F < 1 ( F > 1 ) indicates that, in the time scale of 5 µ s, the sequence of counts is more (less) regular than a sequence of events generated by a HPP process.    Fig. 2 for I dc = 26 mA, we see (as expected) that the blue regions in Fig. 3 correspond to the "plateaus" in Fig. 2 where the ISI distribution is narrow and thus, the spike timing is regular. In Fig. 3 we see that for the pulsed signal there are three blue regions, while for the sinusoidal signal, only two, confirming that 1:1 locking is not generated by small-amplitude sinusoidal modulation.
We interpret the different response of the laser to low-frequency modulation, pulsed or sinusoidal, as due to two reasons. First, the pulsed signal produces small but abrupt variations of the pump current that precisely define the timing of the spikes (the current changes at instants of time, producing spikes at those times); in contrast, the sinusoidal signal varies the laser current gradually, and this variation is unable to precisely define the spike timing. Second, the pulsed signal decreases the current for only short time intervals, while the sinusoidal signal decreases the current during longer intervals, which may lead to a stronger influence of noise.
To analyze the time scale of the spike timing regularity, we study how F depends on the duration of the counting interval, T int , and on the modulation frequency, f mod (in color code). Figure 4 displays results for pulsed modulation, when I dc = 26 mA. Panel (a) displays F vs. T int , and we see that, as expected, F → 1 when T int is short enough, for all f mod . For the longest T int (5 µs), F varies in the range 10 −4 − 10 1 , depending on f mod . After re-scaling the horizontal axis to the period of the signal, panel (b), we see that, for some modulation frequencies, F dips sharply when T int /T mod = n with n = 1, 2, 3, . . . . These minima reveal that, for frequencies that produce locked spikes, the sequence of counts is very regular when the counting interval contains an integer number of periods. In contrast, when the modulation is sinusoidal, Fig. 5a, the first dip occurs for T int /T mod = 2 and no modulation frequency produces a dip at T int /T mod = 1 , confirming that a small-amplitude sinusoidal modulation does not generate harmonic locking.
In Figs. 4 and 5a we also note that for some frequencies F grows with T int as a power law. This behavior reveals spike clustering (bursts of spikes) and has been observed in sequences of events, in different fields [56][57][58] .
To determine whether the F values are due, in part, to "dynamic" properties of the spike sequence (the presence of temporal correlations between the spike times), or they are only due to "static" properties (the shape of the ISI distribution) we shuffled the sequence of ISIs, recalculated the spike times 59 , re-calculated F, and compared the F values of the original and shuffled spike sequence. The presence of long-range order in the timing of the spikes is detected by F values that are significantly higher or lower in the original spike sequence than in the shuffled one. The high or low variability of the sequence of counts reveals the presence of temporal correlations between spikes, which are washed out when we randomly shuffle the values in the ISI sequence.
The Fano factor of the shuffled spikes is displayed in Fig. 5b, where we see that the dips are less pronounced with respect to those in the original spike sequence, Fig. 5a, and we also note that the power law grow disappears in the shuffled data. These differences reveal the presence, in the original spike sequence, of long-range correlations in the spike timing, which are removed after shuffling the ISIs.
The presence of long-range correlations in the spike timing is intriguing at first sight, because near threshold the spiking dynamics of the un-modulated laser is irregular, and we modulate the laser with a small amplitude signal (for I dc = 26 mA the modulation produces a variation of the laser current <3%). However, the delayed feedback introduces a degree of periodicity, and particular combinations of the feedback and modulation parameters can induce spikes that rigidly lock to the modulation with high timing regularity. www.nature.com/scientificreports/ For most frequencies that generate locked spikes, the dips for T int = nT mod gradually disappear when T int increases, revealing that the spike timing regularity does not persist over long intervals, i.e., the spikes do not remain in phase. However, for particular modulation frequencies F decreases with T int as T −1 int , and the dips for T int = nT mod persist for large n. Figure 6 shows two examples of locked spikes, with and without long-range regularity, generated by sinusoidal modulation with frequency 25 MHz and 23 MHz, respectively. As shown in panels (a) and (b), differences can not be distinguished when the intensity dynamics is inspected in a short time interval. In both panels we see that the spikes are periodic but the intensity dynamics is not, because in between spikes, the oscillations are irregular (a similar behavior -regular spikes and irregular oscillations in between spikes-was observed in Fig. 1b). However, if we examine the dynamics over longer time intervals, the spike timing regularity disappears when f mod = 23 MHz, while it persists when f mod = 25 MHz. This difference is seen in the plots of the Fano factor, panels (c) and (d), and also, in the space-time representation of the intensity time series, (e) and (f).
To demonstrate that the regularity of the spike timing generated by sinusoidal modulation with frequency 25 MHz is indeed very long range, we calculate the Fano factor considering longer counting intervals (so far we divided the intensity time series in 1000 intervals of 5 µ s each). To increase the length of the intervals, we need to decrease the number of intervals. The results are shown in Fig. 7a, where we see that F continues decreasing as T −1 int , even when we calculate F using intervals of 500 µ s each (each interval contains 12,500 modulation cycles). In the shuffled spike sequence, Fig. 7b, we see that for large T int the decrease of F saturates and the dips disappear. This confirms a long-range timing regularity in the original spike sequence, which is removed in the shuffled spike sequence, where F values are low just because of the narrow shape of the ISI distribution. This long-range suppression of fluctuations in the timing of the spikes on time scales that contain thousands of modulation cycles, together with the aperiodic nature of the oscillations in between spikes and the lack of 1:1 locking, allows us to interpret this spiking behavior as a non-trivial form of time-crystal dynamics.

Discussion
Summarizing, we have found that near threshold a weak modulation of the laser current induces, for particular modulation conditions, long-range regularity in the timing of the spikes, as revealed by Fano factors that decrease with the length of the spike counting interval as T −1 int . With pulsed modulation we found harmonic and subharmonically locked spikes, but with sinusoidal modulation, we found only subharmonic locking, which is a common feature with time-crystal behavior. We have interpreted the difference in the laser response to pulsed and sinusoidal modulation as due to (1) the existence of the lasing threshold and (2) the fact that harmonic locking occurs for low modulation frequencies. In these conditions the pulsed modulation decreases sharply the laser current, but only for short time intervals, inducing spikes with well-defined periodicity. In contrast, the sinusoidal modulation decreases the current gradually, and this smooth variation brings the laser current close to threshold for longer intervals, allowing for a larger influence of noise and preventing the generation of harmonically locked spikes. Therefore, we have found, in a stochastic dynamical system, discrete time-crystallike behavior generated by the interplay of periodic modulation and delayed feedback.
Regarding the level of noise, compared to other types of lasers, semiconductor lasers are more affected by quantum spontaneous emission because of the phase-amplitude coupling of the optical field (due to the variation of the semiconductor refractive index with the carrier population, an effect that is phenomenological modeled by Henry's alpha factor 38,60 ). Other sources of noise include electric, thermal and mechanical fluctuations. We have not measured the level of noise in our system, but we have performed extensive simulations of the dynamics of the laser with optical feeedback and sinusoidal current modulation using the LK model 37 and found a very   Our findings open a new path for studying the emergence of long range regularity in disordered systems, exploiting the analogy between time-delayed systems and spatially-extended ones. This analogy can be extended to the transverse section of a broad area laser, where, under optical injection, localized structures (that can be regarded as spikes) are created and move as the result of the combined action of nonlinearity, diffraction, carrier diffusion and delayed feedback 63,64 .
The harmonically locked spiking behavior with long-range regularity generated by pulsed modulation may be consider analogous, in the temporal domain, to hyper-uniform states in multi-particle systems, which are characterized by long length scale suppression of density fluctuations 65 . Moreover, because delayed feedback is a popular control technique with multiple applications, we expect that our observations will motivate studies to generate, characterize, and exploit the highly regular oscillations that can be generated by the interplay of nonlinearity, feedback and modulation.

Methods
To compute the Fano Factor we first divide the intensity time trace in N int non-superposing segments of duration T int and count the spikes in each segment; then, we calculate F = σ 2 (N i )/ �N i � , where σ 2 (N i ) and N i are the variance and the mean of the sequence of counts, {N i , i = 1 . . . N int } . F is a function of T int ; if T int is very small, F = 1 because the sequence of counts is a sequence of 0s and 1s, while in time scales over which the spikes are regular, the variance of {N i } is small and F takes low values. Unless specifically stated, to calculate the Fano factor we divide the intensity time-series in N int = 1000 non-overlapping segments of T int = 5 µ s each.
If the spikes are generated by a fully random, homogeneous Poisson point (HPP) process, the probability that N spikes occur in an interval T int is p(N, T int ) = ( T int ) N exp(− T int )/N! where is the average number of spikes per unit time (i.e., the spike rate). For this distribution σ 2 (N i ) = �N i � = T int and thus, F = 1 ∀ T int 56 .

Data availibility
The experimental sequences of inter-spike-intervals are available here https:// doi. org/ 10. 5281/ zenodo. 59135 06. 59. The original spike times are t i = t 0 + i j=1 T j where t 0 is the time where the first spike occurs, and T j = t j − t j−1 is the time interval between spikes j − 1 and j (i.e., the inter-spike interval, ISI). The sum t 0 + N j=1 T j , where N is the total number of spikes, is equal to the recorded time, 5 ms. By shuffling the sequence of ISIs we obtain a new sequence of uncorrelated ISIs, { T ′ j } (with the same distribution of values, that also satisfies t 0 + N j=1 T ′ j = 5 ms), from which we calculate the shuffled spike times as t ′ i = t 0 + i j=1 T ′ j . 60. Henry, C. H. Theory of the linewidth of semiconductor lasers. IEEE J. Quantum Electron. 18, 259 (1982).