Long-lived coherence in driven spin systems: from two to infinite spatial dimensions

Long-lived coherences, emerging under periodic pulse driving in the disordered ensembles of strongly interacting spins, offer immense advantages for future quantum technologies, but the physical origin and the key properties of this phenomenon remain poorly understood. We theoretically investigate this effect in ensembles of different dimensionality, and predict existence of the long-lived coherences in two-dimensional and infinite-dimensional (where every spin is coupled to all others) systems, which are of particular importance for quantum sensing and quantum information processing. We explore the transition from two to infinite dimensions, and show that the long-time coherence dynamics in all dimensionalities is qualitatively similar, although the short-time behavior is drastically different, exhibiting dimensionality-dependent singularity. Our study establishes the common physical origin of the long-lived coherences in different dimensionalities, and suggests that this effect is a generic feature of the strongly coupled spin systems with positional disorder. Our results lay out foundation for utilizing the long-lived coherences in a range of application, from quantum sensing with two-dimensional spin ensembles, to quantum information processing with the infinitely-dimensional spin systems in the cavity-QED settings.


I. INTRODUCTION
Collective quantum coherences of many-spins systems play central role in quantum science and technology. But quantum coherence is fragile, and extending its lifetime is a critical problem. For instance, in spin ensembles the collective coherence (collective transverse polarization) is destroyed by dipolar interactions [1][2][3][4], and the the spin echo signal, which quantifies coherence, quickly decays on the time scale T 2 [5]. Coherence can be preserved e.g. via pulse and/or continuous-wave decoupling that suppresses dipolar interactions [1][2][3]. Recently, an intriguing alternative has attracted much attention: it exploits, rather than fights, the spin-spin interactions. Namely, the unusual many-spin states, which are formed in ensembles of dipolar-coupled spins under periodic driving by π-pulses, exhibit collective spin coherences living up to 10 5 times longer than the T * 2 and about 10 4 time longer than the T 2 time [6][7][8][9]. This phenomenon, along with other similar effects, has been observed in various solid-state nuclear magnetic resonance (NMR) experiments [6,7,[9][10][11][12][13], but still remains poorly understood. The long-lived coherences emerge from the combination of strong dipolar coupling, disorder, and pulse imperfections (understood broadly as deviations of the real π-pulses from perfect instant 180 • rotations) [6,12]. Besides, the long-lived coherences demonstrate subharmonic response, i.e. asymmetry in the magnitudes of even and odd echoes [6,9,12]. Both effects of long coherence lifetime and its subharmonic response are remarkably stable against perturbations and decoherence. a walter.hahn@uibk.ac.at b v.v.dobrovitski@tudelft.nl; author to whom correspondence should be addressed The long-lived coherences could be of great benefit for new quantum technology platforms. For instance, promising platforms for quantum sensing utilize twodimensional systems (d = 2, see Fig. 1), such as surface spins or 2D layers of NV spins [14][15][16][17]), which can be brought close to the system being sensed, thus improving resolution and sensitivity. Employing ensembles of spins boosts the total signal, and thus greatly improves the signal-to-noise ratio, but the collective coherence decays quickly due to dipolar coupling between the spins. Increasing the lifetime of collective coherences would be of enormous benefit for quantum sensing. On the opposite end (d → ∞) are the spin ensembles in a cavity QED-type settings, actively explored for quantum information applications [18,19], where each spin is coupled to all others with a similar strength via collective photonic, phononic, or magnonic mode [20][21][22][23][24][25][26]. Taking full advantage of the long-lived coherences could increase the signal-to-noise ratio in these systems by orders of magnitude. In order to achieve that, detailed understanding the long-lived coherences in systems of different spatial dimensionality d is required. So far, even existence of the long-lived coherences at d = 2 or d → ∞ has remained elusive, and their properties have been unknown.
In this article we predict that the long-lived coherences do exist in these important systems, thus opening the way to employing them in novel quantum information platforms. In order to analyze in detail the transition from d = 2 to d → ∞, and to clarify generic features of the long-lived coherence dynamics, we numerically simulate the dynamics of disordered dipolar-coupled quantum spin ensembles of different dimensionalities, subject to periodic driving by imperfect π-pulses, with the rotation angle slightly deviating from 180 • . For all dimensionalities d studied, we observe the long-lived coherences and see emerging subharmonic response when the time inter- two-dimensional, three-dimensional, and (effectively) infinitedimensional systems. The latter system with all-to-all interactions of similar strength is realized e.g. by coupling the spins to a detuned collective photon/phonon/magnon mode in a cQED-like setting. (b) Schematic representation of the periodic pulse sequence. Imperfect π pulses (P) are applied at times (2n + 1)τ (n ≥ 0), so that a series of echoes (E) is formed at times 2nτ . Each pulse rotates the spins along the x-axis by an angle π(1 + ), where is the pulse imperfection parameter.
val τ between the pulses increases to become comparable to T 2 , see Fig. 1b. Our simulations show that the magnitude of the long-lived coherence decreases at larger d, but still remains quite large even at d → ∞. By analyzing the Floquet operator, we establish the kinematic origin of the long-lived coherences, determine that it is similar for all dimensionalities, and identify the states that are involved in their formation.
Some aspects of the long-lived coherence resemble the time crystal dynamics in periodically driven spin systems [27][28][29][30][31][32][33][34][35], with their characteristic robustness [36][37][38]. Time crystals have been observed in many spin systems, from trapped ions [39] to spin ensembles in diamond [40], including the kind of NMR systems that exhibits the long-lived coherences [35,41,42]. However, the physical origins of the two phenomena are different: the coherences are determined by the transverse magnetization, while the time crystal dynamics refers to the behavior of the longitudinal polarization. They are governed by different processes, and their lifetimes can differ by many orders of magnitude [1,2]; e.g., in the case of perfect pulses (instantaneous 180 • rotations) the time crystal oscillations have the longest lifetime and largest amplitude, while the long-lived coherences completely van-ish [6,9,12]. The relation between the long-lived coherence and the time crystal dynamics is poorly explored, in spite of its fundamental interest. Besides, the studies of time crystals so far have mostly focused on d = 1 [35,37,38,[43][44][45][46][47][48][49][50] and d = 3 systems [39][40][41][42]51], as well as on the systems with d → ∞ [21,36,40,[50][51][52]. For these reasons, we also explore here the time crystal dynamics of the longitudinal spin polarization for different dimensionalities. We find that it also demonstrates time crystal-like oscillations with a very long lifetime for all d, but the physical features are markedly different from those of the long-lived coherences. More detailed exploration of this issue in the future would be of utmost interest, but is beyond the scope of the present paper.
Dimensionality plays a key role in the dynamics of dipolar-coupled positionally disordered spin systems. Such systems exhibit characteristic d-dependent dynamical singularities: the spin dynamics at short times is strongly singular for d = 2, while for d → ∞ the singularity disappears [2,[53][54][55][56]. Besides, dimensionality is well known to be decisive in the context of localization and thermalization dynamics in spin ensembles [28,38,[57][58][59]. Since so many key features of the spin dynamics depend on d, one would expect that the properties of long-lived coherences would also strongly depend on d. Surprisingly, our results demostrate that this expectation is incorrect.

II. QUALITATIVE DISCUSSION OF THE EFFECT
We study an ensemble of N s spins S i = 1/2 (i = 1 . . . N s ) in a standard setting of a magnetic resonancetype experiments. Namely, the spin system is placed in a strong quantizing magnetic field H Q directed along the z-axis [1]; this field induces fast spin precession with Larmor frequency ω Q , which is much larger than all other frequency scales of the problem [60]. Following the standard theory of magnetic resonance, we describe spin dynamics in the coordinate frame that rotates around the z-axis with the circular frequency ω Q , and retain only the secular terms in the system's Hamiltonian, which remain static in the rotating frame [61], or vary slowly in comparison with ω Q [1,2].
Initially, by applying a preparatory π/2 pulse, the spins are prepared in a state weakly polarized along the x-axis of the rotating frame, such that the initial ensemble density matrix is ρ(t = 0) ∝ I − µM x , where I is identity matrix, µ 1 is a parameter determining the absolute polarization of the ensemble, and M x = j S jx is the collective coherence operator. Here and below, S jα with α = {x, y, z} denotes the component of the j-th spin along the rotating-frame axis α. In experiments, the ensemble coherences along the x-and y-directions are quantified by the total transverse magnetizations M x and M y along the corresponding axes, so we use the terms "coherence" and "transverse magnetization" interchangeably. The longitudinal polarization, exhibiting time crystallike behavior (see Sec. V), is quantified by the magnetization M z = j S jz along the z-axis. Note that in this work we vary only the spatial dimensionality d of the ensemble, while spins themselves remain embedded in three dimensions, i.e. have three orthogonal components.
In typical experiments, the spins experience random quasi-static local magnetic fields, described by the Hamiltonian H L = Ns j=1 h j S jz ; everywhere in this article we set = 1 and normalize the spins' gyromagnetic ratio to γ = 1.
[62] These fields cause fast dephasing: the x-component of each spin S jx oscillates at its own rate proportional to h j , and the collective coherence M x (t) = Tr[ρ(t)M x ] vanishes at the timescale T * 2 . This is usually too short for practical needs, and dephasing is suppressed by applying a number of hard π pulses, which reverse the sign of the Hamiltonian H L (ideal, i.e. instantaneous 180 • hard pulse along the x-axis performs rotation S iz → −S iz , S iy → −S iy for all spins at once). A pulse applied at t = τ restores collective coherence, producing Hahn echo signal at t = 2τ [1,2] (see Fig. 1b).
However, the hard π pulses do not affect the spinspin interaction that destroys the ensemble coherence by entangling different spins [1][2][3][4]. In relevant experiments, the dominant interaction is the dipolar coupling, described by the Hamiltonian [1,2] ij is the coupling constant between the spins i and j, r ij = | r ij | is the distance between them, and θ ij is the polar angle of r ij . Under the influence of H I , the Hahn echo signal gradually decays as a function of time t = 2τ at the timescale t ∼ T 2 T * 2 . Note that the dipolar interaction is long-ranged, so that each spin is coupled to all other spins for all systems considered here, for all values of d from 2 to ∞, and the coupling strength J ij decays with distance r ij between the spins i and j in the same way for all d. However, the statistical properties of the values J ij greatly vary with the system's dimensionality, thus leading to dramatically different decay of the Hahn echo signal.
The rate of the Hahn echo decay is governed by the positional disorder of spins. In relevant experiments [9-11, 14-19, 40], the spins are very dilute: they occupy a small fraction of the lattice sites, and are randomly distributed in the sample. As a result, each spin S i has its own set of dipolar couplings J ij to the other spins, such that different spins "feel" different environments made of other spins. In low-dimensional ensembles these spin-tospin variations are very strong [53][54][55], leading to a very fast decay of the collective coherence. In Appendix A we show that, in the limit of T 2 T * 2 , when the flip-flop terms can be omitted and the Hamiltonian (1) acquires the form J ij S iz S jz , the Hahn echo signal M x (t) decays with time as for d ≤ 5, such that for d = 2 the initial decay is infinitely fast. For larger d the fluctuations are not as strong, and the singularity of M x (t) at t = 0 is weak for d = 4 and 5. In the limit d → ∞ each spin is coupled to all others with almost uniform coupling, so the echo decay acquires Gaussian form M x (t) ∝ exp − (t/T 2 ) 2 , without any singularity at all. The total Hamiltonian of the system, taking into account the pulse driving, is where H P (t) describes periodic driving by a train of (generally imperfect) π-pulses, as shown in Fig. 1b. If the pulses were ideal, they would suppress the dephasing term H L (see Appendix A for details), and leave the dipolar interaction H I intact. Thus, the echo signal M x (t) would not depend on the number of pulses applied during the time t, and the echo decay would follow Eq. 2.
Our results show that for non-ideal pulses this is true at short times: the dynamics drastically depends on d, with the initial decay rate varying from infinity at d = 2 to zero at d → ∞. At long times, for spin ensembles of all dimensionalities d, the spin dynamics is controlled by accumulation of the pulse imperfections. This process depends on the specific pulse sequence [6,7,9]; here we focus on the simple and efficient Carr-Purcell-Meiboom-Gill (CPMG) protocol shown in Fig. 1b. Accumulation of the pulse errors, combined with the dipolar spin-spin coupling, leads to long-lived tails in the echo signal, extending far past T 2 time, in the region where the Hahn echo has already vanished. The magnitude of the longlived coherence tails is particularly large when the interpulse time delay τ is short compared to T 2 , as seen in Fig. 2a for d = 2. The effect is remarkably robust, and does not disappear even when τ becomes comparable to T 2 . In that regime, another feature emerges for all values of d: the long-lived coherences exhibit pronounced subharmonic response (Fig. 2b), with even echoes being larger than odd ones. This behavior reflects the fact that the period of the evolution operator for CPMG protocol is twice the period of the Hamiltonian H P (t) of the CPMG pulse train [3] (see Sec. III).
We use direct numerical solution of the many-spin time-dependent Schrödinger equation with the secondorder Suzuki-Trotter [63] or Chebyshev [64] expansion of the evolution operator U (t) for up to N s = 24 spins.
In all situations that we tested, both methods give consistent results. The initial mixed-state density matrix is represented by a random pure state, i.e. as a random unit-norm complex-valued vector of length 2 Ns , uniformly sampled from a sphere S 2 Ns −1 of unit radius. We calculate the experimentally measured normalized magnetization response for initial polarization along the axis α. The spins are randomly placed in a d-dimensional cube, and averaging is performed over 100-360 realizations of the spatial arrangements and values of local fields. The cube's edge length (i.e. the spin density f s ) is adjusted to have T 2 = 1 for each system after averaging (see Appendix A for relation between T 2 and f s ). T 2 is defined as the time when the echo magnitude decreases by the factor e. In experiments the pulse imperfections may arise from deliberate or accidental miscalibration, or also due to local fields and dipolar interactions, which affect the spin rotation during the finite-duration pulses [6,9,12]. For hard short pulses, imperfections are small, and in this article we model nonideal pulses as instantaneous rotations around the x-axis by an angle π(1 + ) with 1.

III. DYNAMICS OF THE LONG-LIVED COHERENCES
The behavior of M x (t) in pulse driven spin ensembles, as described above, is shown in Fig. 2 for d = 2, and in Fig. 3  We consider two different values of the inter-pulse delay τ : short τ ≈ 0.07 T 2 , and long τ ≈ 0.7 T 2 (recall that T 2 = 1). At short times t T 2 , the system's response M x (t) closely follows the Hahn echo, and is in excellent agreement with our analytical predictions (see Appendix A for details), as seen in the figures (the small differences between the simulated Hahn echo decay and the analytical predictions are due to the finite-size effects). At longer times t T 2 , the magnetization response M x (t) exhibits long-time tails for all dimension-alities d considered, and for both short and long τ . These long-time tails extend far beyond the T 2 time, and in experiments are likely to be limited by the spin-lattice relaxation time T 1 . With increasing d, the overall amplitude of the long-time echoes becomes somewhat smaller. Still, even in the limit d → ∞, the long-time tails do not vanish but saturate at a nonzero value. Also, the amplitude of the long-time tails becomes smaller as the inter-pulse delay τ increases.
When the inter-pulse delay τ is large, comparable to T 2 , the long-lived coherences demonstrate pronounced subharmonic dynamics (Fig. 2b), where the period of the magnetization response 4τ is twice longer than the period of driving 2τ . The subharmonic response was observed for all values of d we modeled. This feature can be rationalized with the notion of the cycling period of a pulse sequence [3]: for the control Hamiltonian H P (t), the cycling period t c is defined by the condition of periodicity of the evolution operator, i.e.
It is known that the cycling period of a pulse sequence can differ from the period of the Hamiltonian [3], e.g. for ideal π-pulses, the cycling period of the CPMG protocol is t c = 4τ , and includes two π-pulses [67], i.e. contains two periods of the underlying control Hamiltonian H P (t). We also note that the evenodd asymmetry is present for short τ , but its amplitude is too small to be easily seen.
It is important to point out that the long-lived coherences and the subharmonic response crucially depend on the driving sequence. For instance, for an alternating-phase Carr-Purcell (APCP) sequence, where the direction of driving alternates between the positive and the negative x-direction, no long-lived tails of M x (t) emerge [6].

IV. FLOQUET OPERATOR ANALYSIS OF THE LONG-LIVED COHERENCES
Periodically driven systems, described by a Hamiltonian obeying H(t + 2τ ) = H(t), can be analyzed using Floquet theory. The stroboscopic time evolution of the system's density matrix, considered only at times that are integer multiples m of the driving period 2τ , can formally be written as is the operator of rotation around the x-axis by an angle π(1 + ), and U H (τ ) is the operator of evolution under the action of the system's internal Hamiltonian H I + H L .
As a unitary operator, U F possesses a complete orthonormal set of eigenstates |ψ k , with complex eigenvalues e iφ k of unit modulus:  For short τ (panel a), a large number of large entries concentrate on the diagonal φj ≈ φ k , producing long-lived coherence with noticeable amplitude. For long τ (panel b), the values at semi-diagonals |φj − φ k | ≈ π become comparable to the values on the diagonal, which leads to noticeable subharmonic response. Overall, the values for long τ are smaller than those for short τ . The system size is Ns = 14. Only values larger than 0.25 are included in the figures.
so the magnetization response after m driving periods is The signal M x (2τ m) is therefore mainly determined by two quantities: firstly, by the magnitude of the matrix elements M jk x = | ψ j |M x |ψ k | 2 , and, secondly, by the distribution P (φ j −φ k ) of the quasienergy differences φ j − φ k , i.e. by the number of Floquet eigenstates |ψ j and |ψ k with a given difference in the quasienergies φ j and φ k . The terms with j = k in Eq. (6) do not depend on m, so the long-time response M x (2τ m) is governed  Fig. 4 for short τ (a) and long τ (b) for a two-dimensional system, for one typical realization of the positional disorder. The matrix M jk x for short τ is dominated by large diagonal elements, whereas the off-diagonal elements are almost negligible. This corresponds to the pronounced long-lived coherences in Fig. 2a and Fig. 3, with the almost time-independent amplitude. It is clearly seen that the long-lived coherence contains comparable contributions from a large number of Floquet eigenstates, rather than being confined to a small subset of some special states.  long τ (b). The sharp peaks at ∆φ = 0 and at ∆φ = π correspond to the long-lived coherences and to the subharmonic response, respectively. Specifically, P (∆φ) is the total number of the Floquet eigenstates |ψj and |ψ k having a given difference ∆φ in their quasienergies φj and φ k . The quantity ∆φ is downfolded to the interval ∆φ ∈ [0, π] and binned, so that P (∆φ) includes all states whose quasienergies satisfy the condition |φj − φ k | ∈ [∆φ − β, ∆φ + β] or 2π − |φj − φ k | ∈ [∆φ − β, ∆φ + β], where 2β = π · 10 −5 is the width of a bin. To avoid double counting, only the states with φj ≥ φ k are included in P (∆φ). The results shown are averaged over many realizations of the disorder. The inset in panel (b) shows the magnified view of the peak at ∆φ ≈ π. The simulation parameters are the same as in Fig. 2.
For long τ , the matrix M jk x exhibits large entries both on the diagonal, and on the two semi-diagonal lines corresponding to φ j − φ k ≈ ±π; the latter correspond to emerging even-odd echo asymmetry seen in Fig. 2b. The semi-diagonals also show comparable contributions from a large number of pairs of Floquet eigenstates. In comparison with the case of short τ , the diagonal values M jj x on average are smaller for long τ , corresponding to smaller amplitude of the long-time tails in Fig. 2b.
The matrices M jk x for other d exhibit similar structure. As an example, Fig. 5 exhibits the results for M jk x in the case of d = 5. The results of diagonalization of the Floquet operator for different d evidence that the physical origin of the effect of the long-lived coherences is similar for all dimensionalities.
The other quantity determining the signal M x (2τ m) in Eq. (6) is the distribution P (∆φ) of the quasienergy differences |φ j − φ k |. In order to take into account the symmetries of the summands in Eq. 6, we downfold the quantity ∆φ to the interval [0, π], i.e. we take ∆φ = |φ j − φ k | when |φ j − φ k | ≤ π, and ∆φ = 2π − |φ j − φ k | when |φ j − φ k | > π (so that ∆φ is the smallest angular distance between φ j and φ k on the S 1 circle). The binned distribution P (∆φ) for d = 2 is shown in Fig. 6, the bin width is 2β = π ·10 −5 . A peak at ∆φ ≈ π clearly emerges for long τ . This means that the number of quasienergy pairs with |φ j −φ k | ≈ π becomes larger for long τ . Hence, the subharmonic response emerges not only because the values of M jk x on the semi-diagonals become larger, but also because the total number of non-zero entries on the semi-diagonals increases.

V. LONGITUDINAL MAGNETIZATION AND THE INFINITE-TEMPERATURE TIME CRYSTAL DYNAMICS
Let us now focus on the dynamics of the longitudinal polarization M z (t), which under some circumstances can exhibit robust long-lived infinite-temperature time crystal-like dynamics [27].
Long-lived coherences and time crystal dynamics share a number of similarities: both are induced by strong spinspin interactions, robust to experimental imperfections, and demonstrate subharmonic dynamics under appropriate conditions. However, the underlying physics is totally different. Since the system's internal Hamiltonian H I + H L conserves the total z-magnetization, the signal M z (t) remains constant between the pulses. If the πpulses were ideal, then M z (t) would just switch between +1 and −1. In contrast, the coherence M x (t) would quickly vanish under the action of the dipolar spin-spin coupling, along with the Hahn echo, at the timescale T 2 , without any long-lived tail.
For non-ideal pulses, one would expect the longitudinal polarization to decay with increasing the number of applied pulses. Indeed, for > 0, after an imperfect pulse the absolute value of M z (t) would decrease by a factor cos(π ), while the y-component increases by sin(π ). Accordingly, if the y-component has irreversibly dephased to zero during the time 2τ between pulses, then the absolute value of M z (t) would be expected to decay monotonically as cos m (π ) with increasing the number m of applied pulses [35,41].
The actual behavior of the longitudinal magnetization response M z (t) in a d = 2 disordered spin system with N s = 20 is shown in Fig. 7. Initial decay roughly follows the expected cos m (π ) pattern for both short and long τ (panels (a) and (b), respectively); the additional modulation seen in panel (a) is likely due to incomplete dephasing of the y-component between the pulses due to short τ .
At later times, however, the absolute magnitude of Each π-pulse flips the z-magnetization, so that the subharmonic response with the period 4τ is the dominating component of the system's long-time response. Without dipolar coupling, accumulated pulse error would modulate the magnetization along the z-axis, and Mz(t) would decay as cos m (π ) after m pulses (dashed black line). However, in the presence of the dipolar coupling, Mz(t) exhibits long-time tails, alternating between the directions "up" and "down" after each π-pulse. The system size is Ns = 20, = 0.07, and T * 2 ≈ 0.02.
driving sequence. This long-lived subharmonic response is seen for all dimensionalities we studied: the results for another example, a d = 8 spin ensemble, are presented in Fig. 8, and are very similar, except that the amplitude of the long-time tail is somewhat smaller than in the d = 2 case. This conclusion is also consistent with the experimental evidence reported for d = 3 spin systems [40][41][42]. Note that the initial decay of M z (t) does not exhibit the singularities present in the short-time dynamics of M x (t), because it is governed by a different physical process, by accumulation of the rotation errors. Likewise, the long-time behavior of the longitudinal M z (t) and the transversal M x (t) polarization is also qualitatively different. For instance, with increasing the inter-pulse We have also analyzed the stroboscopic evolution of M z (t) by diagonalizing the Floquet operator; an example of the corresponding matrix M jk z = | ψ j |M z |ψ k | 2 is shown in Fig. 9 for one typical realization of the positional disorder. As anticipated, it exhibits nonzero values only at the semi-diagonals φ j − φ k ≈ ±π, in accordance with the subharmonic response described above. The diagonal elements are negligible. Similar to the case of long-lived coherences, we see again that many states, distributed all over the Hilbert space, contribute to the effect.
For the initial polarization along the y-axis, no longlived magnetization response M y (t) is seen in any dimensionality d, and the elements of the matrix M jk y = | ψ j |M y |ψ k | 2 are also generally small, without any clear structure.

VI. FINAL REMARKS AND CONCLUSIONS
Before concluding, let us make some final remarks. a) Time crystal dynamics and the long-lived coherences exhibit striking similarity: both are induced by spin-spin interactions and demonstrate subharmonic dynamics under appropriate conditions. At the same time, these effects arise in very different regimes, demonstrate very different dynamics, and are differently affected by the pulse imperfections. Understanding the connection between time crystals (in particular, infinite-temperature time crystals [35]) and long-lived coherences is an interesting and important, yet rather unexplored problem.
b) The presence of long-time tails along the z-axis which is perpendicular to the pulse driving field implies that the fundamental process for establishing long-time tails of the magnetization response may not be spin locking [1,2] as may appear in analogy with other similar effects [10,11,68]. c) Our simulations show that even if the direction of the driving during the π-pulses is chosen randomly for each spin, the long-lived coherences emerge and persist for long times, as long as the direction of the driving remains constant in time. This observation suggests that the origin of the long-lived coherences is primarily kinematic; it may arguably also involve many-body localization [28,31], transient prethermal regime [32,69,70], or spin-glass like behaviour.
d) The results shown in the present article were obtained for a fixed pulse imperfection = 0.07. Our numerical simulations with other values of demonstrate that qualitatively the same behavior occurs for other choices of 1. The long-time spin coherences persist even when is chosen randomly for different spins, as long as it remains constant in time. e) We considered in this article spin systems of finite size, with the total number of spins N s ≈ 20. Our simulations demonstrate very modest quantitative changes as the system size has been varied from N s = 10 to N s = 24. Besides, our numerical results for d = 3 agree with the reported experimental results for three-dimensional systems [6,9,40,41]. Summarizing, we have investigated disordered dipolarcoupled quantum many-spin systems of different spatial dimensionality subjected to periodic driving (CPMG protocol). Depending on the dimensionality, such systems exhibit singularities in short-time spin dynamics, which are caused by statistical fluctuations in the dipolar spinspin interaction strength. For all dimensionalities, we observed the long-lived spin polarization along the driving pulses, and along the axis conserved by the internal Hamiltonian (z-axis). The amplitude of the long-lived magnetization depends on the inter-pulse time delay and on dimensionality. The Floquet operator analysis shows that the long-lived coherences M x (t) contain comparable contributions from a large number of Floquet eigenstates. Our results imply that the long-lived coherences and subharmonic response are generic features of dipolar-coupled disordered spin systems, including two-, three-, and infinite-dimensional systems that are particularly relevant for practical applications. Developing specific protocols for such applications is an exciting avenue for further research.
[58] L.  [40,41,51] the role of strong quantizing field is played by strong continuouswave Rabi driving applied at the frequency ωQ. In this case, the dynamics is to be considered in a doubly rotating frame, where the effective quantization axis is directed along the Rabi driving in the primary rotating frame, and the role of the principal Larmor frequency ωQ is played by the Rabi frequency ωR, see e.g. Refs. 1 and 2 for details. The doubly rotating frame in the theory of magnetic resonance is analogous to the dressed state basis in the quantum optics context. [65] For d = 2 the spins are located on a plane, and the zaxis (direction of the quantizing field) can be aligned at different angles with respect to this plane. If the z-axis is normal to the plane, all spin pairs have θij = π/2, while for the in-plane alignment this angle varies between 0 and π. However, the resulting changes in the Hahn echo are very small: for a fixed areal spin density fs, the values of T2 differ by only 4%. The long-lived coherences also do not exhibit any qualitative differences between these two cases. Fig. 2 so that the integral in Eq. (A2) is well defined at large r, where cos b(θ)t/r 3 → 1.
Singularity of the spin dynamics is determined by the competition of two effects: the dipolar coupling constants J ij decrease with increasing r, but the number of spins S j which interact with the given spin S i grows with r as the volume element dv = r d−1 drdA (where dA is the surface element of the (d − 1)-dimensional hypersphere of unit radius). For d < 6, this growth is sufficiently slow, so the integration over r can be extended to infinity, yielding the above-mentioned result (see Eq. (2) of the main text) . The integral over the hypersphere is a numerical factor of order of one, and the quantity (A6) which comes from integration over r, is also of order one, such that T 2 is determined just by the spin density f s . By renormalizing the spin density, we can set T 2 = 1 as explained above.
For d = 2, the singularity of the Hahn echo M H x (t) = exp − |t/T 2 | 2/3 is strong: the initial echo decay is infinitely fast due to strong fluctuations in the positions of the spins at small distances r ij . Although the typical distance between spins is of the order of 1/f a large fraction of spins has many neighbors at much smaller distances. Correspondingly, while the typical dipolar coupling is of the order of one (recall that we normalize f s to yield T 2 = 1), many realizations of the disorder produce very large dipolar couplings J ij .
Of course, in real crystals, at extremely short times the initial decay rate is finite, because the distance between spins is limited by the crystal lattice constant, which in turn limits the maximal dipolar interaction strength. However, the corresponding times are extremely small, orders of magnitude smaller than T * 2 , and are irrelevant for the phenomena considered here.
For d = 3, such fluctuations in J ij are less strong, and the singularity is weaker: M H x (t) = exp (− |t/T 2 |) has a cusp at t = 0, but the initial rate of decay is finite. Still, for both d = 2 and 3, the total spectral power of the resonance line is (formally) infinite. For d = 4 and 5, the fluctuations are less pronounced, the total spectral power of the resonance is finite, and the singularity in M H x (t) is weak.
For d ≥ 6, the integration over r cannot be extended to infinity: the integral Eq. (A6) diverges at small z (which correspond to r → ∞). This divergence means that the number of the spins S j coupled to the given spin S i grows too fast with increasing r. The contribution from the surface of the sample becomes important, dominating at larger d. The form of M H x (t) then depends on the sample shape and size, but in the limit d → ∞ it again acquires a universal shape-independent Gaussian form. For any regular-shaped sample in the limit d → ∞ all spins are located near the surface, and each spin pair is separated by almost the same distance, producing all-to-all interactions with a uniform coupling constant J ij → J. Eq. (A2) also demonstrates almost identical short-time behavior and long-lived coherences for the full and for the reduced models, as shown in Fig. 11a. However, in the case d → ∞ shown in Fig. 11b, the calculated signals for both full and reduced models coincide only at short times, and exhibit quantitative difference later. Although both models clearly demonstrate long-lived coherences, the echo amplitude is approximately twice higher for the reduced model as compared with the full model. A possible reason for this discrepancy is in the geometry of the dipolar couplings. The flip-flop process between two spins is important if the dipolar coupling between the two spins is comparable to, or larger than, the difference in the respective local fields. In the situation considered here, with T 2 T * 2 , the above case can only occur if the local fields of two spins are accidentally similar. In d = 2 systems, these two spins must be close to each other to ensure non-negligible dipolar coupling. In contrast, in d → ∞ systems, these two spins can be at any distance to each other because the dipolar interaction is homogeneous, coupling each spin to all other spins. Thus, the probability for two spins to have accidentally similar local fields, and to undergo a flip-flop process, is much larger in d → ∞ systems than in d = 2 systems.
In this article we use the full model in our simulations, with the exception of the results for the single-pulse Hahn echo shown in Figs. 2 and 3 of the main text. cient equal to | ψ j |M x |ψ k | 2 ; this corresponds to summation along diagonals of Figs. 4 and 5. Two δ-functions in the formula reflect the fact that the quantity ∆φ is downfolded to the interval [0, π] , i.e. we take ∆φ = |φ j − φ k | when |φ j − φ k | ≤ π and ∆φ = 2π − |φ j − φ k | when |φ j − φ k | > π; this downfolding reflects the symmetries of the summands in Eq. 6. Note that the quantity Σ(∆φ) is proportional to the cosine Fourier transform of the CPMG signal M x (2τ m).
To avoid double counting, only the states with φj ≥ φ k are included in Σ(∆φ). The normalization is chosen such that Σ(∆φ)d(∆φ) = 1. All other simulation parameters are the same as in Fig. 2.