⇦Chinese Physics Letters, 2020, Vol. 37, No. 5, Article code 050503Express Letter Imaginary Time Crystal of Thermal Quantum Matter * Zi Cai (蔡子)1,2**, Yizhen Huang (黄易珍)1, W. Vincent Liu (刘文胜)3,4,2** Affiliations 1Wilczek Quantum Center and Key Laboratory of Artificial Structures and Quantum Control, School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240 2Shanghai Research Center for Quantum Sciences, Shanghai 201315 3Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, USA 4Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute, Shanghai Jiao Tong University, Shanghai 200240 Received 8 April 2020, online 25 April 2020 *Supported in part by the National Key Research and Development Program of China (Grant No. 2016YFA0302001), the National Natural Science Foundation of China (Grant Nos. 11674221 and 11745006), the Shanghai Rising-Star Program, Eastern Scholar Professor of Distinguished Appointment Program, the AFOSR (Grant No. FA9550-16-1-0006), the MURI-ARO (Grant No. W911NF-17-1-0323) through UC Santa Barbara, the NSF China Overseas Scholar Collaborative Program (Grant No. 11429402), and the Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01).
**Corresponding authors. Email: zcai@sjtu.edu.cn; w.v.liu@icloud.com Citation Text: Cai Z, Huang Y Z and Liu W S 2020 Chin. Phys. Lett. 37 050503    Abstract Temperature is a fundamental thermodynamic variable for matter. Physical observables are often found to either increase or decrease with it, or show a non-monotonic dependence with peaks signaling underlying phase transitions or anomalies. Statistical field theory has established connection between temperature and time: a quantum ensemble with inverse temperature $\beta$ is formally equivalent to a dynamic system evolving along an imaginary time from 0 to $i\beta$ in the space one dimension higher. Here we report that a gas of hard-core bosons interacting with a thermal bath manifests an unexpected temperature-periodic oscillation of its macroscopic observables, arising from the microscopic origin of space-time locked translational symmetry breaking and crystalline ordering. Such a temperature crystal, supported by quantum Monte Carlo simulation, generalizes the concept of purely spatial density-wave order to the imaginary time axis for Euclidean action. DOI:10.1088/0256-307X/37/5/050503 PACS:05.30.Jp, 02.70.Ss, 03.65.Yz © 2020 Chinese Physics Society Article Text Introduction.—The phases of quantum matter are, in general, categorized according to symmetries and their spontaneous breaking. A striking possibility was put forward,[1] dubbed "time crystal", where the interacting particles spontaneously organize themselves into a "time" coherent state that breaks temporal translational symmetry. Temporal periodicities spontaneously emerge analogous to spatially ordered crystals. This phase has raised considerable interest,[2–17] and was subsequently proven to be forbidden in thermodynamic equilibrium states.[18,19] A variant of the phase was proposed[6–8] by considering certain non-equilibrium condition, and has quickly been observed in experiments.[20,21] Statistical field theory shows that a $d$-dimensional quantum thermal ensemble is equivalently described by the path integral representing the partition function in $d+1$ Euclidean space where time is imaginary. In the field theory dual description, the thermal distribution weighted by the well-known temperature dependent factor is fully captured by integration over auxiliary imaginary time ("iTime") ranging from 0 to $\beta=1/k_{_{\rm B}} T$ ($k_{_{\rm B}}$ the Boltzmann constant). In this manner, finite-temperature physics is mapped to the dynamics of Euclidean action. Motivated by such a profound relation between temperature $T$ and iTime, Wilczek had speculated the novel possibility of spontaneous translational symmetry breaking in iTime in the end of his original paper, dubbed as "imaginary time crystal" (iTC).[1] One may wonder what class of thermal ensembles can show crystalline order in iTime space, whether the "time" analogue of phonons exists and what novel effects should manifest in ground state and thermal response functions. In this Letter, we find a class of open quantum ensembles that shows the iTC phase. A key feature of our mechanism is a non-monotonic bath-induced retarded interaction, whose dependence is on the relative time in two-body scattering. The feature satisfies the two conditions that may previously seem to be conflicting, i.e., being time dependent but fully translational-invariant. As we will show in the following, integrating out certain engineered local thermal baths yields such a retarded interaction, which favors crystalline patterns in iTime with periods determined by the position of interacting potential minimum. Thermodynamic quantities of the iTC exhibit a striking behavior of oscillation in temperature, which provides an experimental observable effect for this novel phase. Model and method.—Our model is hard-core bosons in either one-dimensional (1D) lattice of $L$ sites with system Hamiltonian: $H_{\rm s}=\sum_{\langle i\rangle} -J[a_i^† a_{i+1}+ h.c.], $ where $J$ denotes the hopping amplitude. Focus is on the half-filling case throughout this study. On each site $i$ the density operator of the hard-core boson $n_i$ additionally couples to a local bath. We assume that the total system (system+bath) is in thermodynamic equilibrium with temperature $T$. The partition function takes the form $Z_{\rm tot}={\rm Tr}_{\rm s}{\rm Tr}_{\rm b} e^{-\beta H_{\rm tot}}$ (${\rm Tr}_{\rm s/b}$ denotes tracing over the system/bath degree of freedom). We first integrate out the bath degrees of freedom, which yields a retarded density-density interaction term in iTime. Thus the open system can be described by the reduced density matrix (see Supplemental Material (SM) for details) $$ \hat{\rho}_{\rm s}= \frac{1}{Z_{\rm s}}e^{-\beta H_{\rm s}-S_{\rm R}},~~ \tag {1} $$ where $H_{\rm s}$ is the system Hamiltonian and $Z_{\rm s}={\rm Tr}_{\rm s}\hat{\rho}_{\rm s}$. $S_{\rm R}$ describes the effective action form of the bath induced retarderd interaction which is onsite in space and translationally invariant in iTime: $$ S_{\rm R}\! =\!-\!\!\int_0^\beta \!\!\!\!d\tau \!\!\int_0^\beta \!\!\!\!d\tau' \!\sum_i [n_i (\tau) \!-\! n_0] D(\tau\!-\!\tau') [n_i(\tau') \!- \!n_0] ,~~ \tag {2} $$ where $n_0= 1/2$ is the average particle density, and the site-independent kernel function $D(\tau-\tau')$ depends on the bath properties. For instance, a harmonic oscillator bath with frequency $\omega_0$ gives rise to the kernel function $D(\tau-\tau')\sim e^{-\omega_0|\tau-\tau'|}+e^{-\omega_0(\beta-|\tau-\tau'|)}$ for large $\beta$. We will show later that the iTC phase can exist for a broad class of kernel functions. Here we choose $$\begin{align} D(\tau-\tau') =\,& \alpha [F(\tau-\tau')+F(\beta-|\tau-\tau'|)], \\ F(\tau-\tau') =\,& e^{-\omega_d|\tau-\tau'|}\cos[2\pi\omega_c (\tau-\tau')],~~ \tag {3} \end{align} $$ where $|\tau|\leq\beta$, and $\alpha$ is the retarded interaction strength. The model with the retarded interactions defined in Eq. (3) can be more than a toy model of merely academic interest. In general, a bosonic bath induced interaction is usually attractive, while a fermionic bath-induced interaction is more complex as sometimes it can be oscillating between attraction and repulsion (e.g., the Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction[22–24]). Now we need an RKKY-like interaction, but with oscillating decay in iTime instead of real space. Consider $$ D(\tau)=\sum_m \frac {e^{i\omega_m\tau}}{i\omega_m-\Sigma_b(i\omega_m)} $$ with $\omega_m={2\pi m}/{\beta}$ the Matsubara frequency and $\Sigma_b(i\omega_m)$ the self-energy of the bath. In general, an imaginary part in self-energy leads to an exponential decay of the quasi-particle in real time (finite lifetime), which corresponds to an oscillation in the imaginary time after analytic continuation. Motivated by Kozii and Fu's recent work about an effective non-Hermitian description for fermionic quasiparticles,[25] in the first section of SM, we propose a similar microscopic bath model that can induce the oscillating retarded interaction as expressed in Eq. (3). Moreover, by analogy with interaction in conventional crystals (e.g., Van der Waals force), we should emphasize that the time-oscillating feature in the interaction is not a necessary condition for the iTC phase. In the second section of SM, we have chosen a different form of retarded interaction and numerically shown that the phase can actually exist as long as the retarded interaction is nonlocal and non-monotonic in iTime with at least one minimum, whose position defines the "time" lattice constant of iTC. Path integral quantum Monte Carlo (QMC) simulation is a stochastic numerical method to study the equilibrium properties of quantum many-body systems based on importance sampling of the world line configurations in Euclidean spacetime. Following the algorithm proposed in Ref. [26], we generalize the QMC simulation with worm-type updates[27] to study the open systems whose distribution functions [Eq. (1)] deviate from the Boltzmann distribution. The key point is to evaluate the integrals resulting from the retardation and include them into the QMC acceptance ratio during the updates of sampling. Notice that since the retarded interaction [Eq. (2)] only appears in the diagonal parts (under the Fock basis) of the effective action, they do not cause extra minus sign problem. Thus our QMC simulations can be considered to be numerically exact since our "system" Hamiltonian (of hard-core bosons) is also positive definite. Groundstate properties.—The bosonic models with retarded interactions can be solved numerical exactly by quantum Monte Carlo simulations (see Methods). We first focus on the ground state of the total system (system+bath). In finite size scaling, the inverse temperature is chosen as $\beta=L$, and the ground state in the thermodynamic limit is approached in the limit $L\rightarrow \infty$. In the 1D case without retardation, the ground state is a quasi-superfluid ground state characterized by algebraic decaying correlation functions in both space and iTime. As the retarded interaction is switched on, spatial-temporal configurations with periodically oscillating world lines are favored since the interacting energies are lowered in these patterns. The competition between the hopping and the retarded interaction leads to novel quantum phases. The crystalline order can be identified from the density correlation functions in space and iTime: $$ S(r,\tau)=\frac 1L\sum_i \langle (n_i(0)-n_0)(n_{i+r}(\tau)-n_0)\rangle. $$ We first focus on the onsite correlations ($r=0$). As shown in Fig. 1(a), for small $\alpha$, $S(0,\tau)$ still decays algebraically with $\tau$, but is accompanied by oscillations. For large $\alpha$ a long-range order is built up and a characteristic time period emerges, indicating a spontaneous breaking of translational symmetry in iTime, characterized by the order parameter $m^2_\tau=\frac 1\beta \int_0^\beta d\tau S(0,\tau) \cos[2\pi\omega_c\tau]$. A finite-size scaling of $m_\tau$ is shown in Fig. 1(c), from which we can find that for large $\alpha$, $m_\tau$ extrapolates to a finite value in the thermodynamic limit, indicating a crystalline order in iTime. The existence of a phase transition can be verified by the correlation length $\xi_\tau$ along iTime ($S(0,\tau)\sim e^{-|\tau|/\xi_\tau}$ for sufficiently large $\tau$), which can be derived from the structure factor (Fourier transformation of $S(0,\tau)$). The normalized correlation length $\xi_\tau/\beta$ as a function of $\alpha$ for different system sizes and $\beta$ is plotted in Fig. 1(e), where we find a crossing point at $\alpha_c^\tau=0.23(1)$, indicating a scaling invariant critical point. In the iTC, $\xi_\tau/\beta$ linearly scales with the system size, and extrapolates to a finite value in the thermodynamic limit (inset of Fig. 1(e)).
cpl-37-5-050503-fig1.png
Fig. 1. Temporal (upper panels) and spatial (lower panels) structures of iTC. (a) Onsite density correlation as a function of $\tau$ and (b) equal-time density correlation as a function of $r$ for various retarded interaction strength $\alpha$ in a system with $L=64$; (c) and (d) finite-size scaling of $m_\tau^2$ and $m_{\rm s}^2$ for different $\alpha$; (e) normalized correlation time and (f) normalized correlation length as a function of $\alpha$. Insets in (e) and (f): finite-size scaling of $\xi_\tau/\beta$ and $\xi_{\rm s}/L$ with $\alpha/J=0.28$. For (a)–(f), the scaling relation is $\beta=L$ and the parameters $\omega_c$ and $\omega_d$ are fixed at $\omega_c=J/4$ and $\omega_d=0.1J$.
For sufficiently large $\alpha$, it is interesting to notice that the crystalline pattern emerges not only in iTime, but also in real space, which can be identified by the equal-time correlation function $S(r,0)$ as well as its structure factor $S(Q)=\frac 1L\sum_r e^{iQr}S(r,0)$. In the case of half-filling, the spatial crystalline pattern is a charge-density-wave (CDW) state with a periodicity twice that of the lattice, as shown in Fig. 1(b). A finite size scaling of the order parameter $m_{\rm s}=\sqrt{S(Q=\pi)}$ is shown in Fig. 1(d). Since there is no direct nearest-neighboring repulsive interaction in our system Hamiltonian, the CDW order is formed by an effective repulsive interaction induced by the retardation: the retarded interaction (3) favors oscillating world lines in iTime, and thus requires the hard-core bosons being separated as far as possible to avoid blocking the oscillations of each other. Therefore at half-filling the tendency towards CDW order is expected in real space. Similar to the above analysis, we also plot the normalized correlation length $\xi_{\rm s}/L$ in real space as a function of $\alpha$ in Fig. 1(f), where a continuous phase transition is found at $\alpha=0.245(5)$. Above this value, long-range CDW correlations emerge in real space. To further investigate the properties of the iTC, we calculate the single particle Green's function: $G(r,\tau)=\langle a_i^†(\tau) a_{i+r}(0)\rangle$. The equal-time part $G(r,0)$ is shown in Fig. 2(a), which decays algebraically with $r$ (accompanied by an oscillation) even for large $\alpha$, indicating a quasi-superfluid order in the iTC. This order can be further verified from the superfluid density, which can be calculated as $\rho_{\rm s}=\frac L\beta\langle W^2\rangle$, where $W$ is the winding number denoting the net times by which the world-lines wrap around the 1D lattice. As shown in the inset of Fig. 2(a), $\rho_{\rm s}$ grows monotonically with $\alpha$, indicating that the retarded interaction enhances the mobility of bosons. Now we turn to the temporal part of $G(r,\tau)$, and focus on its $k=0$ component: $G(k=0,\tau)=\frac {1}{L} \sum_r G(r,\tau)$. Figure 2(b) shows that $G(k=0,\tau)$ decays exponentially with $\tau$, indicating a finite single-particle charge gap. This agrees with the results of the compressibility $\kappa=\frac\beta L(\langle N^2\rangle-\langle N\rangle^2)$ with $N$ the total particle number. Figure 2(c) shows that $\kappa$ decays to zero for large $\alpha$, thus the iTC is incompressible. In summary, the iTC phase induced by the retarded interaction (3) spontaneously breaks the translational symmetries in both real space and iTime. Besides these intriguing properties, this phase is unique in the sense that it is incompressible but possesses quasi-long range superfluid order, so it is fundamentally different from the conventional bosonic ground states in any closed quantum many-body systems: e.g., the (quasi-)superfluidity, Mott insulator, supersolid, Bose glass or Bose metal. The nature of this phase can be qualitatively understood in the strong coupling limit $\alpha\gg J$, where the retarded interaction dominates in the partition function. For each site, the retarded interaction (3) favors configurations with alternating $n(\tau)$ that oscillates between 0 and 1 with a frequency $2\pi\omega_c$, which is reminiscent of the soliton solutions in the BCS pairing model.[28,29] Since a world line needs to be continuous, the temporal-spatial configurations minimizing the retarded interacting action (2) are highly degenerate. As shown in Fig. 2(d), these configurations can be understood as a synchronous movement of the bosons in iTime, in the sense that different bosons hop simultaneously towards the same direction. These synchronous movements can reduce the possibility of the collisions of world lines of different hard-core bosons, and at the same time, minimize the retarded interacting energy. In this hugely degenerate manifold, those configurations with nonzero winding number contribute to the nonzero superfluid density (e.g., the 2nd figure in Fig. 2(d)). The incompressible nature can also be understood in such a picture: an extra particle added into the half-filled system will inevitably block or impede the synchronous movement of other bosons due to its hard-core nature, and thus will increase the retarded interacting energy and give rise to a finite charge gap. Symmetries and excitations.—The iTC phase simultaneously breaks the continuous and discrete translational symmetry in temporal and spatial directions. The spacial and temporal spontaneous symmetry breakings lock with each other, leaving an unbroken subgroup of intertwined space-iTime translational symmetry. The unit cell of the iTC is spanned by the primitive basis vectors ${\boldsymbol b}_{1/2}$ with a rhombic symmetry as shown in Fig. 2(e), as opposed to a "decoupled" direct product of vectors in $\tau$ and $x$ directions with a rectangular symmetry. Due to the fact that the time continuum is locked with the discrete spatial lattice we started with, small translational fluctuations of the atomic positions along the primitive basis vectors do not generate gapless "phonons", but excitations with finite energy cost due to the discreteness of the spatial translational symmetry in our model. The space-iTime locking avoids the dangerous gapless modes in renormalization, so remarkably stabilizes the iTC phase. Temperature crystal.—The most striking effect of iTC can be found at finite temperature. In conventional systems, the temperature dependence of a physical quantity could be either monotonic or non-monotonic with one or several peaks. In contrast, an iTC could exhibit exotic thermal behavior that its physical observable depends on temperature in an oscillating way. This intriguing property can be understood as follows: the crystal structure in iTime will perfectly fit at certain temperatures satisfying the condition that $\beta$ is a multiple of the period of iTC; otherwise, for general $\beta$ the mismatch between them will introduce defects or distortions, which impede the synchronous movement of the bosons, thus are detrimental to crystalline order in both iTime and space. As an example, in Fig. 2(f), we plot the CDW ordering parameter as a function of $\beta$, which shows an oscillating behavior with peaks at the positions of integer multiples of the iTC period. This abnormal thermal behavior can be considered as an experimental diagnostic of the iTC.
cpl-37-5-050503-fig2.png
Fig. 2. (a) Equal-time Green's function $G(r,0)$ for various $\alpha$ in the iTC phase (inset: the superfluid density as a function of $\alpha$) with $L=64$. (b) Zero-momentum component of the unequal-time Green's function $G(k=0,\tau)$ in the iTC phase and various system sizes (inset: the finite-size scaling of the exponent in the exponential decay of $G(0,\tau)\sim e^{-\Delta\tau}$). (c) The compressibility as a function of $\alpha$. (d) Typical temporal-spatial configurations that equally contribute to the partition function in the strong coupling limit $J/\alpha\rightarrow 0$. (e) Sketch of temporal-spatial crystalline structure of the iTC phase and its unit cell (the grey rhombic area spanned by primitive basis vectors ${\boldsymbol b}_{1/2}$). (f) The (inverse) temperature dependence of $m_{\rm s}^2$ in the iTC phase and system size $L=32$. $L=\beta$ for (a)–(c) and $\alpha=0.3J$ for (b) and (f). Here $\omega_c$ and $\omega_d$ are the same as those in Fig. 1.
Experimental realization.—The crucial ingredient to realize the proposed idea is to look for a system that has retarded interaction and allows the system-bath coupling to be easily tuned. Let us elaborate on what was expressed when introducing the model Hamiltonian. Such an interaction can be induced by a fermionic bath with finite life time, which generally arises from fermion self-scattering or from the scattering of bath fermions with impurities or phonons. Instead of the non-Hermitian mechanism as proposed in Ref. [25] (where the effective interaction is also induced by the scattering with phonons albeit they had a completely different purpose), here we propose a strongly interacting ultracold fermionic bath, whose self-energy naturally contains an imaginary part (finite line width in spectrum). Recent breakthrough in an ultracold atomic system[30] has shed light on realization of such an unconventional interaction. In Ref. [30], fermion-mediated interactions have been observed in a bose-fermi mixture in the regime that the dynamics of fermions is much faster than that of the bosons. The difference needed for our model appears to adjust the characteristic time scales of the bosons and fermions to be compatible, so that the time-delayed feature of the bath induced interaction can not be neglected. Conclusion and outlook.—We have introduced a class of quantum many-body systems with oscillating retarded interactions in iTime and performed numerically exact QMC simulations. Our study finds an iTC phase with exotic zero-temperature properties and abnormal thermal behavior. To conclude this study, we wish to outline some differences between our results and other related studies recently. Firstly, our model is conveniently presented by a Euclidean action, so it is beyond the scope of Hamiltonian problems and does not contradict the "no time crystal ground state" theorem.[18,19] Secondly, even though it corresponds to a "time"-dependent problem, the time dependence enters our model through the relative time difference in terms of two-body interactions, in contrast to the early Floquet models whose time-dependence is explicitly through overall time.[6,8] This difference fundamentally changes the symmetry and a continuous symmetry model opens the wonder of what the analogue of Goldstone modes in spatial crystals, i.e, the ubiquitous and all important phonons, will be after the spontaneous breakdown of time translational symmetry. The potential of discovering the time phonons is now open for further work. Finally, a thermodynamic physical system is found to support our "theoretical" time-crystal model. Standard statistical physics is applied to show the correspondence of our imaginary-time quantum mechanical model with retarded interaction (dependent on relative time!) to the ensemble of hard-core bosons embedded in engineered baths. This report therefore extends the study of time crystal to the Euclidean spacetime corresponding to thermal ensemble systems. Our method can be straightforwardly generalized to other bosonic and spin systems with retarded interactions, in one or higher dimensions. The competition between time-oscillating interaction and quantum-thermal fluctuations is expected to give rise to other forms of nontrivial temporal-thermal orderings that are completely absent in the conventional Hamiltonian systems. A generalization to fermionic models may be even more interesting since it opens up new possibilities of finding novel "non-fermi liquid". It is also interesting to study the real time counterpart of our model, which may be derived by integrating out a non-equilibrium bath. This may provide a new perspective for current time crystal researches since such an effective action has the full continuous (rather than discrete) symmetry in temporal translation. Acknowledgment.—We appreciate insightful discussions with F. Wilczek and V. Galitski. References Quantum Time CrystalsClassical Time CrystalsSpace-Time Crystals of Trapped IonsSuperfluidity and Space-Time Translation Symmetry BreakingModeling spontaneous breaking of time-translation symmetryFloquet Time CrystalsPhase Structure of Driven Quantum SystemsDiscrete Time Crystals: Rigidity, Criticality, and RealizationsTime Crystal Behavior of Excited EigenstatesDefining time crystals via representation theoryFloquet time crystal in the Lipkin-Meshkov-Glick modelDiscrete Time-Crystalline Order in Cavity and Circuit QED SystemsClean Floquet Time Crystals: Models and Realizations in Cold AtomsTime crystals: a reviewBoundary Time CrystalsQuantum Time Crystals from Hamiltonians with Long-Range InteractionsTime-crystalline Topological SuperconductorsImpossibility of Spontaneously Rotating Time Crystals: A No-Go TheoremAbsence of Quantum Time CrystalsObservation of discrete time-crystalline order in a disordered dipolar many-body systemObservation of a discrete time crystalIndirect Exchange Coupling of Nuclear Magnetic Moments by Conduction ElectronsA Theory of Metallic Ferro- and Antiferromagnetism on Zener's ModelMagnetic Properties of Cu-Mn AlloysNon-Hermitian Topological Theory of Finite-Lifetime Quasiparticles: Prediction of Bulk Fermi Arc Due to Exceptional PointIdentifying a Bath-Induced Bose Liquid in Interacting Spin-Boson Models“Worm” algorithm in quantum Monte Carlo simulationsInstanton Sector of Correlated Electron Systems as the Origin of Populated Pseudo-gap and Flat “Band” Behavior: Analytic SolutionNonperturbative quantum dynamics of the order parameter in the BCS pairing modelObservation of fermion-mediated interactions between bosonic atoms
[1] Wilczek F 2012 Phys. Rev. Lett. 109 160401
[2] Shapere A and Wilczek F 2012 Phys. Rev. Lett. 109 160402
[3] Li T, Gong Z X, Yin Z Q, Quan H T, Yin X, Zhang P, Duan L M and Zhang X 2012 Phys. Rev. Lett. 109 163001
[4] Wilczek F 2013 Phys. Rev. Lett. 111 250402
[5] Sacha K 2015 Phys. Rev. A 91 033617
[6] Else D V, Bauer B and Nayak C 2016 Phys. Rev. Lett. 117 090402
[7] Khemani V, Lazarides A, Moessner R and Sondhi S L 2016 Phys. Rev. Lett. 116 250401
[8] Yao N Y, Potter A C, Potirniche I D and Vishwanath A 2017 Phys. Rev. Lett. 118 030401
[9] Syrwid A, Zakrzewski J and Sacha K 2017 Phys. Rev. Lett. 119 250602
[10] Khemani V, von Keyserlingk C W and Sondhi S L 2017 Phys. Rev. B 96 115127
[11] Russomanno A, Iemini F, Dalmonte M and Fazio R 2017 Phys. Rev. B 95 214307
[12] Gong Z, Hamazaki R and Ueda M 2018 Phys. Rev. Lett. 120 040404
[13] Huang B, Wu Y H and Liu W V 2018 Phys. Rev. Lett. 120 110603
[14] Sacha K and Zakrzewski J 2018 Rep. Prog. Phys. 81 016401
[15] Iemini F, Russomanno A, Keeling J, Schirò M, Dalmonte M and Fazio R 2018 Phys. Rev. Lett. 121 035301
[16] Kozin V K and Kyriienko O 2019 Phys. Rev. Lett. 123 210602
[17] Chew A, Mross D F and Alicea J 2019 arXiv:1907.12570 [cond-mat.mes-hall]
[18] Bruno P 2013 Phys. Rev. Lett. 111 070402
[19] Watanabe H and Oshikawa M 2015 Phys. Rev. Lett. 114 251603
[20] Choi S, Landig R, Kucsko G, Zhou H, Isoya J, Jelezko F, Onoda S, Sumiya H, Khemani V, von Keyserlingk C 2017 Nature 543 221
[21] Zhang J, Hess P W, Kyprianidis A, Becker P, Lee A, Smith J, Pagano G, Potirniche I D, Potter A C, Vishwanath A 2017 Nature 543 217
[22] Ruderman M A and Kittel C 1954 Phys. Rev. 96 99
[23] Kasuya T 1956 Prog. Theor. Phys. 16 45
[24] Yosida K 1957 Phys. Rev. 106 893
[25] Kozii V and Fu L 2017 arXiv:1708.05841 [cond-mat.mes-hall]
[26] Cai Z, Schollwöck U and Pollet L 2014 Phys. Rev. Lett. 113 260403
[27] Prokof'ev N V, Svistunov B V and Tupitsyn I S 1998 Phys. Lett. A 238 253
[28] Mukhin S 2009 J. Supercond. Novel Magn. 22 75
[29] Galitski V 2010 Phys. Rev. B 82 054511
[30] DeSalvo B, Patel K, Cai G and Chin C 2019 Nature 568 61