Article https://doi.org/10.1038/s41467-025-58400-6 Exploring nonlinear dynamics in periodically driven time crystal from synchronization to chaotic motion Alex Greilich 1,4 , Nataliia E. Kopteva 1,4 , Vladimir L. Korenev2, Philipp A. Haude1 & Manfred Bayer 1,3 The coupled electron-nuclear spin system in an InGaAs semiconductor as testbed of nonlinear dynamics can develop auto-oscillations, resembling time- crystalline behavior, when continuously excited by a circularly polarized laser. We expose this system to deviations from continuous driving by periodic modulation of the excitation polarization, revealing a plethora of nonlinear phenomena that depend on modulation frequency and depth. We find ranges in which the system’s oscillations are entrained with the modulation fre- quency. The width of these ranges depends on the polarization modulation depth, resulting in an Arnold tongue pattern. Outside the tongue, the system shows a variety of fractional subharmonic responses connected through bifurcation jets when varying the modulation frequency. Here, each branch in the frequency spectrum forms a devil’s staircase. When an entrainment range is approached by going through an increasing order of bifurcations, chaotic behavior emerges. These findings can be described by an advanced model of the periodically pumped electron-nuclear spin system. We discuss the con- nection of the obtained results to different phases of time matter. Nonlinear systems are characterized by significant interactions between their constituents, resulting in collective dynamics1. A pro- minent example of nonlinear effects is frequency synchronization, a phenomenon observed for living organisms, like for heartbeats or cricket chirping, as well as for nonliving systems, such as clocks and generators. Synchronization refers to the adjustment of the fre- quencies of autonomous systems, which periodically oscillate for constant external excitation, due to their weak interaction. Particularly interesting is the synchronization of the system’s auto-oscillations with an external periodic source, the frequency of which is close to one of the frequencies in the unperturbed har- monic spectrum. Then, even for a weak influence, significant chan- ges occur, such as the adjustment of the circadian rhythms in organisms to the day-night cycle or the control of the pacemakers regulating heart rhythms. Despite their seeming differences, these synchronization phenomena can be understood using the common framework established for the nonlinear dynamics of complex systems1–3. Implementing synchronization in semiconductors is an intriguing possibility because of the application relevanceof thesematerials, e.g., in modern electronics. Previously, we revealed highly robust, non- decaying auto-oscillations in the electron-nuclear spin system (ENSS) of a tailored semiconductor4 as a clear signature for a continuous time crystal. The concept of time crystals was introduced in 2012 by Frank Wilczek5,6 for a closed many-body system in thermodynamic equili- brium. While this original time crystal phase is prohibited7–9, time crystalline behavior demonstrating spontaneous breaking of the time translation symmetry was predicted to be feasible in non-equilibrium. Such systems offer an attractive avenue for high-precision exploration of non-equilibrium states of matter. Received: 23 July 2024 Accepted: 19 March 2025 Check for updates 1Experimentelle Physik 2, Technische Universität Dortmund, Dortmund, Germany. 2Ioffe Institute, St. Petersburg, Russia. 3Research Center FEMS, Technische Universität Dortmund, Dortmund, Germany. 4These authors contributed equally: Alex Greilich, Nataliia E. Kopteva. e-mail: alex.greilich@tu-dortmund.de; natalia.kopteva@tu-dortmund.de Nature Communications | (2025) 16:2936 1 12 34 56 78 9 0 () :,; 12 34 56 78 9 0 () :,; http://orcid.org/0000-0001-7813-682X http://orcid.org/0000-0001-7813-682X http://orcid.org/0000-0001-7813-682X http://orcid.org/0000-0001-7813-682X http://orcid.org/0000-0001-7813-682X http://orcid.org/0000-0003-0865-0393 http://orcid.org/0000-0003-0865-0393 http://orcid.org/0000-0003-0865-0393 http://orcid.org/0000-0003-0865-0393 http://orcid.org/0000-0003-0865-0393 http://orcid.org/0000-0002-0893-5949 http://orcid.org/0000-0002-0893-5949 http://orcid.org/0000-0002-0893-5949 http://orcid.org/0000-0002-0893-5949 http://orcid.org/0000-0002-0893-5949 http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58400-6&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58400-6&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58400-6&domain=pdf http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58400-6&domain=pdf mailto:alex.greilich@tu-dortmund.de mailto:natalia.kopteva@tu-dortmund.de www.nature.com/naturecommunications During the last decade, two distinct scenarios have been devel- oped: continuous (CTC) and discrete (DTC) time crystals, where the external source is acting in a steady or modulated way, respectively10–12. For excitation without modulation, dissipative CTCs become possible, where the continuous energy flow from an external source is transformed into oscillatorymotion in so-called autonomous systems. This type of behavior has been discovered across manifold different physical systems ranging from Bose-Einstein-Condensates13, Rubidium atom condensates14, Erbium ion ensembles15, strongly interacting Rydberg gases16, light fields in photonic metamaterials17, polariton condensates in semiconductors18, to electrons in electro- statically defined double quantum dots19 and mesoscopic fractional quantum Hall devices20,21. For modulated excitation, DTCs were experimentally confirmed in various systems by demonstrating sub- harmonic responses to periodic excitation.22–31. Here, we use the robust CTC platform implemented in a semi- conductor ENSS to experimentally and theoretically explore the non- linear dynamics achieved by deviating from continuous driving through periodic modulation of the exciting laser polarization. The deviation causes a wide variety of phenomena ranging from synchro- nization to chaoticmotion that depend on the frequency and depth of the modulation. The synchronization becomes evident through fre- quency entrainment, where the ENSS synchronizes with the modula- tion frequency or its rational fractions. The width of the entrainment frequency bands depends on the modulation depth, leading to an Arnold tongue pattern32,33. Furthermore, we observe frequency bifur- cation jets at the entrainment edges. Across the different ranges, each frequency branch in the system’s response corresponds to a devil’s staircase32. The entirety of the branches forms a fractal pattern show- ing self-similarity, confirmed by simulations. Going through an increasing number of bifurcations, we see a transition to chaotic behavior before reaching entrainment again. The entirety of different phenomena observed in a single system underscores that this system is an ideal testbed for the nonlinear dynamics of nonequilibrium sys- tems, particularly for understanding synchronization phenomena. It is essential to highlight the difference between the synchroni- zation observed in nonlinear autonomous systems and another well- known phenomena, resonance and parametric resonance. The latter can be observed in periodically driven systems that do not demon- strate self-sustained oscillations. So, there is no synchronization, as the system does not have its own rhythm, and if the periodic driving is interrupted, the oscillation stops after some transition time2. The effect of parametric resonance or subharmonic oscillations, initially observed by Faraday in periodically driven granular material34,35, represents a fundamental mechanism where a system oscillates at a fraction of the driving frequency. Landau and Lifshitz later formalized this phenomenon in the context of parametric resonance36, highlighting its universal occurrence. Unlike the classical Faraday instability, our results arise in a nonlinear auto-oscillating system, demonstrating a synchronizationof higher orders andoffering new insights into subharmonic dynamics. Results Periodic auto-oscillations For our studies, we have used the same ternary semiconductor struc- ture, with which the CTCwas implemented in a reservoir of In, Ga, and As nuclear spins4. The structural parameters are: A 10μmthick epilayer of In0.03Ga0.97As is doped with Si donors providing an electron con- centration of 3.9 ×1016 cm−3 (see inset in Fig. 1a andMethods section for further details). The In atoms are incorporated to induce isotropic strain that causes a nuclear quadrupole splitting required for CTC operation. The reservoir also requires optical pumping because of its open system character. The resulting losses of angular momentum, e.g., due to nuclear dipole-dipole interaction, can be exactly com- pensated by circularly polarized laser excitation, which orients the localized electron spins of the donors, fromwhich access is granted to the nuclear reservoir via the hyperfine interaction. This interaction within the electron localization volume maintains the nuclear spin polarization via flip-flop processes, establishing a robust CTC in the illuminated volume. The photoluminescence of this structure at low temperature (T = 6 K) is governed by the emission from free and donor-bound excitons; see the blue curve in Fig. 1a. Using circular polarization of a pump laser at Epu = 1.579 eV photon energy, non-zero electron spin polarization is created, which can be assessed by measuring the Faraday rotation (FR) of a linearly polarized continuous wave probe laser. The FR spectral dependence is shown by the red trace in Fig. 1a. In all further experiments, the probe energy is fixed at Epr = 1.454 eV, corresponding to optimized experimental conditions: The FR has a local maximum with weak probe laser absorption. Measuring the Faraday rotation by the electron spins also gives us access to the nuclear polarization, which is imprinted in the electron spins through the effective nuclear magnetic field (Over- hauser field). Probe 1.43 1.44 1.45 1.46 1.47 0 1 Energy (eV) P ho to lu m in es ce nc e F ar ad ay r ot at io n T = 6 K a FR PL GaAs substrate 10 µm In0.03Ga0.97As ne = 3.9 × 1016 cm-3 b 0 50 100 150 200 1 2 Time (s) F R ( m ra d) 0 0.2 0.4 0.6 0.8 1 1.2 0 1 Frequency (Hz) |F F T | x zBext c f0 nat 2f0 nat 3f0 nat 4f0 nat 5f0 nat 6f0 nat Fig. 1 | Optical properties and auto-oscillations. a Photoluminescence spectrum at T = 6 K, excited by a continuouswave diode laser at Epu = 1.579 eV photon energy (blue line). The red trace shows the Faraday rotation spectrum,measured by tuning the photon energy of the continuous wave probe laser. The inset shows a sketch of the Si-doped In0.03Ga0.97As layer with resident electrons concentration of ne = 3.9 × 1016 cm−3 on top of a GaAs substrate. b Periodic auto-oscillations mea- sured in a tilted magnetic field with components Bx = − 1 mT and Bz = 0.176 mT, corresponding to the tilt angle α = 10∘. The time trace recording started after five minutes of pumping for transients tohave faded. Pumpandprobephoton energies: Epu = 1.579 eV, Epr = 1.454 eV. Pump and probe powers: Ppu = 0.3 mW, Ppr = 1 mW. c Fast Fourier transform of the signal in panel (b), recorded for 10 min. The inset shows a sketch of the experimental geometry with the magnetic field orientation relative to the sample. Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 2 www.nature.com/naturecommunications For specific pumping conditions (continuous, non-modulated excitation with, for example, magnetic field components Bx = − 1 mT and Bz = 0.176mT, laser powers Ppu = 0.3 mW and Ppr = 1 mW, T = 6 K), the ENSS demonstrates robust non-decaying auto-oscillations - a clear CTC manifestation, see Fig. 1b, as studied in detail in ref. 4. The fast Fourier-transformed (FFT) spectrum of these oscillations is shown in Fig. 1c, corresponding to a CTC structure analysis like for a space crystal, where X-ray inspection gives the Fourier transform of the electron charge density in a crystal lattice. The spectrum is harmonic with the ground (natural) frequency f nat0 =0:1607 Hz and the higher harmonics nf nat0 , where n is an integer. Periodic modulation of excitation polarization We apply the modulation to the pump laser polarization at the fre- quency fm, and keep the excitation power constant. The degree of circular polarization is varied in the following way: ρc = 1 2 ½2� ε 1� cosð2πfmtÞ � ��, with ε being the modulation depth. For ε=0, the light polarization isfixed at circular; for ε= 1, there is variation between circularly and linearly polarized pumping, while for ε < 1, the light is modulated between circular and elliptic polarization of a cer- tain degree, see Supplementary Fig. 1b. Figure 2a demonstrates exemplary time traces for fm varied starting from proximity to the f nat0 of the unperturbed CTC case to higher frequencies with ε = 1, where otherwise the same experimental conditions were applied as for the auto-oscillations in Fig. 1b. The corresponding FFT spectra in Fig. 2b, taken for time traces with 10- minute recording time, reveal a non-trivial behavior. The red trace demonstrates the case for fm = 0.1607 Hz, close to the basic CTC harmonic of the autonomous ENSS. The spectral positions of the FFT peaks coincide with those of the unperturbed case shown by the gray spectrum behind the red one in Fig. 2b, but the amplitudes differ. When fm increases, the FFT harmonics shift rigidly with the applied modulation frequency. This behavior is well known from periodically driven dissipative systems, titled frequency locking37 or entrainment2,38. In our case, the whole frequency spectrum f exp of the system is entrained with the modulation frequency fm. This situation continues up to the frequencyof fm=0.197Hz, where rather abruptly a multitude of subharmonics appears; see the bottom blue spectrum in Fig. 2b. The frequency range across which the spectrum is locked to the modulation frequency shift from the unperturbed spectrum without any appearance of additional subharmonics is called the entrainment range.Wenormalize the FFT spectraby fm tohighlight the frequency entrainment. Then, the peaks do not change their spectral position across the entrainment range, see Fig. 2c. It is important to note that f nat0 is the system’s natural frequency without modulation. If modulation is applied, the pump’s average helicity is changed, renormalizing the natural frequency to f0. In the Supplementary Fig. 2, we additionally show time traces and corresponding FFT spectra for fm that match several rational fractions of f0 to fm in the range from f0/fm = 1/2 to f0/fm = 1. The FFT spectra demonstrate subharmonic frequencies, which divide the spectra into several equidistant intervals. The number of these intervals for f exp=fm ≤ 1 is equal to the denominator of f0/fm. Close to these rational fractions, one observes frequency entrainment. As an example, Fig. 3a presents the FFT spectra for f0/fm varied around 1/2, from 0.477 to 0.539, for ε = 1, evidencing a range of entrainment with two FFT peaks in the shown frequency interval: at the modulation frequency fm � 2f nat0 , close to the second harmonic of the unmodulated ENSS, and at half of the modulation frequency representing a subharmonic 0 50 100 150 0 1 2 3 4 5 Time (s) F R ( m ra d) 0.197 a b c 0 0.1 0.2 0.3 0.4 0.5 Frequency, fexp (Hz) |F F T | fm (Hz) 0.190 0.172 0.1607 no mod = 1 0 0.5 1 1.5 2 2.5 fexp / fm |F F T | f0 / fm 1.082 1.011 0.915 0.883 Fig. 2 | Periodically driven ENSS. a FR signal of periodically driven oscillations at frequencies deviating increasingly from thebasic eigenfrequency,first startingwith this frequency fm = 0.1607 Hz (red), and then increasing to fm = 0.172 Hz (green), fm = 0.190 Hz (brown), fm = 0.197 Hz (blue), all for ε = 1. The other experimental parameters are as in Fig. 1. The trace recording starts after five minutes of waiting time with the pump on to have transient effects suppressed. b FFT spectra of the signals from panel (a), recorded for 10minutes. The gray trace reproduces the FFT of the unmodulated signal with the basic harmonic f nat0 = 0:1607 Hz. c Same FFT spectra as in panel (b) with the frequency axis normalized by fm, so that the FFT peaks are held at constant positions in the case of entrainment. The labels at the curves give f0/fm. Note that f0 is affected by the average helicity value. 0.48 0.5 0.52 0.54 0 0.25 0.5 0.75 1 f0 / fm 0 0.5 1 fexp / fm |F F T | a b 0.477 0.484 0.490 0.497 0.526 0.539 f0 / fm= = 1 Fig. 3 | Entrainment range andArnold tongue. a FFT spectra for fm varied around twice the ENSS basic harmonic (2f nat0 ) with full modulation depth ε = 1. Here, the frequency scale is also normalized to the modulation frequency fm. b Points mark the borders of the entrainment range (blue area) in dependence on themodulation depth, forming an Arnold tongue. The borders are set at the half-step size between the synchronized case with clean FFT peaks and the unsynchronized case with multiple side peaks in units of f0/fm. The error bars are correspondingly given by half-step sizes in the same units. Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 3 www.nature.com/naturecommunications response. As discussed in the introductory paragraphs, observing a stable subharmonic response with entrainment to the modulation frequency is generally considered a signature of a DTC state, as will be discussed further below. The blue FFTs with multiple subharmonics in Fig. 3a indicate the boundaries of the entrainment range. At the edges of all entrainment ranges, the FFT spectra consist of multiple peaks around the main harmonics. Exemplarily, we discuss the edge at f0/fm = 0.883, shown in Fig. 2c, and analyze the spectra taken in small frequency steps to characterize the transition across the edge carefully, see Fig. 4a. We chose this edge due to a broader range on the low-frequency side of modulation that is not obscured by additional subharmonics in the entrained state which reduces the amplitudes of contributing frequencies. Here, at the very edge of the transition to entrainment, the number of FFT peaks drastically increases, see the spectrum for f0/ fm = 0.886 with the corresponding time trace in Fig. 4b. The merging FFT peaks are a signature of chaotic oscillations. To confirm this, we perform several chaos tests established in the literature. To that end, a 30-min recording time trace is measured and ana- lyzed. The auto-correlation function shows slowly decaying beats with increasing delay. The beats appear due to the superposition of oscil- lating contributions with frequencies f0, fm, and commensurate har- monics of both of these39. The slow decay is an indication of deviation from periodic behavior. A nonlinear time series analysis gives the positive maximal Lyapunov exponent λmax =0:06 40, evidencing an exponential deviation of closely spaced initial states with time, and a non-integer correlation dimension ofD2 = 2.541, confirming the chaotic behavior42 (see Supplementary Fig. 3a–c). For comparison,we alsogive as an example the same quantities for the synchronization plateau at f0/fm = 2/3, resulting in λmax =0 and D2 = 1, see the Supplementary Fig. 3d–f, confirming a near ideally periodic behavior. Additionally, there are ranges between the synchronization plateaus where the correlation dimension has the value of D2 = 2, which represents quasi- periodicity with an irrational ratio among the observed frequencies (see the Supplementary Fig. 3g–i) and corresponds to an intermediate situation where the system is neither periodic nor chaotic42. Finally, to reveal the complete picture of the ENSS response to frequencymodulation, we provide FFT spectra for a range of fm varied from the basic harmonic f nat0 up to slightly more than twice this fre- quency. Here, we use the inverted frequency scale again so that the abscissa interval covers 0.45 ≤ f0/fm ≤ 1, for ε = 1. The resulting mea- surement series is shown as a color map in Fig. 5a (and in Supple- mentary Fig. 7 in stretched format for more details) and demonstrates a multitude of frequencies in the subharmonic regime below the applied modulation frequency fm. We want to highlight some points introduced before to interpret the observed map. The horizontal plateaus of finite width seen at rational values of f0/fm (see top axis in Fig. 5a) represent entrainment ranges. The length of the plateaus decreases with increasing denomi- nator. So, the plateaus at the ratios 1/2, 3/5, 2/3, 3/4, and 4/5 have lengths of 0.059, 0.012, 0.018, 0.013, and 0.010 ± 0.001 in units of f0/ fm, correspondingly. Leaving an entrainment range by moving to smaller fm, we observe a splitting of each frequency into many bran- ches, comprising bifurcation jets38. Moving further and approaching the next entrainment range, the number of bifurcations continuously increases and becomes so high that eventually, chaotic behavior emerges through the coupling between the different resonances caused by the nonlinearities32. Model of periodically modulated ENSS To obtain additional insight into the experimentally discovered behavior, we have generalized the model described in refs. 4,43 for an autonomous system by adding the periodic force variation provided by polarization modulation at frequencies close to the intrinsic auto- oscillation. We briefly repeat the basic features: The circularly polar- ized pump excitation orients the donor electron spins, which subse- quently polarize the nuclear spin system via the hyperfine interaction44. TheOverhauser field of the polarized nuclear spinsBN, in general, is oriented not parallel to the average electron spin S, so that an electron spin precesses about BN, causing a variation of S. Thus, in the strongly coupled nonlinear ENSS system, the electron spins and the Overhauser field are mutually interdependent via their magnitude and direction. Then, for continuous pumping, the dynamic regime of self-sustained auto-oscillations appears under specific conditions. Due to the short electron spin lifetime (Ts ~ 1 μs) compared to the longitudinal nuclear spin relaxation time (TN ~ 1 s), the electron spin (S) is described by the solution of the stationary Bloch equation accounting for the sum of the external magnetic field (Bext) and the Overhauser field (BN): S=S0 + μBgT s _ ðBext +BNÞ× S: ð1Þ Here, S0 is the average electron spin polarization induced by the pump without magnetic field, μB is the Bohr magneton, ℏ is the reduced Planck constant, and g is the electron g-factor. 0 0.5 1 1.5 fexp / fm |F F T | 0 50 100 150 200 250 300 1 1.5 Time (s) F R ( m ra d) f0 / fm = 0.880 0.885 0.889 0.899 f0 / fm = 0.886 a c b 0 100 200 300 400 500 −0.5 0 0.5 1 Delay (s) A C 0.886 Fig. 4 | Chaos at edge of entrainment range. a FFT spectra for fm varied close to the ENSS basic harmonic f nat0 for full modulation depth ε = 1. Again, the abscissa is normalized to the modulation frequency. b Faraday rotation signal of periodically driven oscillations with frequency f0/fm = 0.886. c Autocorrelation (AC) function calculated for the signal in (b). The green dashed line highlights the AC decay. Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 4 www.nature.com/naturecommunications The precession of the electron spin polarization about the total magnetic field (Bext + BN) changes the Overhauser field in time according to44,45: dBN dt = � 1 TN BN � âS � � , ð2Þ where â is a second-rank tensor describing the process of dynamic nuclear polarization. Themodel suggests that âS is a linear function of S. The tensor has been simplified in the lowest order approximation (for details, see Methods section, ref. 4 and corresponding Supple- mentary Material). This allowed us to simulate the CTC auto- oscillations, as shown in ref. 4. For the calculations in this paper, we use the same experimental parameters as in the previous work: α = 10∘, Bx = − 1mT, effectivemagneticfields ofaN = 20mT and bN = 21mT, and nuclear spin relaxation time of TN = 0.5 s. We extend the model by implementing the modulation of laser polarization, leading to the periodicity of S0 with the same frequency: S0 follows the degree of circular polarization ρc so that we substitute it in Eq. (1) by: S0,m = S0 2 ½2� εsimð1� cosð2πfmtÞÞ�, ð3Þ where εsim is the modulation depth in the simulation, allowing a deviation from the experimental ε. The analysis of the dynamical Eqs. (1)–(3) reveals several new results, which are described below. We simulate the electron spin in the presence of the nuclear field toobtain themeasured signal for varyingmodulation frequency fmand calculate the dependence of the FFT using 20-minute time traces, which evolve to formboth frequency entrainment ranges and complex bifurcation patterns between them, as presented in Fig. 5b. The Sup- plementary Fig. 4 additionally presents the simulated maps for dif- ferent depths of modulation εsim. Overall, we find good agreement between calculation and experiment, but with an adjusted value of εsim =0:5 compared to ε = 1 in the experiment, which can be explained by several factors, including strongly non-resonant pump excitation, resulting in partial electron spin relaxation. Again, we plot the spectra in a contour as a function of the modulation rate f0/fm33. As mentioned, using this coordinate, the entrainment branches are identified as horizontal lines. An additional advantage of this presentation is the possibility to introduce a char- acteristic parameter, the so-called effective winding number w = fcirc/ fm= Tm/Tcirc, which is the ratio of themodulation period Tm to the time Tcirc that the ENSS takes to return to its initial point duringmotion on a closedperiodic trajectory (limit cycle). This quantitywas introduced as a simplified version of the winding number in the circle map repre- sentation of refs. 33,42. To understand its connection with the observed spectra, let us consider the entrainment plateau in Fig. 5b around f0/fm = 1 on the x-axis and at f exp=fm = 1 on the y-axis. In this case, the system under- goes one full revolution along the limit cycle during one period Tm of modulation, giving the effective winding numberw = 1. Further, for the entrainment plateau around f0/fm = 1/2 on the x-axis and at f exp=fm = 1=2 on the y-axis the system undergoes one full revolution along the limit cycle during two periods Tm of modulation, giving the effective winding number w = 1/2. Larger values of f0/fm > 1/2 asso- ciated with entrainment plateaus have multiple frequency peaks for f exp=fm ≤ 1, where only the lowest frequency corresponds to w = fcirc/fm. Our calculation shows that near resonances f0/fm = M/N the effective winding number w = 1/N, where M, N are mutually prime (coprime) natural numbers33. The Supplementary Fig. 2a demonstrates the period Tcirc for several cases, representing the lowest harmonic frequency in the corresponding FFT spectra in the Supplementary Fig. 2b. We additionally show several examples of simulated spin tra- jectories giving all vector components, which highlight the spin evo- lution for the cases Tcirc = Tm, 2Tm, 4Tm, see the Supplementary Fig. 5a–c, respectively. In themeasured spectra of Fig. 5a, there are also other harmonics, so that the whole series of observed frequencies is given by f exp=fm =K=N, for K = 1, 2, . . . . Remarkably, the entrainment ranges associatedwith the numbers f exp=fm =M=N = 1=N or M/N = (N − 1)/N in dependence of f0/fm form devil’s staircases46,47. The step width in each staircase becomes the smaller, the larger the denominator N32, following the so-called Farey tree sequence, eventually leading to a chaotic state42. Our simulation shows that such ranges with a very high number of bifurcations appear close to each transition to synchronization, see Fig. 5b, c. Due to the redistribution of the frequency components between multiple sub- harmonics, these are much harder to observe experimentally. 0.5 0.6 0.7 0.8 0.9 1 1/2 3/5 2/3 3/4 1/14/5 1/1 1/2 1/3 2/3 1/4 3/4 1/5 2/5 3/5 4/5 f0 / fm 0.61 0.62 0.63 0.64 0 0.1 0.2 0.3 f0 / fm simulation a b c f0 / fm f0 / fm 0 0.5 1 0.5 0.6 0.7 0.8 0.9 1 1/2 2/3 3/43/5 4/5 0 0.2 0.4 0.6 0.8 1 f0 / fm f ex p / f m experiment Fig. 5 | Entrainment, bifurcations, and fractal structure. a Contour plot of experimental FFT spectra with the frequency axis f exp normalized by the modula- tion frequency fm as a function of inverse modulation frequency, multiplied by the basic harmonic frequency. ε = 1. For each FFT, time traces of 10 minutes of recording time are used.b Sameas (a), butwith FFT spectra fromsimulations, using the parameters from Fig. 2c and εsim =0:5. To highlight the fractal nature of the signal, panel (c) shows a zoomof the areamarkedby the light bluebox shown in (b). The scale on the right shows the color scheme for the normalized amplitude of the contour maps. Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 5 www.nature.com/naturecommunications The observed devil’s staircases evidence the self-similar character of fractal structures, i.e., each detailed section, no matter how small, contains the same features as the whole picture. Figure 5c shows a zoom of the simulated part marked by the light blue-dashed box in Fig. 5b, basically reproducing the large area and confirming the self- similarity. A similar behavior, in terms of devil’s staircase structures, was recently demonstrated for frequency-locked breathers in ultrafast lasers48 and theoretically proposed for the optomechanical locking in driven coupled polariton condensates49, confirming the universal nature of the observed non-linear phenomena. Furthermore, reducing the modulation depth ε leads to the nar- rowing of the observed entrainment ranges. To confirm this, we determine the edges of the experimentally observed entrainment range close to f0/fm = 1/2 and plot it as a function of ε by the symbols in the phase diagram in Fig. 3b. Below ε = 0.25, the entrainment vanishes so that the FFT spectra resemble that of the unmodulated ENSS, as the deviations from the circular pumping become too small. This demonstrates the remarkable stability of the ENSS for CTC operation with respect to variations in the laser helicity (ε <0.2). The dependence of the width of the entrainment range on themodulation depth results in the blue shaded area in Fig. 3b, which is known as Arnold tongue2,50. Additionally, the Supplementary Fig. 6 shows the dependence of the auto-oscillations on the temporally constant degree of pump polar- ization, having maximal amplitude for circular polarization while dropping to zero for linear polarization. In ref. 32, the authors studied the scaling law underlying the Arnold tongues. Similar to their results, the observed fractal structure, in our case, is related to the structure of the self-similar Cantor set. Using the experimentally observed Farey tree sequence with major plateaus and gaps between them, we arrive at the fractal dimension value of D0 =0:852±0:021, while the simulated data, with the possi- bility to zoom into the staircase structure, delivers D0 =0:853 ±0:002, seeMethods. These values are close to 0.87 expected for the complete devil’s staircase46 and align with the universal properties of mode- locking transitions and the devil’s staircase structure seen in dis- sipative systems51 as well as in refs. 48,52,53. This suggests that the ENSS likely belongs to a broader universality class, sharing con- vergence and scaling properties akin to those found in circlemaps and forced oscillators. To round out the picture, the Supplementary Fig. 4 shows simu- lated colormaps for decreasing εsim, demonstrating how in parallel the entrainment ranges narrow and the bifurcation structures diminish, leading to the disappearance of the devil’s staircases. They also con- firm the Arnold tongue structure in Fig. 3b. On the other hand, with increasing modulation depth (εsim>0:5), the entrainment ranges (or neighboring Arnold tongues) start to overlap, also leaving no space for bifurcations. This situation cannot be reached for pump polarization modulation in the case of non-resonant excitation, but can be realized for modulation of another pump parameter, namely the pump power, as will be analyzed elsewhere. Discussion Experimental results are well described within the framework of non- linear physics, which can and needs to be connected further to the physics of time crystals in nonlinear systems. As discussed in the introductory paragraphs, observing a stable subharmonic response with entrainment to the modulation across a finite frequency range, like for the ENSS here, is generally related to a DTC state. The entrainment range represents the area of stability for such a state. Further, the emergence of multiple subharmonic resonances at spe- cific fractions f0/fm, as observed in the Supplementary Fig. 2b, has been discussed in termsof fractional andhigher-orderDTC states54–59. In this regard, the CTC state, in our case, undergoes the transition to a DTC phase,where the subharmonic responseoccurs on rational fractions of the modulation frequency fm, as well described by our model. Recently, a similar transition was demonstrated for an atom-cavity system60. Considering all these findings, our DTC phase belongs to a specific class of nonlinear dynamic systems. The discussion in terms of time crystals is generally applicable to strictly periodic phenomena in autonomous systems on a limit cycle and non-autonomous systems in the synchronization regime. Our system also features aperiodic modes of nonlinear systems, such as chaotic oscillations and fractal structures with quasi-periodic oscilla- tions in the subharmonic range. The former system phase may be calledmelted time crystal or time glass and the latter phase timequasi- crystal, similar to recent experimental work on a strongly interacting spin ensemble in diamond61. However, so far, the definition of time matter is not fully clear-cut within the general frame of dynamics of nonlinear autonomous and non-autonomous oscillating systems or whether it even goes beyond that frame. The connection between our findings and those in Rydberg gases57–59, Bose-Einstein-Condensates of Rubidium atoms60, exciton- polariton systems49, and spins in diamond61, underscores the universal mechanisms underlying time-crystalline behavior across different physical systems. Our solid-state platformoffers newopportunities for further studies of nonlinear dynamics. The phenomena studied in this work are primarily driven by the local dynamics arising from nonlinear feedback between the electron and nuclear spins within an ensemble of donors. However, the spatial coupling between donors would introduce a degree of nonlocality, which can lead to spatial synchronization across neighboring donors. While our current analysis focuses on temporal dynamics, future work will explore potential spatial patterns in terms of partial synchronization62–64 or instabilities in the chaotic regime, similar to those described in spatially extended systems65,66. Furthermore, the use of a semiclassical model in our case is jus- tified by the large ensemble of nuclear spins (more than 105 per donor) and multiple orders of magnitude of time-scale separation between the electron and nuclear spin dynamics. However, one could consider detecting quantum correlations among separate single donors, which would require high spatial resolution and a diluted nature of Si doping for their experimental investigation. Finally, the initially observed, unmodulated CTC state of the ENSS represents a limit cycle in the phase space. The feedback strength in this nonlinear system, which keeps it on a limit cycle, determines the oscillation stability (i.e., the quality factor). The results presented here can also be seen from a different perspective, where the external quality factor is mapped on the ENSS through synchronization. One may envision a situation where a highly stable external frequency generator stabilizes the ENSS oscillations, as is the case of quartz resonators in atomic clocks. Such an ultra-stable macroscopic state may inspire new applications in nonlinear physics and quantum technology. Methods Setup Reference 4 and the corresponding supplementary information describe our sample and the experimental setup in all details. Here, we reiterate themost relevant parts.Weusea continuouswave laser diode emitting at 1.579 eV photon energy (785 nm wavelength) as a pump laser, which is then routed through an electro-optical modulator. It changes the phase between the two orthogonally linear polarized light components of the pump in a sinusoidal way, resulting in the desired change of light polarization. The linearly polarized probe laser is cre- ated by a continuous wave Ti:Sapphire ring-laser and is fixed at 1.454 eV photon energy (852.63 nm wavelength). Both lasers are combined on a non-polarizing beam splitter and focused by a single lens onto the sample. The pump laser is completely absorbed by the GaAs substrate of the sample, while the probe laser is transmitted through it and then analyzed by a polarization bridge, which consists of a half-wave plate Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 6 www.nature.com/naturecommunications and a Wollaston prism. A balanced photodiode is used to measure the Faraday rotation of the probe’s linear polarized plane, see Supple- mentary Fig. 1a. The sample is mounted in a helium flow cryostat at temperature T = 6K in the center of two orthogonal pairs of electro- magnet coils, generating the magnetic field components Bx and Bz. Details on data processing To detect the position of f0 in the modulation case, we set the mod- ulation frequency fm to about 3:5f nat0 . In this way, the position of f0 is unaffected by the proximity of the modulating frequency. Details on simulation model The quadrupole unperturbed nuclear spin sublevels create the Over- hauser field B0 N as in pure, unstressed GaAs. B0 N is aligned along the external magnetic field and compensates for the Zeeman splitting of the electrons in the external magnetic field. Due to the strong defor- mation causedby the indium incorporation, the spinof the i-th nucleus is oriented along the main local axis ni of the tensor describing the quadrupole interaction rather than along the external magnetic field. The contribution of these nuclei to the total Overhauser field is BQ = ∑iai(Sni)ni, where the summation is carried out over all quadru- pole perturbed nuclei within the electron localization volume around a donor. For an isotropic distribution of the axes, the field can bewritten as BQ = aNS. Therefore, â can be reduced to the simplified form: âS=BQ +B0 N =aNS +bNðShÞh, where bN is the parameter of the hyper- fine interaction between the electrons and the nuclei, and h is the unit vector of the externally applied magnetic field. The tensor compo- nents are: α̂αβ =aNδαβ +bNhαhβ, with α, β = x, y, z coordinates. Fractal dimension We select out two arbitrary neighboring rotation numbersp0=q0 and p″/ q″ and measure the length of the gap S in between. We know from the Farey tree structure that the largest synchronization interval in between belongs to the rotation number ðp0 +p00Þ=ðq0 +q00Þ. The gaps between the new interval and the twoprevious ones are denoted S0 and S″. According to ref. 67, an approximation D0 for the fractal dimension D can be determined from the relation: ðS0=SÞD 0 + ðS00=SÞD 0 = 1. A better approximation is given for the dimensions by zooming into the steps structure, which is possible using the simulations. In this case, using three successive plateaus, we determine the Sn, S0n, S00n, and use: D0 = limn!1Dn with ðS0n=SnÞ D0 n + ðS00n=SnÞ D0 n = 168. Data availability The data on which the plots in this paper are based and other findings of this study are available from the corresponding authors upon request. Code availability The code onwhich the calculations within this paper are based, as well as other findings of this study, are available from the corresponding authors upon request. References 1. Scott, A. C. The Nonlinear Universe, The Frontiers Collection (Springer, 2007). 2. Pikovsky, A., Rosenblum,M.&Kurths, J.Synchronization. A universal concept in nonlinear sciences (Cambridge University Press, 2010). 3. Schuster, H. G. and Just, W. Deterministic chaos: An Introduction (Wiley-VCH Verlag GmbH&Co. KGaA, 2005). 4. Greilich, A. et al. Robust continuous time crystal in an electron–nuclear spin system. Nat. Phys. 20, 631 (2024). 5. Wilczek, F. Quantum time crystals. Phys. Rev. Lett. 109, 160401 (2012). 6. Shapere, A. & Wilczek, F. Classical time crystals. Phys. Rev. Lett. 109, 160402 (2012). 7. Nozières, P. Time crystals: Can diamagnetic currents drive a charge density wave into rotation? Europhys. Lett. 103, 57008 (2013). 8. Bruno, P. Impossibility of spontaneously rotating time crystals: A no-go theorem. Phys. Rev. Lett. 111, 070402 (2013). 9. Watanabe, H. & Oshikawa, M. Absence of quantum time crystals. Phys. Rev. Lett. 114, 251603 (2015). 10. Sacha, K. & Zakrzewski, J. Time crystals: a review. Rep. Progr. Phys. 81, 016401 (2017). 11. Sacha, K. Time Crystals, Springer Series on Atomic, Optical, and Plasma Physics (Springer Cham, 2020). 12. Zaletel, M. P. et al. Colloquium: Quantum and classical discrete time crystals. Rev. Mod. Phys. 95, 031001 (2023). 13. Autti, S., Eltsov, V. B. & Volovik, G. E. Observation of a time quasi- crystal and its transition to a superfluid time crystal. Phys. Rev. Lett. 120, 215301 (2018). 14. Kongkhambut, P. et al. Observation of a continuous time crystal. Science 377, 670 (2022). 15. Chen, Y.-H. & Zhang, X. Realization of an inherent time crystal in a dissipative many-body system. Nat. Commun. 14, 6161 (2023). 16. Wu, X. et al. Dissipative timecrystal in a strongly interacting rydberg gas. Nat. Phys. 20, 1389 (2024). 17. Liu, T., Ou, J.-Y., MacDonald, K. F. & Zheludev, N. I. Photonic metamaterial analogue of a continuous time crystal. Nat. Phys. 19, 986 (2023). 18. Carraro-Haddad, I. et al. Solid-state continuous time crystal in a polariton condensate with a built-in mechanical clock. Science 384, 995 (2024). 19. Ono, K. & Tarucha, S. Nuclear-spin-induced oscillatory current in spin-blockaded quantum dots. Phys. Rev. Lett. 92, 256803 (2004). 20. Yusa, G., Hashimoto, K., Muraki, K., Saku, T. & Hirayama, Y. Self- sustaining resistance oscillations: Electron-nuclear spin coupling in mesoscopic quantum hall devices. Phys. Rev. B 69, 161302 (2004). 21. Hennel, S. et al. Nonlocal polarization feedback in a fractional quantum hall ferromagnet. Phys. Rev. Lett. 116, 136804 (2016). 22. Zhang, J. et al. Observation of a discrete time crystal. Nature 543, 217 (2017). 23. Choi, S. et al. Observation of discrete time-crystalline order in a disordered dipolar many-body system. Nature 543, 221 (2017). 24. Randall, J. et al. Many-body-localized discrete time crystal with a programmable spin-based quantum simulator. Science 374, 1474 (2021). 25. Smits, J., Liao, L., Stoof, H. T. C. & van der Straten, P. Observation of a space-time crystal in a superfluid quantum gas. Phys. Rev. Lett. 121, 185301 (2018). 26. Kyprianidis, A. et al. Observation of a prethermal discrete time crystal. Science 372, 1192 (2021). 27. Keßler, H. et al. Observation of a dissipative time crystal. Phys. Rev. Lett. 127, 043602 (2021). 28. Kongkhambut, P. et al. Realization of a periodically driven open three-level dicke model. Phys. Rev. Lett. 127, 253601 (2021). 29. Zhu, B., Marino, J., Yao, N. Y., Lukin, M. D. & Demler, E. A. Dicke time crystals in driven-dissipative quantum many-body systems. N. J. Phys. 21, 073028 (2019). 30. Skulte, J. et al. Parametrically driven dissipative three-level dicke model. Phys. Rev. A 104, 063705 (2021). 31. Taheri, H., Matsko, A. B., Maleki, L. & Sacha, K. All-optical dissipative discrete time crystals. Nat. Commun. 13, 848 (2022). 32. Jensen, M. H., Bak, P. & Bohr, T. Transition to chaos by interaction of resonances in dissipative systems. I. Circle maps. Phys. Rev. A 30, 1960 (1984). 33. Pikovsky, A., Rosenblum,M. & Kurths, J. Synchronization of periodic oscillators by periodic external action. In: Synchronization: A Uni- versal Concept in Nonlinear Sciences, Cambridge Nonlinear Sci- ence Series p. 175–221. (Cambridge University Press, 2001). Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 7 www.nature.com/naturecommunications 34. Faraday, M. Xvii. on a peculiar class of acoustical figures; and on certain forms assumed by groupsof particles upon vibrating elastic surfaces. Phil. Trans. R. Soc. 121, 299–340 (1831). 35. Goldstein, R. E. Coffee stains, cell receptors, and time crystals: Lessons from the old literature. Physics Today 71, 32 (2018). 36. Landau, L. D. and Lifshitz, E. M. Mechanics, Third Edition: Volume 1 (Course of Theoretical Physics), 3rd ed. (Butterworth-Heine- mann, 1976). 37. Adler, R. A study of locking phenomena in oscillators. Proc. IRE 34, 351 (1946). 38. Strogatz, S. Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Westview Press, 2001). 39. John Argyris, M. H. R. F. Gunter Faust An exploration of dynamical systems and chaos (Springer, 2015) Chap. 5 Dynamical Systems with Dissipation, pp. 189–298. 40. Rosenstein,M. T., Collins, J. J. &DeLuca,C. J. Apracticalmethod for calculating largest Lyapunov exponents from small data sets. Phys. D Nonlinear Phenomena 65, 117 (1993). 41. Hegger, R., Kantz, H. & Schreiber, T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos 9, 413 (1999). 42. John Argyris, M. H. R. F. Gunter Faust An exploration of dynamical systems and chaos (Springer, 2015) pp. 473–592. 43. D’yakonov, M. I., Merkulov, I. A. & Perel’, V. I. Optical-orientation anisotropy produced in semiconductors by quadrupole splitting of the spin levels of the lattice nuclei. Sov. Phys. JETP 49, 160 (1979). 44. Fleisher, V. G. andMerkulov, I. A.Optical orientation (North-Holland, Amsterdam, 1984) Chap. 5. 45. D’yakonov, M. I., Merkulov, I. A. & Perel’, V. I. Instabilities in the spin system of optically oriented electrons and nuclei in semi- conductors. Sov. Phys. JETP 51, 175 (1980). 46. Jensen, M. H., Bak, P. & Bohr, T. Complete Devil’s staircase, fractal dimension, and universality of mode- locking structure in the circle map. Phys. Rev. Lett. 50, 1637 (1983). 47. Bak, P. The devil’s staircase. Physics Today 39, 38 (1986). 48. Wu, X. et al. Farey tree and devil’s staircase of frequency-locked breathers in ultrafast lasers. Nat. Commun. 13, 5784 (2022). 49. Ramos-Pérez, I. A. et al. Theory of optomechanical locking in driven- dissipative coupled polariton condensates. Phys. Rev. B 109, 165305 (2024). 50. Arnold, V. Small denominators i, mappings of the circumference onto itself. Trans. Am. Math. Soc. Ser. 2 46, 213–284 (1965). 51. Bak, P., Bohr, T. & Jensen, M. H. Mode-locking and the transition to chaos in dissipative systems. Phys. Scripta 1985, 50 (1985). 52. Brown, S. E., Mozurkewich, G. & Grüner, G. Subharmonic shapiro steps and devil’s-staircase behavior in driven charge-density-wave systems. Phys. Rev. Lett. 52, 2277 (1984). 53. Baums, D., Elsässer, W. & Göbel, E. O. Farey tree and devil’s stair- case of a modulated external-cavity semiconductor laser. Phys. Rev. Lett. 63, 155 (1989). 54. Matus, P. & Sacha, K. Fractional time crystals. Phys. Rev. A 99, 033626 (2019). 55. Pizzi, A., Knolle, J. &Nunnenkamp, A. Period-ndiscrete time crystals and quasicrystals with ultracold bosons. Phys. Rev. Lett. 123, 150601 (2019). 56. Pizzi, A., Knolle, J. & Nunnenkamp, A. Higher-order and fractional discrete time crystals in clean long-range interacting systems. Nat. Commun. 12, 2341 (2021). 57. Jiao, Y. et al. Observation of a time crystal comb in a driven- dissipative system with rydberg gas. arXiv, https://arxiv.org/abs/ 2402.13112 (2024). 58. Liu, B. et al. Bifurcation of time crystals in driven and dissipative rydberg atomic gas. arXiv, https://arxiv.org/abs/2402. 13644 (2024). 59. Liu, B. et al. Higher-order and fractional discrete time crystals in floquet-driven rydberg atoms. arXiv, https://arxiv.org/abs/2402. 13657 (2024). 60. Kongkhambut, P. et al. Observation of a phase transition from a continuous to a discrete time crystal. Rep. Progr. Phys. 87, 080502 (2024). 61. He, G. et al. Experimental realization of discrete time quasicrystals. Phys. Rev. X 15, 011055 (2025). 62. Kuramoto, Y. & Battogtokh, D. Coexistence of coherence and incoherence in nonlocally coupled phase oscillators. arXiv, https:// arxiv.org/abs/cond-mat/0210694 (2002). 63. Abrams, D. M. & Strogatz, S. H. Chimera states for coupled oscil- lators. Phys. Rev. Lett. 93, 174102 (2004). 64. Panaggio, M. J. & Abrams, D. M. Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators. Nonlinearity 28, R67 (2015). 65. Clerc, M. G. & Verschueren, N. Quasiperiodicity route to spatio- temporal chaos in one-dimensional pattern-forming systems. Phys. Rev. E 88, 052916 (2013). 66. Russomanno, A. Spatiotemporally ordered patterns in a chain of coupled dissipative kicked rotors. Phys. Rev. B 108, 094305 (2023). 67. Hentschel, H. & Procaccia, I. The infinite number of generalized dimensions of fractals and strange attractors. Phys. D Nonlinear Phenomena 8, 435 (1983). 68. Glazier, J. & Libchaber, A. Quasi-periodicity anddynamical systems: An experimentalist’s view. IEEE Trans. Circ. Syst. 35, 790 (1988). Acknowledgements The authors are thankful to D. R. Yakovlev for fruitful discussions. A.G. and M.B. acknowledge support by the BMBF project QR.X (Contract No.16KISQ011). The Resource Center “Nanophotonics" of Saint- Petersburg State University provided the epilayer sample. Author contributions A.G. and N.E.K. contributed equally to this paper. A.G. built the experi- mental apparatus and performed themeasurements. P.A.H. contributed to the fine-step measurements of the contour map in Fig. 5a. N.E.K. and A.G. analyzed the data. N.E.K. and V.L.K. provided the theoretical description. All authors contributed to the interpretation of the data. N.E.K., V.L.K., and A.G. wrote the manuscript in close consultation with M.B. Funding Open Access funding enabled and organized by Projekt DEAL. Competing interests The authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-58400-6. Correspondence and requests for materials should be addressed to Alex Greilich or Nataliia E. Kopteva. Peer review information Nature Communications thanks Marcel G. Clerc, Leon Zaporski, and the other, anonymous, reviewer for their con- tribution to the peer review of this work. A peer review file is available. Reprints and permissions information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jur- isdictional claims in published maps and institutional affiliations. Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 8 https://arxiv.org/abs/2402.13112 https://arxiv.org/abs/2402.13112 https://arxiv.org/abs/2402.13644 https://arxiv.org/abs/2402.13644 https://arxiv.org/abs/2402.13657 https://arxiv.org/abs/2402.13657 https://arxiv.org/abs/cond-mat/0210694 https://arxiv.org/abs/cond-mat/0210694 https://doi.org/10.1038/s41467-025-58400-6 http://www.nature.com/reprints www.nature.com/naturecommunications Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2025 Article https://doi.org/10.1038/s41467-025-58400-6 Nature Communications | (2025) 16:2936 9 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ www.nature.com/naturecommunications Exploring nonlinear dynamics in periodically driven time crystal from synchronization to chaotic motion Results Periodic auto-oscillations Periodic modulation of excitation polarization Model of periodically modulated ENSS Discussion Methods Setup Details on data processing Details on simulation model Fractal dimension Data availability Code availability References Acknowledgements Author contributions Funding Competing interests Additional information