Persistent dark states in anisotropic central spin models.

Long-lived dark states, in which an experimentally accessible qubit is not in thermal equilibrium with a surrounding spin bath, are pervasive in solid-state systems. We explain the ubiquity of dark states in a large class of inhomogeneous central spin models using the proximity to integrable lines with exact dark eigenstates. At numerically accessible sizes, dark states persist as eigenstates at large deviations from integrability, and the qubit retains memory of its initial polarization at long times. Although the eigenstates of the system are chaotic, exhibiting exponential sensitivity to small perturbations, they do not satisfy the eigenstate thermalization hypothesis. Rather, we predict long relaxation times that increase exponentially with system size. We propose that this intermediate chaotic but non-ergodic regime characterizes mesoscopic quantum dot and diamond defect systems, as we see no numerical tendency towards conventional thermalization with a finite relaxation time.

Introduction

State-of-the-art quantum technologies can control and coherently manipulate qubit systems with exquisite precision[1-4]. The surrounding environment of the qubit however eventually decoheres the qubit and limits quantum applications[5-7]. An efficient way of extending coherence times is to prepare the system in so-called dark states, in which the qubit is effectively decoupled from the bath[8-10]. Dark states have been identified in several integrable central spin models[8,10], and are central to quantum computing[11,12], metrology[13,14] and control[15,16] applications in a variety of experimental qubit systems, including nitrogen vacancy centers in diamond[17,18] and semiconducting quantum dots[19,20]. A central goal of this work is to show that dark states can persist in experimentally relevant non-integrable central spin models. At numerically accessible system sizes, they exist as exact eigenstates. In the thermodynamic limit, the qubit could eventually thermalize but only after long times.

Central spin systems are typically described by a spin-1/2 Hamiltonian ():where is a local magnetic field on the central qubit, is a uniform magnetic field on the bath spins, sets the anisotropy of the qubit-bath interaction, and sets the strength of the interaction between the central qubit and the ith bath spin for . Experimentally, the interaction strengths are inhomogeneous because of the randomness in the positions of the bath spins and/or the geometrical factors in dipolar interactions. For simplicity, we model the inhomogeneity as uncorrelated disorder, and take the to be independently and identically distributed uniformly in the interval with setting the disorder strength. Moreover, we study the model near resonance where qubit-bath interactions are enhanced. Since H conserves total spin magnetization , we set without loss of generality. Fig. 1 shows a schematic of the model.

Figure 1

Schematic of the spin-1/2 anisotropic central spin model. A central qubit in a magnetic field of strength interacts with an environment of spin-1/2 particles in a uniform magnetic field of strength with interaction strengths , .

The Hamiltonian H has three known integrable families. The first is the fully isotropic XXX model ( and arbitrary ), which describes systems with contact interactions such as quantum dots in s-type semiconductor bands[19,20]. This model belongs to the class of integrable XXX Richardson-Gaudin models[21-23]. The second is the fully anisotropic XX model ( and arbitrary ), which describes resonant exchange interactions in dipolar spin systems in rotating frames[24-27]. It was only recently established that the XX model is integrable, arising as a singular limit of the class of hyperbolic XXZ Richardson-Gaudin models[10]. The third is the homogeneous XXZ model ( and arbitrary ), which describes effective two-body interactions between the qubit and the collective spin of the bath [28]. Fig. 2 shows the three integrable families in a broader phase diagram. The integrability of these models has enabled analytical and numerical studies of experimentally relevant systems using a variety of integrability-based techniques[20,29-37].

Figure 2

Finite size crossover diagram. Integrable lines with product dark eigenstates are shown as dashed black lines, while the integrable XXX line with no product dark eigenstates is shown in blue. In the chaotic regime between integrable lines, we show a color plot of the central spin polarization , averaged over eigenstates with the smallest values of and disorder samples in a fixed sector of total magnetization at resonance . In this regime, non-thermal persistent dark states with (violet region) coexist with bright states whose central spin polarization is close to zero. While this crossover diagram shows a fixed system size , we find no significant dependence on L for system sizes amenable to numerical simulation.

A remarkable feature of H is that it exhibits dark eigenstates of a particularly simple product form when either or [10,28]. These product dark states are states which exhibit no qubit-bath entanglement. Namely, they have a product state structure or , in which the central spin is polarized along the z-direction, and the state of the bath satisfiesThese states are furthermore independent of (and ), making the qubit state insensitive to external axial fields, in addition to bath fluctuations. More generally, we define dark states as states in which the qubit is nearly z-polarized and is not in thermal equilibrium with the surrounding bath. Dark states allow the surrounding spin bath to be used as a robust quantum memory[8,9,38,39]. Moreover, they pose limitations for dynamical nuclear polarization (DNP) experiments which attempt to polarize a mesoscopic bath by repeated qubit polarization and transfer[28,40-42]. DNP protocols eventually prepare the system in a statistical mixture of dark states, producing effectively isolated qubits for decoherence-free quantum computation[43].

In this work, we first establish that dark states are robust to integrability-breaking perturbations that tune the anisotropy and the disorder strength . Specifically, we show that dark states persist as exact eigenstates of the Hamiltonian, that only perturbatively mix with bright states (i.e. not dark states) over a broad range of values for and , at system sizes amenable to numerical simulation. This perturbative mixing only slightly reduces the polarization of the qubit along the z-direction. Remarkably, while H is non-integrable/chaotic away from its integrable lines, dark states are well protected due to the presence of quasi-conserved charges.

To test the stability of dark states, we apply a recently developed exponentially sensitive probe for chaos, based on the scaling of the norm of the adiabatic gauge potential (AGP)[44]. The AGP is defined as the operator which generates continuous adiabatic transformations between eigenstates and measures their sensitivity to perturbations of the underlying Hamiltonian[45-49]. Its norm is closely related to the quantum geometric tensor and the fidelity susceptibility[45,50,51]. The norm of the AGP was found to scale exponentially with system size for chaotic perturbations, in contrast to integrable perturbations leading to polynomial scaling[44].

In our present context, the AGP norm grows exponentially in accordance with quantum chaos , but interestingly, the growth rate of the logarithm of this norm is twice the rate expected for ergodic systems satisfying the eigenstate thermalization hypothesis . This rate saturates the upper bound for eigenstate sensitivity to perturbations[44]. It reflects a very strong mixing between neighbouring eigenstates of the system and leads to ultra-slow (exponentially long in system size) relaxation dynamics , reminiscent of the Arnold diffusion in classical near-integrable systems[52]. While a similar behavior of the AGP norm was found in previous work[44] for spin chains with weak integrability breaking perturbations, here we find that this chaotic but non-ergodic (CNE) regime extends to large perturbation strengths, even for the largest system size that we are able to simulate.

Results

Persistent dark states

Away from the integrable lines in Fig. 2, the eigenstates of H no longer admit exact product dark states. Nevertheless, we can identify persistent dark eigenstates with approximate product form using (i) the eigenstate expectation value of the central qubit z-projection , or (ii) the eigenstate entanglement entropy of the qubit. The latter is defined aswhere is the reduced density matrix for the qubit obtained by tracing out the bath degrees of freedom, and is the density matrix of eigenstate with energy .

Consider for reference the XX model () at resonance () in a sector with negative net magnetization ()[10]. For any product dark eigenstate of this model, and . On the other hand the bright eigenstates of the model satisfy and at resonance (far from resonance, the central spin is nearly polarized even in the bright states: and ).

For , we find eigenstates of H that are adiabatically connected to by a unitary transformation generated by the adiabatic gauge potential (AGP) :In these eigenstates, both the z-projection and entanglement entropy of the qubit will deviate from their values. The question becomes whether these deviations are perturbatively small in , and how this depends on the system size L. In chaotic systems, the AGP is generally a highly non-local many-body operator with an exponentially large norm[44,45,53]. In the present context, the parameter breaks the integrability of the system (see section below on the chaotic but non-ergodic regime). Naively, we expect qubit observables in to perturbatively connect to their values in only for that is exponentially small in the system size.

Remarkably, at numerically accessible system sizes, we find that can be well-approximated by few-body operators and that perturbation theory works exceedingly well to characterize qubit observables in . To illustrate, consider the perturbative expansion of the expectation to leading order in :The leading term is of order , as the coefficient of the linear in term, , vanishes because . Fig. 3 numerically demonstrates that the left-hand side scales as for a subset of the eigenstates over several orders of magnitude of and .

Figure 3

Dark states persist away from the integrable lines at finite size. (a) Upper panel: Expectation value of the central spin z-projection for every eigenstate of H in a typical sample as a function of . Persistent dark (black circles) and bright (red diamonds) states are easily distinguished by their value of . Dotted lines (gray) show scaling, while the horizontal dashed line (red) indicates . Inset: System size dependence of averaged over disorder samples and the eigenstates with smallest z-projection: . (b) Lower panel: (markers) as a function of for several values of . The solid lines plot the perturbative prediction of Eq. (5). Inset: Upon re-scaling the vertical axis by , the curves collapse onto a single curve. Parameters: , , , , and in (a) .

In more detail, when , the resonance condition in a given polarization sector is shifted by the mean anisotropy to (see Supplemental Information). For concreteness, we focus on a single polarization sector such that the resonance occurs at ; in the Supplemental Information, we generalize to other polarization sectors. Fig. 3(a) shows numerical computations of the expectation value in every eigenstate of H at moderate disorder strength () over several orders of magnitude in . Persistent dark states (black/dark circles) are easily identifiable, as they connect smoothly to as . The deviation from scales as , consistent with Eq. (5) (dotted lines). The bright states (red/light diamonds) are similarly perturbed around their value at resonance (dashed red horizontal line). As , dark and bright states attain comparable central spin projections. The inset of Fig. 3(a) shows the system-size dependence of the averaged expectation value , where denotes an average over all persistent dark states, and denotes an average over disorder samples. In eigenstates that satisfy the ETH, the expectation value approaches zero with increasing L, in a sector with fixed magnetization. However, we find that approaches its thermal value only very slowly with system size L, suggesting that dark state properties persist to system sizes much larger than we probe here. From the current analysis we cannot conclude whether or not they survive the thermodynamic limit.

Fig. 3(b) shows for varying disorder strengths . The markers show numerical data and the solid lines show the analytic predictions given by Eq. (5) up to . Again, we find leading order perturbation theory to be in excellent agreement with numerical simulations for and the entire range between the integrable points and . As perturbation theory begins to break down (see line in plot). When , perturbation theory breaks down much faster at (see Supplemental Information).

Persistent dark states are well captured by perturbation theory due to the quasi-locality of at numerically accessible system sizes. To see this, we decompose the AGP into k-body operators:Here with denote the Pauli basis operators on site , where for every . Throughout this work, we define the norm of any operator by its normalized Frobenius norm:We find has non-zero weight only for k-body operators with . Moreover, the total weight of k-body operators decays as for , so that is well-approximated by 3-body operators. In the Supplemental Information, we showcase the quasi-locality of and estimate .

One can similarly find persistent dark states based on translations in the -parameter space, as in Eq. (4), with a different adiabatic gauge potential . We find perturbatively accessible persistent dark states away from the line at numerically accessible system sizes. The inset of Fig. 3(b) shows the re-scaled averaged expectation value vs at resonance. The data collapse for shows that:at small , in perfect agreement with the perturbative result. At larger values of (), we see deviations from the perturbative result as . Persistent dark states were previously found by mapping exact product dark eigenstates from the homogeneous isotropic limit to the inhomogeneous isotropic regime [8,41]. Our results perturbatively extend dark states into a broader region of parameter space at finite size (see Fig. 2).

Persistent dark states can also be identified by their low central spin entanglement entropy (see Fig. 4). For moderate to low disorder (upper panels with ), low entanglement (dark) states persist and do not fully mix with high entanglement (bright) states, even as . At sufficiently large disorder (lower panels with ), the persistent dark state picture breaks down as , since most states acquire large central spin entanglement .

Figure 4

Low central spin entanglement entropy reveals persistent dark states. The entanglement entropy of the central spin in eigenstates vs the re-scaled energy in typical samples at low (top panels) and large (bottom panels) values of . In the top panels, the persistent dark states have low central-spin entanglement at all values. At larger disorder (bottom panels), the entropy approaches for all eigenstates as . Parameters: , , , is the maximum energy of H in the given total magnetization sector.

Quasi-conserved operators

The question of whether and how systems thermalize is a fundamental one in quantum statistical mechanics. Steady states of integrable systems typically have non-thermal correlations due to the presence of extensively many conserved quantities, and are described by Generalized Gibbs Ensembles (GGEs) that account for these conserved quantities[54-56]. Generic integrability-breaking perturbations usually yield Hamiltonians which are chaotic and satisfy the Eigenstate Thermalization Hypothesis (ETH)[57,58]. Nevertheless, the integrable Hamiltonian can control the approach to a long-lived pre-thermal state when the strength of the integrability-breaking perturbation is sufficiently small[59].

In this section, we establish that the central spin model in Eq. (1) has long-lived non-thermal states controlled by the XX and XXZ integrable lines at accessible system sizes. Specifically, we show that H has approximate conservation laws that persist away from the integrable lines, giving rise to non-thermal correlations in local qubit observables. A simple way to detect the non-thermal correlations in quench experiments is through observables, such as , that differentiate between dark and bright states. As takes non-thermal values (close to ) in the dark state manifold, we find quenched steady states that retain memory of the initial z-polarization of the central spin.

The integrable lines of the model (, , and ) constitute families of Richardson-Gaudin models with extensive numbers of bilinear two-body conserved charges [10,60]. Upon breaking integrability, there no longer exists an extensive number of exactly conserved charges. Instead we find an extensive number of quasi-conserved charges, which very nearly commute with H. To find such quasi-conserved charges, we numerically construct an exhaustive set of few-body operators . These operators are conserved iff . The quasi-conserved charges are those operators with very small ratio: , where sets a typical energy scale.

We construct using the ansatz:and restrict to a complete set of m trace-orthogonal one and two-body spin-1/2 operators with unit norm. To avoid cumbersome notation, we leave the dependence of and on k implicit. We further set . To determine the coefficients , we solve the eigenvalue problem:where and we take . The eigenvectors of M then yield through Eq. (9) a set of orthogonal and bilinear operators with known decay properties. Specifically, the eigenvalue equals the norm of the commutator:To connect to the operator decay rate, consider the short-time expansion of the symmetrized unequal time correlator of at infinite temperature[61,62]:The correlator’s decay rate is thus given by . Hence, if , the unequal-time correlator equals one for all t. If instead , where sets a typical decay rate, then the correlator is close to one for a long time and is approximately conserved up to this time.

Away from the integrable lines, we generally find three kinds of eigenvalues : a few zero eigenvalues, large eigenvalues, and an extensive number of eigenvalues with . Zero eigenvalues correspond to exactly conserved charges related to known conservation laws, while the extensive number of small positive eigenvalues can be identified with quasi-conserved charges.

Fig. 5 shows the smallest eigenvalues (re-scaled by ) obtained by numerically solving Eq. (10) as a function of the index k of the corresponding operator . The right panel shows the eigenvalues at several system sizes () for a fixed disorder strength . We find 4 eigenvalues which are zero within numerical accuracy (shaded gray region), corresponding to exactly conserved charges; namely H, , , . We also see a cluster of intermediate eigenvalues corresponding to quasi-conserved charges, which are separated by a gap from a set of larger eigenvalues (only a small fraction of this set is shown). The vertical dashed lines mark the indices of the quasi-conserved operators with largest eigenvalue for each system size. These maximal indices increase in precise proportion to L, showing that the number of quasi-conserved charges is extensive. Furthermore, the eigenvalues themselves remain of the same order of magnitude for the different system sizes, indicating that the operator lifetimes have no significant system size dependence up to .

Figure 5

An extensive number of quasi-conserved quadratic charges persist upon breaking integrability. The smallest eigenvalues of Eq. (10) in dimensionless units at fixed for a typical sample at different disorder strengths (left) and at different L for fixed (right). Values in the shaded region are zero within numerical accuracy. The left panel shows that the dimensionless decay rate of the quasi-conserved operators increases with , while the right panel shows that the number of quasi-conserved operators is extensive (vertical lines denote the largest quasi-conserved index for each L). Parameters: , (left), (right).

The left panel of Fig. 5 shows the smallest eigenvalues of Eq. (10) for several values of disorder strength at . As expected, we find exactly conserved charges as we approach the integrable line (see data within gray region). On increasing , only 4 charges remain exactly conserved, while the remaining charges become quasi-conserved. The lifetime () of the quasi-conserved charges furthermore systematically decreases with increasing . Previous studies have found similar long-lived quasi-conserved charges in a family of near-(Richardson-Gaudin)-integrable spin models with all-to-all interactions[63].

The extensively many two-body quasi-conserved charges control long-lived non-thermal states in quench experiments. The top panel of Fig. 6 shows the relaxation of to a non-thermal value in a typical sample at moderate and large disorder strength. The system is initialized in a polarization sector far from resonance () in the mixed statewith the bath spins at infinite temperature and the central spin maximally polarized along . The top panel plots following a quench to resonance (). We observe a fast decay to a positive value that is different from the thermal value . Thus, retains memory of its initial condition at these system sizes. This memory is a consequence of the weight of on the persistent dark state manifold.

Figure 6

Central spin z-polarization retains memory in the pre-thermal state. Top: time evolution of for a quench to resonance () from the initial density matrix (Eq. (13)) in a typical sample at two different disorder strengths. The late time values (horizontal black dash-dotted lines) differ from the thermal value (horizontal dashed red line). Bottom: The ratio of the disorder-averaged long-time value to the long-time value on the integrable line vs L. The ratio is close to one with no significant finite-size flow at moderate , but decreases with increasing L at large . Parameters: (top), (bottom), , .

Two comments are in order. First, the hybridization between the dark and bright state manifolds increases with disorder strength. Consequently, at any given , in the persistent dark state manifold decreases with increasing (see Fig. 3). This explains why decreases with increasing in Fig. 6. Next, the initial decay in the top panel of Fig. 6 is a consequence of dephasing between the perturbed bright states. As in each perturbed bright eigenstate, the weight of on the perturbed bright states does not contribute to the non-zero value of .

The lower panel of Fig. 6 plots the re-scaled and disorder-averaged long time value with L. The re-scaling factor is the expected polarization of the central spin due to the weight of on the dark manifold,where and are respectively the number of dark and bright states with the central spin pointing along in the appropriate polarization sector at large . On the integrable XX line at , we expect that . The blue (filled) curve in the lower panel of Fig. 6 is perturbatively accessible from the integrable line at the numerically accessible system sizes, and thus we find that is close to one. At larger disorder strength however, the hybridization between the dark and bright states increases with L at the accessible sizes. The long time value is thus smaller than the long time value at , with the discrepancy growing with L (see orange curve with open markers). At large disorder strength, shows a trend toward the thermal value with increasing L. This thermal value follows from and as in the sector with total magnetization; given a total magnetization density , as . The ratio is less than one for all . At the given system sizes, we cannot determine with certainty whether the qubit will saturate at or before it reaches its thermal value, even in the presence of strong disorder.

We conclude this section with two remarks. First, in addition to , generic two-body observables with significant overlap on the quasi-conserved charges are expected to exhibit similar non-thermalizing behavior and non-thermal eigenstate expectation values. While it is possible that in the thermodynamic limit all such observables will thermalize, there is no indication that this will happen from available data at small or intermediate disorder. Even if it happens, the non-thermal state after the quench could crossover to an extremely long-lived and very stable prethermal regime.

Chaotic but non-ergodic regime

In previous sections, we established that dark eigenstates persist on adding putative integrability-breaking perturbations to Eq. (1) at numerically accessible system sizes. At these sizes, the model thus does not satisfy the ETH and few-body observables do not thermalize in isolation. However, we expect the eigenstate and dynamical behavior to change with increasing L. In this section, we provide evidence that the model is in a chaotic non-ergodic regime (CNE) characterized by an exponential sensitivity of eigenstates to small perturbations and the presence of relaxation times that are exponentially long in the system size.

Energy level statistics are a widely used tool to diagnose chaos and predict thermalization[64-66]. Integrable systems generally follow Poisson level statistics, while chaotic systems exhibit Wigner-Dyson statistics due to level repulsion in accordance with random matrix theory[57]. The use of level statistics to diagnose chaos is limited to relatively small system sizes where exact diagonalization can be reasonably implemented. For our present model, level statistics show weak-to-negligible level repulsion, thus proving insufficient to establish chaos.

Fig. 7 shows the distribution P(r) of the ratio r of consecutive energy level spacings in a sector with fixed polarization. The ratio for the trio of energy levels with energies is defined as , where and the energy levels are ordered [66,67]. The left panel shows that P(r) agrees with that expected for a Poisson spectrum near/on the integrable lines points with , , and . The center panel shows that the Poisson behavior persists in the presence of moderate anisotropy () and disorder () at the largest size we access numerically. Only when the disorder strength is much larger than one () do we see some level repulsion with a weak trend towards Wigner-Dyson statistics with increasing L (right panel). Fig. 7 therefore shows no tendency of the model to become chaotic with increasing L at moderate values of and .

Figure 7

Level-spacing ratio distributions. Left: In the vicinity of the integrable lines, P(r) agrees with that expected for a Poisson spectrum (red line). Center: Distributions remain indistinguishable from the Poisson one (red line) at moderate and when the system is no longer expected to be integrable. Right: At large disorder strength, we see level repulsion and a weak trend toward the Wigner-Dyson distribution (black line) with increasing L. Parameters: , , , (left, center), (right), (right). Energies are sampled in the middle two quartiles of the spectrum.

Recently, Ref.[44] proposed the norm of the AGP as a highly sensitive probe for chaos. Related measures were also proposed earlier in the context of many-body localization, see e.g. Refs.[68,69]. Chaos manifests in the exponential scaling of the Frobenius norm of the AGP with system size, which can be interpreted as an exponential sensitivity of the eigenstates to perturbations of the Hamiltonian. In contrast, integrable perturbations show polynomial scaling[44].

For any Hamiltonian , the AGP operator can be represented as[44,70]:where . In the energy eigenbasis of , the off-diagonal matrix elements of read:where . Note that the diagonal matrix elements , which is a gauge choice[45]. The (scaled) Frobenius norm is then given by:In chaotic systems, fluctuates wildly with L when we take the limit because the terms with the smallest energy differences dominate the sum in the norm. This is a standard manifestation of the problem of small denominators[71]. Instead of taking the limit , it is convenient to set as a regulator. This regulator provides a two-fold advantage: (i) it suppresses the wild fluctuations of the norm with system size, and (ii) it allows us to retain the exponential sensitivity of the AGP norm to small perturbations if we pick , where c is a system-size-independent constant. The regulator is thus parametrically larger than the level spacing, while maintaining a small deviation from the exact AGP[44]. Physically, plays the role of a cutoff time for operator growth in Eq. (15). By picking this time to be exponentially large in L, we probe the sensitivity of eigenstates to infinitesimal perturbations.

For systems satisfying the ETH[57] the states at energy density corresponding to infinite temperature satisfywhere is a random variable with zero mean and unit variance, and is a smooth function that is proportional to the Fourier transform of the correlation function at the frequency . In general, the function f also depends on the average energy . However, as the summation in Eq. (17) is dominated by the eigenstates corresponding to infinite temperature, we suppress this additional dependence. Typically, increases as decreases until becomes comparable to the slowest relaxation time scale in the system (such as the Thouless time) and saturates for smaller values of .

Interestingly, it was recently observed that the function can be defined and remains smooth even in generic integrable systems[72]. Then vanishes as for deformations of the Hamiltonian along integrable directions[44,73-75].

Combining Eq. (17) at finite and Eq. (18) at exponentially small scales , we get the following estimate for the AGP norm:Thus if our system were to satisfy ETH,up to polynomial corrections. In contrast, we would expect only a polynomial scaling of with the system size in an (interacting) integrable system for deformations of the Hamiltonian along an integrable direction[44].

Figure 8 shows the exponential divergence of the disorder-averaged norm of the AGP corresponding to perturbations of the anisotropy parameter at moderate disorder strength in the polarization sector. At , scales polynomially with system size L. Away from the integrable point, the scaling of is polynomial until a critical length , which marks the onset of exponential growth and thus chaos (see vertical dash-dotted lines). The critical length increases with decreasing as:such that becomes infinite in the integrable limit . The power is found using linear regression (see inset).

Figure 8

The exponential divergence of the adiabatic gauge potential norm shows signatures of chaos. Plot shows the disorder-averaged norm as a function of system size L. Dotted lines show the scaling behavior in the chaotic non-ergodic regime . Vertical dashed-dot lines mark the onset of exponential growth at . Inset: vs. , and a regression line whose slope is numerically found to be . Parameters: , , , , , with .

For L larger than the critical length , the norm of the AGP in Fig. 8 grows exponentially at twice the rate predicted by the ETH:From Eq. (19), we obtain . It follows from Eq. (18) that the off-diagonal matrix elements are not exponentially suppressed with system size in the narrow energy interval , in contrast with the ETH prediction. The absence of an exponential suppression in the off-diagonal matrix elements is shown in Fig. 9. The figure shows nearest neighbor (in energy) matrix elements for , averaged over disorder samples and eigenstates ; we denote this averaging by the double-overline . These matrix elements differ by several orders of magnitude between low () and high () disorder. However, no exponential decay with L is observed at either disorder strength at numerically accessible systems sizes. This absence of exponential suppression can only persist when polynomially many nearby eigenstates mix. In contrast, ETH would require that a given eigenstate mix equally with exponentially many nearby eigenstates upon perturbing the system.

Figure 9

Nearest neighbor matrix elements of do not decay with system size L. Plot shows nearest neighbor (NN), next nearest neighbor (NNN), and next next nearest neigbor (NNNN) matrix elements (squared) averaged over eigenstates and disorder. Green filled markers () and orange open markers () show no decay with system size in the chaotic non-ergodic regime. Parameters: , , , . Eigenstates are sampled in the energy window , where is the maximum energy of H in the given magnetization sector.

The behavior of the chaotic non-ergodic regime is manifest separately within the dark and bright manifolds, and jointly in the interactions between these manifolds. We plot the various contributions to the disorder-averaged AGP norm (AGP) in Eq. (17) from dark and bright classes of eigenstates in Fig. 10. The DD (BB) contribution comes from terms involving matrix elements between dark eigenstates (bright eigenstates in the same band; the bright states in the XX model come in Landau-Zener (LZ) pairs that can be continuously followed as a function of in each magnetization sector. These bright states form two bands, consisting of the positive and negative energy states of each LZ pair, respectively. See the Supplemental Information). The matrix elements between the dark and bright states contribute the DB piece, while the matrix elements between the two bright state bands contribute the Landau-Zener (LZ) piece. At and , the sum in the AGP norm is dominated by the intra-band bright-bright (BB) matrix elements. Consequently, the off-diagonal matrix elements between neighboring bright states are not exponentially suppressed (see the discussion below Eq. (22)). Dark-dark contributions also exponentially increase with L; however their total value is many orders of magnitude smaller than the BB contribution at these sizes. The DB contributions show a striking difference between the left and right panels of Fig. 10 at the accessible sizes. The DB contributions only grow exponentially with L in the right panel with , reflecting the strong hybridization between neighboring dark and bright states in the spectrum on perturbing (see Supplemental Information). For , on the other hand, the dark and bright state manifolds are separated in energy at the accessible sizes. This limits the hybridization between the two manifolds. However, we expect that the DB contribution diverges exponentially with L at sufficiently large sizes at any . We remind the reader that perturbation theory worked exceedingly well to characterize qubit observables in dark states, even in parameter regimes where the DB contribution exponentially increases with L. This suggests that the strong DB mixing should primarily affect bath observables in persistent dark states at these sizes. In the next section, we discuss the potential implications of this behavior in the thermodynamic limit.

Figure 10

Decomposition of the adiabatic gauge norm into classes based on the nature of the eigenstates. The disorder-averaged AGP norm (AGP) is composed of contributions from matrix elements between dark states (DD), bright states in the same band (BB), bright states in different bands (LZ), and dark and bright states (DB). Left: At , BB and DD show onset of non-ergodic chaos. The AGP norm is dominated by BB, whereas DD is negligible. As DB is relatively constant with L, dark states do not progressively hybridize with bright states at accessible sizes. Right: At , the AGP norm is still dominated by BB. However, dark states hybridize strongly with nearby bright states due to the exponential divergence of DB contribution. Parameters: , , , , cutoff with (left) and (right). The values of c in the two panels are chosen such that .

In the Supplemental Information, we show that the analysis in this section extends to a polarization sector with non-extremal density ,. In particular, we find no deviation from the exponential scaling in Eq. (22).

In sum, the exponential divergence of provides strong evidence that arbitrarily small in perturbations are integrability-breaking, with a growth rate that is twice that predicted by the ETH at numerically accessible system sizes. In this chaotic non-ergodic regime, eigenstates in exponentially small shells of order hybridize, with interaction matrix elements that show no exponential suppression with system size.

Possible fates of the CNE regime in the thermodynamic limit

We have provided evidence that the family of central spin models in Eq. (1) is chaotic, but non-ergodic at moderate disorder strengths and numerically accessible system sizes. Eigenstate chaos manifests in the exponential scaling of the AGP norm with L, while the non-ergodicity is manifest in both the non-thermal value of the central spin polarization in persistent dark states and the non-ETH scaling of the AGP norm.

Ref.[44] argued that the non-ETH scaling of the AGP norm in Eq. (22) indicates slow relaxation with exponentially long in L relaxation times. We repeat the argument for completeness. From Eq. (19) we find:As , is proportional to the relaxation time of the operator . Using that and Eq. (21), we find that for :where and the constant C can have a weak (power law) dependence on L. As the DD, DB and BB components of the AGP norm exhibit non-ETH scaling in the right panel of Fig. 10, we expect that Eq. (24) characterizes certain relaxation processes in both the dark and bright sectors.

In the dark eigensector, the exponentially long relaxation times largely arise from the dark-bright mixing (cf. Fig. 10) and coexist with a robust non-thermal value of the central spin magnetization. This suggests the following ‘cartoon’ for the decomposition of the persistent dark states in the eigenbasis of the XX integrable model,with a non-zero central spin residue . Above, we use the eigenbasis of the XX model at resonance () in a magnetization sector with fixed positive value, is a dark states satisfying Eq. (2), and is a normalized superposition of bright states with the property . Note: In the XX model, only connects pairs of bright states with equal and opposite energy[10]. As the bright states hybridize in (exponentially small in L) energy shells due to the -perturbation, it is plausible that statistically does not involve these pairs. From this property, we obtain:In contrast, the central spin polarization equals in an infinite temperature eigenstate of the thermalizing system, which vanishes as in a fixed magnetization sector.

Fig. 3 indicates that at intermediate disorder , the residue Z is well-captured by perturbation theory with no noticeable system size dependence. For stronger disorder, we however find that Z slightly decreases with L (this can be inferred from the bottom panel of Fig. 6). This suggests that the (exponentially) long relaxation times are associated with slow dynamics of the bath spins adjusting to long time (likely non-thermal) configurations. Likewise, the exponential increase of the AGP norm in the dark sector is mostly due to their mixing with bright states (cf. Fig. 10), in turn implying that the part of the wave function in Eq. (25) is chaotic, i.e. exponentially sensitive to infinitesimal perturbations.

From the presented data, it is not possible to predict what happens as the system size L increases beyond . We propose two possible distinct possibilities: One can also imagine other more exotic scenarios, where, for example, the residue Z remains finite at while the bath becomes ergodic, or conversely, but the whole system remains non-ergodic. We do not discuss these further as we see no indications that they could be realized.

KAM-type: the residue Z remains finite, the bath remains non-ergodic, and persists as .

Ergodic: the residue Z ultimately vanishes as and the AGP norm crossovers to the ETH scaling, .

Figure 11 schematically depicts the possible scenarios of (i) KAM and (ii) Ergodic behavior. The figure shows the system size dependence of the (log) relaxation time (upper panel) and the disorder averaged central spin magnetization in the dark manifold (lower panel). At system sizes smaller than the critical size , we have a region where the dynamics of the system are dominated by the integrable lines and the system quickly relaxes to a non-thermal steady state, with at most polynomial dependence of on L. In the chaotic non-ergodic (CNE) regime, eigenstate mixing gives rise to an exponentially increasing relaxation time for the bath (cf. Eq. (24)), while persistent dark states maintain a non-thermal qubit polarization. As we approach the thermodynamic limit , scenario (i) would result in a continuation of the CNE regime (see dotted curves), while (ii) would show a second critical size marking the onset of ergodic dynamics. For , saturates and the system always reaches local thermal equilibrium under its own isolated dynamics.

Figure 11

Schematic of possible scenarios in the thermodynamic limit. Upper panel shows (log) relaxation time vs L, while lower panel shows disordered averaged central spin magnetization in the dark manifold vs L. The integrable regime (left) is controlled by the integrable lines. In the CNE regime (center), the relaxation time in the bath is exponentially long in L, despite the persistence of qubit polarization in dark states. As (right), two scenarios are shown: (i) KAM (dotted curves), where the CNE regime persists at all system sizes and (ii) Ergodic (solid curve), where the CNE regime crosses over to the normal ETH ergodic behavior.

Both possibilities outlined above are very interesting and have nontrivial implications. If the KAM scenario (i) is realized, then there is a true non-ergodic phase in the thermodynamic limit. In this case, both the dark and the bright sectors behave non-ergodically at all system sizes, violating the ETH. They are characterized by bath-spin relaxation times that are exponentially long in L. These violations will necessarily lead to a breakdown of various thermodynamic relations such as the fluctuation-dissipation theorem, which heavily rely on ETH[57]. If the Ergodic scenario (ii) is realized, then the system will eventually relax in a finite time to a thermal steady state with a thermal value of the central-spin magnetization. Even if this scenario is realized, according to our numerical results, this can only happen at extremely large system sizes (cf. the lower panel in Fig. 6). As the relaxation time scales exponentially with L, it could be astronomically large at before saturation. This suggests that the dark states, while not exact eigenstates in scenario (ii), will be extremely stable and long lived.

Our numerical results do not predict which of the two possibilities is realized. Based on the available data the KAM scenario (i) seems to be the most likely, at least for moderate and , as there are no visible deviations from (Fig. 3) and (Fig. 8). The absence of deviations from these scalings is especially remarkable since there are no small parameters in the system, so there is no obvious estimate for the length scale . Nevertheless, a more careful analysis is needed to reach a definite conclusion.

Discussion

Physically, dark states can be realized in several qubit systems with mesoscopic environments. For example, in diamond systems, a nitrogen vacancy (NV) center serves as a qubit and the electronic spins on the surface act as a bath. In a suitable rotating frame, the Hamiltonian is well approximated by the XX central spin model. Furthermore, as the qubit-bath interactions are dipolar, they decay sufficiently rapidly as the distance between the NV and a surface spin grows that only a handful of surface spins can be experimentally accessed[17,18]. Our results imply that such small NV-surface spin systems exhibit dark states that are robust to the presence of moderate anisotropy and disorder.

One potential avenue for applications involves quantum information processing in the manifold of persistent dark states. To initialize the system in the persistent dark state manifold, one can implement dynamic nuclear polarization (DNP) to repeatedly polarize the central spin and transfer its polarization to the bath. DNP works by harnessing the flip-flop interactions () already present in the Hamiltonian H. This transfer can be achieved with several methods, e.g. tuning external fields to resonant Hartmann-Hahn conditions where flip-flop interactions dominate[24,25]. By Eq. (2), dark states are unaffected by flip-flop interactions, and therefore only bright state populations will be continually transferred to dark states and other bright states. Repeating the process, the system tends to a statistical mixture of persistent dark states[41,76]. The low qubit-bath entanglement of the resulting manifold ensures robust qubit states with large decoherence times for high-fidelity quantum computing.

A closely related application of DNP is to fully polarize a mesoscopic bath. It has long been known that DNP protocols populate dark states in which the bath is only partially polarized, preventing complete bath polarization and severely limiting this goal[77]. Methods to overcome these limitations have been proposed[28,40]. Our results extend these limitations to mesoscopic central spin systems with moderate anisotropy and disorder. A promising avenue for future research would be to characterize experimentally relevant integrability-breaking perturbations which destroy mesoscopic persistent dark states.

Our results also extend the class of mesoscopic systems relevant for applications to quantum memory. Dark states have been proposed for the storage and retrieval of qubit states[8,38]. In one scheme, the qubit is initialized in an arbitrary state, which can be expressed as a superposition of bright and dark states. By controlling the external field which does not couple to dark states, the information about the qubit state can be completely transferred to the surrounding bath state and retrieved at a later time[8]. The scheme immediately generalizes to persistent dark states with moderate anisotropy and disorder, opening avenues for quantum memory in new systems.

To conclude, we investigated the robustness of dark states in a family of central spin models with anisotropic and inhomogenous qubit-bath interactions. The model is integrable along three lines in parameter space, two of which exhibit exact product dark eigenstates in which the central spin is unentangled with its environment. At moderate deviations away from these exact lines, we found persistent dark states whose central spin polarization and entanglement entropy are well-described by perturbation theory at numerically accessible system sizes. We furthermore showed that the extensive set of conserved operators at the integrable lines morph into an extensive set of quasi-conserved operators away from the integrable lines. In quench experiments, these quasi-conserved operators result in non-thermal correlations in a long-lived non-thermal state. To address the possibility of chaotic behavior at larger system sizes than numerically accessible, we investigated the scaling behavior of the norm of the generator of adiabatic deformations of eigenstates with system size. Although the scaling predicts the onset of chaos at any non-zero strength of the (integrability-breaking) perturbation, the ETH is not obeyed and the relaxation time of the system diverges exponentially with system size. While these effects may disappear in thermodynamically large systems, we see no evidence for this at the numerically-accessible system sizes.

Supplementary Information.