Integrability and quench dynamics in the spin-1 central spin XX model

Central spin models provide an idealized description of interactions between a central degree of freedom and a mesoscopic environment of surrounding spins. We show that the family of models with a spin-1 at the center and XX interactions of arbitrary strength with surrounding spins is integrable. Specifically, we derive an extensive set of conserved quantities and obtain the exact eigenstates using the Bethe ansatz. As in the homogenous limit, the states divide into two exponentially large classes: bright states, in which the spin-1 is entangled with its surroundings, and dark states, in which it is not. On resonance, the bright states further break up into two classes depending on their weight on states with central spin polarization zero. These classes are probed in quench dynamics wherein they prevent the central spin from reaching thermal equilibrium. In the single spin-flip sector we explicitly construct the bright states and show that the central spin exhibits oscillatory dynamics as a consequence of the semilocalization of these eigenstates. We relate the integrability to the closely related class of integrable Richardson-Gaudin models, and conjecture that the spin-$s$ central spin XX model is integrable for any $s$.


Introduction
Central spin models provide a minimal description for a central degree of freedom interacting with an environment of surrounding spins. They are ubiquitous in physics, and have recently gained increased attention with advances in quantum metrology and sensing [1][2][3][4][5][6][7][8][9][10][11][12]. In such setups the central degree of freedom is typically well controlled and can be used to sense or influence the environment. In solid-state quantum computing platforms, the central degree of freedom could be the spin associated with an electron (hole) in a quantum dot or that associated with a defect center in diamond, while the environment is composed of nuclear spins [13][14][15][16][17][18]. In cavity-QED systems on the other hand, the cavity acts as the central degree of freedom and the many atoms it interacts with form the environment [19][20][21][22][23][24][25].
On the theory side, central spin models have been widely investigated because of their underlying integrability. Integrability guarantees an extensive set of conserved quantities and allows all eigenstates to be exactly obtained using Bethe ansatz techniques, which has led to various studies of the equilibrium and dynamical properties in these models [10,[26][27][28][29][30][31][32]. For XXX interactions (S x 0 S x j + S y 0 S y j + S z 0 S z j ) between the central spin (at site 0) and the environment spins (at sites j), the central spin model belongs to the class of Richardson-Gaudin models [33][34][35][36]. Such models are integrable for any value of the central spin and their exact solution has long been established.
In this work, we focus on a system where the central spin interacts with its environment through XX interactions. Such spin-flip terms (S x 0 S x j + S y 0 S y j ) ∝ (S + 0 S − j + S − 0 S + j ) naturally arise from dipolar couplings in nuclear magnetic resonance experiments [37][38][39], in nitrogen vacancy (NV) centers [40], and in certain quantum dots [41]. Some of the authors recently showed that the XX model with a central spin-1/2 particle is integrable with two classes of eigenstates: dark states, in which the central spin is maximally polarized along the z-axis and is unentangled with the environment, and bright states, in which the central spin is entangled with the environment [31]. Subsequent work [42] showed that the spin-1/2 XX model remains integrable in the presence of an arbitrarily oriented magnetic field with emergent dark states (building on results in Refs. [31,43,44]). However, in all cases the exact solution strongly depends on the central spin being a spin-1/2 particle.
We here consider the case where a central spin-1 particle interacts with its environment through XX interactions. By explicitly constructing an extensive set of conserved charges and exact Bethe eigenstates, we establish the integrability of the spin-1 model. While the eigenstate structure is different from that of the spin-1/2 model, we can again identify different classes of bright and dark states. The bright and dark states have striking consequences on the dynamics of the central spin and prevent the central spin from equilibrating with its environment. This work is structured as follows. In Sec. 2 we present an overview of the integrability of the spin-1 central spin model, detailing its conserved charges and eigenstates. Its integrability can be closely connected to the integrability of XXZ Richardson-Gaudin models, and in Sec. 3 we review some relevant results. These are then used in Sec. 4 to construct the conserved charges and discuss several simple limits. The exact eigenstates are constructed in Sec. 5. These eigenstates are similar to the eigenstates of the homogeneous model (where all couplings are set to be equal), which can be solved in terms of collective spin operators [45]. We therefore present the eigenspectrum for the homogeneous model before moving on to the inhomogeneous model. Following this theoretical analysis, we probe the (semi)localization properties of the eigenstates and the effect on quench dynamics in Sec. 6 in the limiting case of a single spin-flip excitation above the polarized ground state. Dynamics for quenches to resonance from a maximally mixed and unpolarized environment is presented in Sec. 7. We combine the known structure of the eigenstates and the exact solution at the homogeneous point to make predictions for the long-time values of central spin polarization and show that they retain memory of the initial state. Sec. 8 is reserved for conclusions.
While our construction now explicitly depends on the central spin being a spin-1 particle rather than spin-1/2, the integrability of both models suggests that the central spin model with XX interactions is integrable for arbitrary spin values s 0 . We present three pieces of evidence in support of this conclusion in Appendix C. The first is a numerical calculation of the level spacing ratio distribution in the spin-3/2 model, exhibiting the Poissonian statistics expected in integrable models. The second is the integrability of the effective Hamiltonian in the limit of a large z-field on the central spin. Specifically, up to second order in the strength of the XX interactions, the effective Hamiltonian obtained by a Schrieffer-Wolff transformation is integrable for any value of the central spin. Third, numerical investigations of the corresponding classical (large s 0 ) model-which can be simulated efficiently-show features of integrability. Integrability of the classical model would imply integrability at smaller s 0 within a truncated Wigner approximation [46,47].
Despite this evidence, proving integrability beyond the spin-1/2 and spin-1 cases remains an outstanding challenge. Establishing integrability in these models would be particularly interesting since, apart from the classical central spin model, another particular limit of this family is given by the inhomogeneous Tavis-Cummings model [19,48], which is prevalent in cavity-and circuit-QED.

Overview of main results
The focus of this work is the central spin Hamiltonian describing a central spin-1 particle interacting with an environment of L surrounding spins through an inhomogeneous XX interaction. Both the interaction strengths g j and the spin quantum numbers s j of the surrounding environment particles can be chosen freely. The central spin and the environment spins are subject to external fields along the z-direction with strength ω ′ 0 and Ω respectively. However, since the total z magnetization is conserved, only the detuning ω 0 = ω ′ 0 − Ω governs the structure of the eigenstates. Indeed, a rotating frame transformation by e −iΩS z tot t takes ω ′ 0 → ω 0 and Ω → 0, while leaving the eigenstates (which may be chosen to be eigenstates of S z tot ) unaltered. Without loss of generality, we work within this rotating frame, where the Hamiltonian takes the form For convenience, we introduce the environment spin effective raising/lowering operators This model is illustrated in Fig. 1.
We establish the integrability of the Hamiltonian (1) by constructing both an extensive set of conserved charges and the exact eigenstates. For each environment spin S j there is an associated conserved charge given bỹ in which P 0 = 1 − (S z 0 ) 2 is a projector on central spin |0⟩ 0 , {·, ·} the anticommutator, and These conserved charges mutually commute and commute with the central spin Hamiltonian. Note that these charges consist of up to 4-body operators, whereas the conserved charges of the spin-1/2 central spin Hamiltonian consist of up to 2-body operators [31].
Two different classes of exact eigenstates with a fixed number of spin excitations can be constructed by adding excitations to the vacuum state, which is defined as the fully polarized state |∅⟩ = ⊗ L j=1 |−s j ⟩ if the environment spin at site j has total spin s j . Environment states are then given by unnormalized Bethe states parametrized by N (possibly complex) variables v 1 , v 2 , . . . , v N . These variables are also known as rapidities 1 . The number of excitations is bounded by N ≤ 2 L j=1 s j , where equality corresponds to the fully polarized state ⊗ L j=1 |s j ⟩. First, the Hamiltonian has a class of degenerate dark eigenstates, where the central spin is maximally polarized along either the positive or negative z-direction. For M < 0, all dark states have central spin down, reading with rapidities satisfying a set of Bethe equations such that G − |v 1 , . . . , v N ⟩ = 0. Note that the rapidities, and thus the dark states |D⟩, are independent of the magnetic field strength ω 0 . The dark states span a degenerate manifold of energy E = −ω 0 : Above half-filling, the dark states have central spin polarization |+1⟩ 0 . The environment states are annihilated by G + and can be similarly obtained by spin inversion.
The second class of eigenstates are bright states, in which the central spin is entangled with the environment states. Bright states can be parametrized as satisfying the eigenvalue equation provided the rapidities satisfy the set of Bethe equations These bright states contain N + 1 spin excitations on top of the vacuum state |−1⟩ 0 ⊗ |∅⟩ and have total spin magnetization M = N − L j=1 s j . While completeness of the Bethe ansatz is typically not easy to establish, in Sec. 5 we argue that these bright and dark states exhaust all possible eigenstates, such that the Bethe ansatz is complete for this model.
Beyond these results on integrability, we also characterize the dynamical properties of the central spin-1 model in several relevant limits. In the single-excitation sector a more detailed analysis of the eigenstates is possible, even in the L → ∞ limit. The model continues to support both dark and bright states, but the spin-flip excitation is neither localized nor delocalized, but rather semilocalized [49] -the central spin has an O(1) probability of carrying the excitation, to be contrasted with the O(1/L) probability for the environment spins. This feature can be directly observed in quench dynamics, where it implies that initial states where the central spin carries the excitation exhibit a nonvanishing oscillation in ⟨S z 0 (t)⟩, even in the L → ∞ limit. Consequences of integrability remain visible with an unpolarized environment. We show numerically that the remanent central spin magnetization in a quench to resonance (ω 0 = 0) differs from the Gibbs ensemble prediction. While the remanent magnetization is determined by the diagonal ensemble corresponding to the full Hamiltonian, we show that its value is well-approximated by a homogeneous dephasing approximation (HDA). Within this approximation we assume that the matrix elements of S z 0 in the inhomogeneous model can be accurately approximated by the matrix elements in the homogeneous model, while keeping the energies different. The latter leads to dephasing at long times, which is absent in the homogeneous model, and an approach to the diagonal ensemble prediction in the inhomogeneous model.

Factorizable Richardson-Gaudin Hamiltonians
In this section, we review various properties of the class of factorizable Richardson-Gaudin Hamiltonians [33][34][35][36] that will be useful in establishing the conserved charges and eigenstates of the spin-1 central spin Hamiltonian. The integrability and eigenstates of the spin-1/2 model were similarly obtained using the properties of these models in Ref. [31], but the construction for the spin-1 model is more involved and cannot be seen as a direct generalization of the spin-1/2 model. The family of factorizable Hamiltonians can be written as The Hamiltonian H(α) is integrable for every choice of α, and results for its eigenstates and conserved charges can be found in, for example, Refs. [36,50]. For reference, the conserved quantities are: These satisfy [H(α), Q j (α)] = [Q j (α), Q k (α)] = 0, for all j, k = 1, . . . , L. In the conserved charges of the central spin model we have Q j ≡ Q j (α = 0). Note that different equivalent expressions for these conserved charges appear in the literature: the asymmetric expressions used here [50][51][52][53][54][55][56] as well as more symmetric ones [34,35,[57][58][59].
All eigenstates can be written as Bethe states, where a Bethe state with N spin excitations on top of the vacuum state |∅⟩ = ⊗ L j=1 |−s j ⟩ is defined as expressed in terms of generalized spin raising operators that depend on the parameters v 1 , v 2 , . . . v N . These Bethe states are eigenstates provided these rapidities satisfy a set of Bethe equations resulting in eigenvalue equations for the conserved charges, as well as the Hamiltonian, Note that the eigenstates and eigenvalues have an implicit dependence on α through the Bethe equations (17). As apparent from Eq. (19), the rapidities can be given an interpretation as spin excitation energies on top of a vacuum energy, with the Bethe equations acting as a set of self-consistency equations. The model exhibits a quantum phase transition at |α| = 1. Consider, for example, α = 1, at which point the Hamiltonian reduces to a positive semi-definite Hamiltonian G + G − . The ground states have energy zero, are necessarily annihilated by G − , and are highly degenerate (see, for example, Ref. [56]). The Bethe ground states are parametrized by N rapidities satisfying All excited states have strictly positive energy and correspond to Bethe states with a single diverging rapidity v → ∞, and the remaining N − 1 (finite) rapidities satisfy the set of Bethe equations leading to a strictly positive (energy) eigenvalue L j=1 2s j g 2 j − 2 N −1 a=1 v a for the Hamiltonian G + G − . The ground and excited states result in dark and bright states respectively in Sec. 5.

Conserved charges
The general conserved charges of the central spin Hamiltonian are most easily derived in a 3 × 3 block-matrix representation, in which the Hamiltonian (1) is given by The different blocks correspond to different eigenvalues of the central spin polarization S z 0 , here ordered as {+1, 0, −1}, and every matrix element acts on the L environment spins. The diagonal terms correspond to ω 0 S z 0 and are proportional to the identity within each block. The off-diagonal factors of √ 2 arise from the action of S ± 0 connecting different blocks.
The L corresponding conserved charges (5) establishing integrability can be expressed in block-matrix form as These satisfy [Q j , H] = [Q j ,Q k ] = 0, for all j, k = 1, . . . , L. These properties are checked by direct calculation below using properties of the Q j as defined in Eq. (6). The different terms in Eq. (23) can first be motivated by considering two simplifying limits. Far away from resonance. Close to the limit ω 0 → ∞ we can perform a Schrieffer-Wolff transformation [60,61] to obtain an effective Hamiltonian which has block-matrix representation All diagonal elements correspond to Richardson-Gaudin integrable Hamiltonians for the environment from Sec. 3, such that the effective Hamiltonian itself is also integrable. The diagonal elements of Eq. (23) are the dominant terms for ω 0 → ∞ and correspond exactly to the conserved charges of the Hamiltonian in Eq. (25). At resonance. In the opposite limit where ω 0 = 0, i.e. at resonance, the diagonal elements of both the Hamiltonian (22) and the conserved charges (23) vanish. In this limit the commutator of the Hamiltonian with the conserved charges can be directly evaluated as This expression for the commutator is independent of the choice of Q j . The single nontrivial element now vanishes since G + G − + G − G + is again a Richardson-Gaudin integrable Hamiltonian, with conserved charges Q j . An alternative way of obtaining the conserved charge in this limit is by noting that H 2 commutes with (S z 0 ) 2 . If we consider the component of H 2 in the space with zero spin polarization, we have that returning the factorizable Hamiltonian with conserved charges Q j for the environment space. The Hamiltonian squared then has conserved charges P 0 Q j , such that at resonance the Hamiltonian itself has conserved charges {H, P 0 Q j }. Expressing these charges as a block matrix then returns the conserved charges (23) with ω 0 = 0. General. The general conserved charges (23) are linear in ω 0 , such that they interpolate between the two limiting cases. The commutation relation at arbitrary values of ω 0 can be checked by evaluating the commutator [H,Q j ], which reads where we have introduced the shorthand Q No properties of the Q j have been used yet. The diagonal element again vanishes since Q j are the conserved charges (15)]. The off-diagonal elements can be shown to vanish by noting that G + Q It is an open question how the integrability of this model and the construction of the conserved charges can be incorporated in the general framework of (Richardson-Gaudin) integrability such as generalized Gaudin algebras [58], constructions based on solutions to the Yang-Baxter equation such as the algebraic Bethe ansatz [62,63], or based on solutions to the "generalized" classical Yang-Baxter equation [44,64] and corresponding modified Bethe ansatz [65].

Eigenstates
Central spin XX models generally support two different classes of eigenstates: bright and dark states. All such states can be obtained explicitly and systematically. We note that the construction of dark states does not depend on the value of the central spin (see, for example, Refs. [45,[66][67][68][69]), such that the construction for dark states in the spin-1/2 model immediately extends to the current model. The construction of the bright states, however, is particular to this model, and these exhibit a richer behavior as compared to the spin-1/2 case.

Homogeneous limit
It is instructive to first examine the special case of homogeneous couplings, g j = g for all j. In this case G ± is proportional to the total spin raising/lowering operator on the environment and we can write The eigenstates of the model (29) can be found without resorting to Bethe ansatz machinery. The homogeneous limit has symmetries where J 2 = 1 2 (J + J − + J − J + ) + (J z ) 2 is the total spin operator for the environment. The central spin Hamiltonian can be represented as a block-diagonal matrix in a fixed (M, J) sector of linear dimension 1, 2 or 3, depending on the relation between J and M . All such matrices can be explicitly diagonalized to return the spectrum of the homogeneous central spin model.
Bright states. We first consider the case J > |M |. The Hamiltonian has contributions from 3 × 3 blocks spanned by where |m 0 ⟩ 0 denotes the eigenstates of S z 0 and |J, M J ⟩ is a simultaneous eigenstate of J 2 and J z with eigenvalues J(J + 1) and M J , respectively. In this block H takes the form The resulting states generally depend strongly on ω 0 and are known as bright states-to be contrasted with the dark states that will be introduced later in this section.
In the special case of resonance (ω 0 = 0) the eigenstates can be analytically constructed and used to further divide the bright states in two subclasses. The Hamiltonian matrix reduces to A single eigenstate can be constructed as with |B| 2 = (c + JM ) 2 + (c − JM ) 2 and energy E 0 = 0. A notable feature of this class is that these states have no weight on |0⟩ 0 . This feature will persist to the inhomogeneous case, and we will refer to these states as double states.
The remaining two eigenstates follow as with degenerate energy eigenvalues E ± = √ 2gB. The persistent feature of this class is that exactly half the weight of the state is on |0⟩ 0 . These states form pairs and since they always have support on all three central spin states we will refer to these as triple states.
Next, we can consider the case J = |M | with M ̸ = 0. The Hamiltonian now only couples two different states, depending on the sign of M , Constructing the central spin Hamiltonian in this basis leads to a 2×2 matrix. For example, for M > 0, which can be explicitly diagonalized to return a pair of eigenvalues These states have common properties of both double and triple states. At resonance the corresponding eigenstates are supported on two states, as with double states, but half the weight of the state is on |0⟩ 0 , as with the triple states. Since we will be interested in quench dynamics of the central spin magnetization, we also refer to these as triple states. Dark states. If we consider a total environment spin where J = |M | − 1, the blocks reduce to 1 × 1 blocks. This condition enforces that the environment is in a state |J, J⟩ or |J, −J⟩, and the corresponding states can now only take the form Crucially, these states have the property that they are annihilated by both S + 0 J − and S − 0 J + , the interaction terms in the Hamiltonian, such that |D⟩ = |±1⟩ 0 ⊗ |J, ±J⟩ is an eigenstate of H with eigenvalue ±ω 0 . These product eigenstates |D⟩ are called dark states and are independent of ω 0 . Note that dark states with central spin state |+1⟩ 0 only appear for M > 0, whereas dark states with central spin state |−1⟩ 0 only appear for M < 0. In the specific case where J = M = 0, the homogeneous model supports additional dark states of the form |0⟩ 0 ⊗ |0, 0⟩, which are eigenstates of the Hamiltonian with zero eigenvalue.
In Appendix A, we provide the counting of each class of states, assuming each environmental spin is spin-1/2. In the limit of large L, there are twice as many triple states as double states, while the ratio of the number of dark states to that of triple states scales as |M |/L. We use these results to predict the late-time expectation value of central spin projectors in Sec. 7.

The inhomogeneous model
The eigenstates of the inhomogeneous model can be constructed as Bethe states that are a direct generalization of the eigenstates of the homogeneous model.

Dark states
At every value of the magnetic field ω 0 , the central spin Hamiltonian in Eq. (1) supports a set of dark eigenstates in which the central spin is not entangled with the environment spins. These dark states |D⟩ are product states of the form |−1⟩ 0 ⊗|D − ⟩ or |+1⟩ 0 ⊗|D + ⟩, where the environment state is adiabatically connected to the states |J, ±J⟩ and satisfies [31,45,[66][67][68] G ± |D ± ⟩ = 0.
This condition guarantees that such dark states are annihilated by the (inhomogeneous) interaction part of the Hamiltonian (1), as well as being eigenstates of the central spin term S z 0 with eigenvalues ±1. As such, dark states are independent of the central spin field ω 0 and form degenerate manifolds with energy ±ω 0 .
As outlined in Ref. [31], the environment states correspond to the ground states of the factorizable Hamiltonians G + G − and G − G + , which can be expressed as Bethe states satisfying the Bethe equations (20). As reviewed in Sec. 3, for M < 0 these dark states can be written as with rapidities satisfying the Bethe equations

Bright states
The bright eigenstates can similarly be related to the eigenstates of the factorizable Richardson-Gaudin models, albeit in a more involved way. Specifically, we consider an ansatz expressed in terms of a single environment state |ψ⟩ and a free parameter κ, writing The above state is an eigenstate of the Hamiltonian (1) with eigenvalue κ provided the environment state |ψ⟩ satisfies the (self-consistent) eigenvalue equation This equation is self-consistent because the Hamiltonian H κ depends on the eigenvalue, but, crucially, is Richardson-Gaudin integrable for every choice of κ. As such, its eigenstates can be exactly constructed as Bethe states for every choice of κ. The Hamiltonian is (up to a prefactor) the Hamiltonian from Eq. (14) from Sec. 3, where the parameter α can be determined as ω 0 /κ. In order to find a set of Bethe equations we can express the eigenvalue κ in terms of rapidities, and now the Bethe equations for these rapidities need to be modified to take into account the self-consistency. This approach directly returns the equations (12).

Spectrum at resonance
The eigenstates of the inhomogeneous model at resonance can again be compared with the eigenstates at resonance in the homogeneous limit, recovering the double, triple and dark states. At resonance, ω 0 = 0, the self-consistency equation for the bright states (46) reduces to a regular eigenvalue equation. In this limit the self-consistent equation can be rewritten as such that the environment states will correspond to the eigenstates of the above (integrable) Hamiltonian. Since the Hamiltonian G + G − + G − G + is positive definite, κ 2 is always positive. For a given eigenstate of this Hamiltonian with eigenvalue κ 2 /2, the central spin Hamiltonian has two corresponding eigenstates with eigenvalue ±κ, given by These are the (normalized) triple states identified previously in the homogeneous limit, and continue to have exactly half their weight on |0⟩ 0 . The double states, with zero energy and vanishing weight on |0⟩ 0 , can be constructed in an alternative way (since in this limit the corresponding Bethe equations become singular): where the inverse should be interpreted as a pseudo-inverse, with the condition that the pseudo-inverse should act as the actual inverse on the environment states |ψ⟩, i.e.
This condition can be satisfied if we consider an initial state that has vanishing overlap with the dark states, since the dark states lie in the kernel of either G + G − or G − G + . For example, if we consider M < 0, all dark states are annihilated by G + G − whereas G − G + has no dark states, such that the inverse of G − G + is well defined. Taking the state |ψ⟩ to be an excited state of G + G − , every excited state gives rise to a well defined double state.

Completeness
It is possible to count the total number of dark and bright states and show that they exhaust all eigenstates of the central spin Hamiltonian (1). The number of these states depends on the choice of environment spins, and for concreteness we here focus on the case where each environmental spin is spin-1/2 and M < 0. The argument for completeness does not depend on the specific choice of spins or total magnetization.
The total number of dark states is set by the number of solutions to For a total magnetization M and a dark state |−1⟩ 0 ⊗ |D − ⟩, the state |D − ⟩ has a magnetization M + 1 and the state G − |D − ⟩ has magnetization M . The total number of solutions to the above equation is given by the dimension of the former Hilbert space (fixing the number of variables) minus the dimension of the latter (fixing the number of constraints), and we find 2 The same result can be recovered in the homogeneous case (see Eq. (84)).
For the bright states, the total number of solutions to the self-consistent equation can be found by plotting the spectrum of the Hamiltonian as a function of κ, as illustrated in Fig. 3 for generic choices of the interaction strengths and ω 0 . Our conclusions do not depend on any specific choice of the parameters. Any intersection between this spectrum and the dashed red line denoting κ determines a solution to the self-consistent equation (46). The number of solutions can now be directly related to the number of bright states: the different lines in Fig. 3 correspond to the different eigenstates of the Hamiltonian (53). As κ → ±∞ all eigenvalues go to zero. At intermediate values of κ all eigenvalues are monotonously decreasing, as follows from the Hellmann-Feynman theorem: Here we have made use of the positive semi-definiteness of G ± G ∓ . Since the Hamiltonian diverges at κ = ±ω 0 , there are now two options for every eigenvalue E(κ): either these diverge at κ = ±ω 0 and the eigenvalue has two vertical asymptotes, or the corresponding eigenstate is annihilated by the residue of H κ at κ = ω 0 or κ = −ω 0 and the eigenvalue has a single vertical asymptote (the state cannot be annihilated by both, as will be made apparent shortly). In the former case the eigenvalue has 3 intercepts with the diagonal line, leading to 3 bright state solutions per environment state. In the latter case the eigenvalue has 2 intercepts with the diagonal line, leading to 2 bright states per environment state. Crucially, the number of states that are annihilated by the residue is exactly equal to the number of dark states, since these are the states that are annihilated by G + G − (or G − G + for M > 0). As such, the total number of bright state solutions equals three times the environment space dimension minus the number of dark states, which combined with the total number of dark states returns the full dimension of the spin-

Single excitation
In the case of a single spin excitation the eigenstates are amenable to a more detailed analytical treatment, even away from resonance and in the limit of an infinite environment L → ∞. In the following, we first analyze the localization properties of the (bright) eigenstates and then present exact results for quench dynamics starting from a product state.

Multifractality and semilocalization
In the case of a single excitation, the Hamiltonian (1) can be written as a so-called arrowhead matrix. Such models generally support both dark and bright eigenstates, and these states have recently gained attention [49] in the context of semilocalization, being neither fully localized nor fully delocalized. Calculations of the inverse participation ratio (IPR) instead indicated a multifractal behavior.
The IPR is defined as with |ψ j | 2 the component of the (normalized) wave function where the excitation is located on spin j. For a delocalized eigenstate, all components are on the order 1/L, resulting in an IPR scaling with L as P(q) = O L 1−q , whereas a localized eigenstate has a few components O(1), resulting in a scaling of P(q) = O(1). A change of scaling as q is varied is a signature of multifractality in the eigenstate [49,[72][73][74]. For a single excitation in the central spin model, there are L − 1 dark states and 2 bright states. In the construction of the bright states (45), the environment state |ψ⟩ is necessarily the vacuum state |∅⟩. The wave function reads with κ satisfying As G − |∅⟩ = 0 and G − G + |∅⟩ = [G − , G + ] |∅⟩ = 2 L j=1 s j g 2 j |∅⟩, the self-consistent eigenvalue equation simplifies to a quadratic equation for κ. This quadratic equation can be explicitly solved to return the two bright states.
In order to have a finite κ value when the number of environment sites L goes to infinity, we will consider a distribution of interaction strengths g j =g j / √ L withg j distributed in some fixed interval. The quadratic equation returns two solutions corresponding to two bright states with Crucially, κ stays finite in the limit L → ∞, resulting in both a finite eigenvalue and nonzero components in the wave function (56). The normalized components of the wave function immediately follow as The IPR can be calculated from these components as Assuming a uniform distribution for 2s jg 2 j in a finite interval [0, 2g 2 ], the sum can be explicitly evaluated to return The scaling of the IPR with L for a fixed q results in This quantifies what is already apparent from the parametrizations (56) and (59)

Quench dynamics
The effect of semilocalization can be directly observed in quench dynamics. We consider quenches where the system is initially prepared in a product state, with the single excitation localized either on the central spin or on one of the environment spins, and is subsequently evolved using the central spin Hamiltonian.
For simplicity, we focus on the dynamics of the central spin magnetization ⟨S z 0 (t)⟩ starting from a general initial state |ψ 0 ⟩. Since the dark states are eigenstates of S z 0 all nontrivial dynamics is due to the two bright states, and we can write for an excitation initially localized on site j. Parameters: L = 12, ω 0 = 2, andg j is uniformly distributed in the interval [1,2] for environment spins with s j = 1/2.
where we have labeled the two bright states by their eigenvalues κ ± = −ω 0 /2± ω 2 0 /4 + 2g 2 and the (L − 1) dark states as D. The central spin magnetization oscillates with a single frequency κ + − κ − = 2 ω 2 0 /4 + 2g 2 , and both the amplitude of the oscillations and their average value are determined by the overlaps with the bright states.
Consider first the case where the initial state consists of an excitation on the central spin, i.e. |ψ 0 ⟩ = |0⟩ 0 ⊗ |∅⟩. This state has a vanishing overlap with the dark states, and the contribution from the bright states can be calculated using the explicit parametrization (56) as which holds for both bright states |κ ± ⟩. The resulting central spin dynamics immediately follows as This result is illustrated in Fig. 4 and is identical to the dynamics in the homogeneous model (29) with interaction strength g = g/ √ 2J and environment spin J = L j=1 s j . Crucially, the amplitude of the central spin oscillation remains finite in the limit L → ∞ provided g 2 remains finite. The nonvanishing amplitude of the oscillation is a direct consequence of the semilocalized nature of the bright states: the overlap between the initial state and the two bright states remains O(1) in this limit.
This finite amplitude oscillation can be contrasted with the central spin dynamics for an initial product state localized on an environment spin. The amplitude of the central spin oscillations is set by ⟨ψ 0 |κ − ⟩⟨κ + |ψ 0 ⟩ and hence by the component ψ j from Eq. (59) for an initial excitation localized on spin j. The individual overlaps scale as O 1/ √ L , such that the total amplitude of the spin oscillations will scale as O(1/L). As illustrated in Fig. 4(a) and Fig. 5(a), central spin oscillations are indeed suppressed for states initially localized in the environment. A similar behavior is observed in the dynamics of the environment spin polarizations ⟨S z j (t)⟩ for an excitation initially localized on spin i. Since the dark states are no longer eigenstates of the observable these will also contribute to the dynamics, and the spins will oscillate with three frequencies The dynamics is illustrated in Fig. 4(b). Note that the frequencies are generally not commensurate, but the spin dynamics may exhibit approximate revivals whenever the two frequencies κ ± + ω are close to being commensurate (in which case the third frequency κ + − κ − = (κ + + ω) − (κ − + ω) is also close to being commensurate), as apparent in Fig. 4(b). The revivals become exact in the special case of quenches to resonance (ω 0 = 0), in which case the three frequencies reduce to two commensurate frequencies, 2g 2 and 2 2g 2 . In this scenario the environment spins oscillate periodically with a frequency that is half the oscillation frequency of the central spin, as illustrated in Fig. 5(b). In all cases, the amplitude of these oscillations scales as O(1/L). This scaling can be understood by noting that L j=1 ⟨ψ 0 |S z j (t)|ψ 0 ⟩ = 1 − ⟨ψ 0 |S z 0 (t)|ψ 0 ⟩, with the right-hand side being O(1). Because of the delocalization in the environment all contributions to the summation are on the same order, leading to the observed O(1/L) scaling.
To summarize, we note that (semi)localization of the eigenstates can be observed through measurements of the central spin polarization. In an initial state that is supported on the localized component of the bright state (i.e. on the central spin) ⟨S z 0 ⟩ will oscillate with an amplitude that does not vanish as the system size goes to infinity. This highly nonthermal behavior can be contrasted with the behavior of initial states that are supported on the delocalized environment spins, which exhibit oscillations that vanish as L → ∞.

Quenches to resonance with an unpolarized environment
In this section we consider generic quenches to resonance and use the known structure of the eigenstates to show that central spin observables do not relax to thermal equilibrium. Specifically, we consider with the latter two corresponding to projectors on central spin states |0⟩ 0 and |±1⟩ 0 respectively. While exact predictions in the inhomogeneous model are currently out of reach, we show that, for not too strong inhomogeneities, the late-time expectation values in the inhomogeneous model are well approximated by the diagonal ensemble expectation values in the homogeneous model. We refer to this approximation as the homogeneous dephasing approximation (HDA). Inhomogeneity in the couplings breaks the degeneracy of states in the homogeneous model, causing dephasing between formerly degenerate eigenstates. Meanwhile, the matrix elements of central spin observables between eigenstates are not significantly affected. A similar approximation was used in Ref. [75] for a nonintegrable Ising model in a many-particle dephasing regime. The HDA ignores any change to said matrix elements, and only accounts for dephasing.
Depending on the initial state and measured observable we can systematically probe the effect of bright states, both triple and double, as well as dark states. Specifically, we consider an initial state where the central spin is polarized in the state |m 0 ⟩ 0 , the environment (which we again take to consist of spin-1/2 particles) is at infinite temperature in a fixed magnetization sector M E = M − m 0 , and the total magnetization 3 M = 1. The initial density matrix can be written as Here 1 M E acts as the identity on states with magnetization M E and as zero everywhere else. For convenience we take L, the number of environment spins, to be even.

Prediction of late-time values under the HDA
The HDA simplifies the prediction for late-time values of central spin observables by circumventing the use of the Bethe Ansatz solution, which is mathematically cumbersome. Within the diagonal ensemble, the late-time values of all projectors are determined by their overlaps with the eigenstates in the inhomogeneous model. Under the HDA, these overlaps are approximated by those with the corresponding states in the homogeneous model. This approximation can be justified by numerically comparing the expectation values of projectors in the homogeneous model with those of the inhomogeneous model. Fig. 6 shows that the eigenstate expectation values of P 0 , P 1 , S z 0 and J 2 in the inhomogeneous model have small spread around the homogeneous limit.
Initial state m 0 = 0. We first consider the case where m 0 = 0 and hence M E = M . The initial density matrix is only nonvanishing in the triple bright state manifold (35). The diagonal values are equal ⟨B|ρ(t = 0)|B⟩ = 1/2, as the triple states have half their weight on states with central spin value 0. The diagonal ensemble (obtained by setting the off-diagonal elements in the triple bright basis to zero) is thus The property ⟨B|P 0 |B⟩ = 1/2 also implies that the long-time expectation value of P 0 is where we further used that the total number of triple bright states is 2Z E . Under the HDA, we find for the projectors on central spin ±1 that (see Appendix B) where is the degeneracy of environment states with spin J, and Eq. (72) follows from Stirling's approximation in the limit L ≫ |M | (far away from the single-excitation limit).
Initial state m 0 = −1. For M > 0 (such that all dark state have central spin |+1⟩ 0 ), the initial state has overlap with both the double and triple bright states, but not with the dark states. In this case, the HDA gives: where the approximation on the second line again holds in the limit L ≫ |M |. The additional exponential factor arises from the environment sector M E = M + 1.
For the projectors on the central spin states |±1⟩ 0 we similarly find and Initial state m 0 = +1. For M > 0, the dark states contribute to the quench dynamics.
Following the same steps as above, the late-time values are, for L ≫ |M |, The expressions above demonstrate that under the HDA, the late-time values of central spin observables retain memory of the initial state m 0 . In contrast, the maximally mixed Gibbs ensemble in a fixed M sector predicts the same late-time values for each central spin projector P 0,±1 , regardless of initial state 4 (up to O(|M |/L) corrections). For instance, in the limit L ≫ |M |, we find that Tr [ρ(t → ∞)P 0 ] approaches 1/4 for m 0 = ±1 and 1/2 for m 0 = 0, clearly differing from the Gibbs prediction of 1/3.

Comparison with the inhomogeneous model
We can now compare the theoretical predictions in Sec. 7.1 with numerical results for the inhomogeneous model. In the case where m 0 = 0, the HDA prediction (70) applies exactly. This is because the initial density matrix is only non-vanishing in the triple state manifold (48), which has exactly the same weight on |0⟩ 0 and counting as the triple states in the homogeneous model. Otherwise, the late-time values in the diagonal ensemble are different between the inhomogeneous and homogeneous models. Nonetheless, late-time values of central spin observables are well approximated by the HDA. These approximations are compared with numerical results for the inhomogeneous model in Fig. 7. In all cases the diagonal ensemble from the homogeneous model accurately reproduces the steady-state value of the inhomogeneous model. Furthermore, the late-time values for different initial states m 0 are clearly different. Fig. 8(a) shows the corresponding dynamics of ⟨S z 0 (t)⟩. Crucially, in a given total magnetization sector M , the late-time expectation values heavily depend on the initial value of the central spin m 0 and differ from the Gibbs ensemble prediction, which can be calculated as This scaling with L can be contrasted with the numerically observed scaling of ⟨S z 0 ⟩ HDA , as illustrated in Fig. 8 (b). While the m 0 = +1 curve shows the 1/L scaling from the Gibbs predictions, the m 0 = 0 and m 0 = −1 curves show different scaling exponents, approximately given by L −0.8 and L −0.7 respectively.
In the case where the initial environment does not have a fixed magnetization and is at infinite temperature, i.e. ρ E ∝ 1, one can also perform similar calculations for all M sectors and perform a weighted average. The resulting long-time values for the polarization follow as ⟨S z 0 ⟩ HDA = O 1/ √ L for L → ∞ as illustrated in Fig. 9(d). We note that this result is consistent with numerical results in the classical model (Fig. 11 in Appendix C.2), and inconsistent with the Gibbs prediction (82). These results show that even in highly excited states, the integrability of the inhomogeneous model can be detected by the remnant memory of the initial state m 0 in late-time central spin observables.

Conclusion and discussion
We have established the integrability of the spin-1 central spin XX model by providing the exact construction of Bethe eigenstates and the extensive set of conserved charges, extending the results in Ref. [31] for the spin-1/2 central spin XX model. Like the spin-1/2 model, the eigenstates in the spin-1 model can be broadly classified into dark states and bright states, with the bright states showing semilocalization [49,74] in the single-excitation sector. Unlike the spin-1/2 model, bright states in the spin-1 model on resonance can further be classified into double or triple bright states.
The eigenstate structure of the spin-1 model prevents the central spin from reaching thermal equilibrium in quenches to resonance. In particular, for weakly inhomogeneous couplings, the late-time values of central spin observables approach diagonal ensemble expectation values in the homogeneous models. We expect this can be observed experimentally. Based on our numerical results, the time required to reach these late-time values is on the order of 10τ , where τ = ( L j=1 g 2 j ) −1/2 . This is within the range of the spin relaxation time T 1 in NV systems, which limits the measurement of central spin observables in the z basis. Indeed, the relaxation time we estimate as 10τ is essentially the dephasing time T 2 , which can be much shorter than T 1 . For instance, at room temperature, T 1 has been observed to exceed 1 ms [76], while Ref. [77] recently measured T 2 ≈ 1 µs.
In Appendix C, we have provided three pieces of evidence that support the integrability of the XX central spin model for any value of central spin s 0 . First, numerical calculations of level spacing ratios in the spin-3/2 model show Poisson level statistics, as expected of integrable models. Second, the effective Hamiltonian at large ω 0 for arbitrary s 0 is integrable. Third, numerical simulations of the fully classical model (s 0 , s j → ∞) show residual memory in the late-time central spin magnetization, supporting integrability of the classical equations of motion. Within the truncated Wigner approximation, said classical equations govern dynamics for any value of s 0 [46,47]. An obvious future direction is to rigorously establish integrability at any s 0 .
Other technical questions remain open. While the conserved charges in the spin-1 central spin XX model are closely related those of the XXZ Richardson-Gaudin integrable models, it is unclear how to incorporate them into the general framework of Richardson-Gaudin integrability. Another challenge remains to directly understand in which way the conserved charges constrain, for example, late-time observables in dynamical experiments.

A Counting of states in the resonant homogeneous model
In this section, we provide the counting of the three classes of states (dark, double and triple bright states) in the homogeneous model on resonance.
The degeneracy of every eigenvalue is set by the total number of ways in which the L environment spins can be combined to form a total spin J. We now focus on the case where each environmental spin is spin-1/2. For a given M , J can then take integer values ranging from J = J min = max(0, |M | − 1) to a maximal value of L/2. Each of the spin-J irreducible representations has multiplicities given by entries in Catalan's triangle Dark states reside in the J min sector. Thus, for a fixed total magnetization M ̸ = 0, the degeneracy of the dark states immediately follows as For bright states, the total number of double states is given by in M ̸ = 0 sectors. The triple states are allowed in J ≥ |M | for |M | ̸ = 0, leading to It is easily checked that the total number of dark and (double and triple) bright states leads to the expected number of eigenstates in each magnetization sector.
For a given M ̸ = 0 sector, the ratio of the number of double states to that of the triple states is given by remaining finite in the limit of large L. The ratio of the number of dark states to that of the triple states is given by In all cases, each class of states spans a nonvanishing fraction of the Hilbert space in the thermodynamic limit L → ∞ provided |M | scales with L. Keeping |M | fixed and increasing L, the fraction of dark states vanishes as O(1/L). For M = 0, the number of J = 0 dark states is given by C(L/2, L/2). The number of double states is given by in the M = 0 sector. There are triple states in the J ≥ 1 sectors, resulting in The ratio N double /N triple in this sector is exactly 1/2, while B Expectation values at large system size under the HDA In this section we detail the approximations used in evaluating the summations for the diagonal ensemble expectation values in the homogeneous model. In the limit where J ≪ L, Stirling's approximation gives where Similarly, when M E ≪ L, the same approximation can be used for the ratio where p is a non-negative integer. In all these cases, the summands are dominated by the smallest J, i.e. M 0 . The summation may be approximated by substituting J = M 0 into the expression for the summand and multiply it by an O(1) factor to capture the entire sum. This factor can be extracted by comparing the sum to the summand for specific values of M , M 0 , M E and (large) L. We find that multiplying the summand by a factor of 4 gives the overall best fit to the exact sum. This results in

C Integrability in higher spin models
Motivated by the results from the main text, we conjecture that the central spin Hamiltonian is integrable for any central spin quantum number s 0 . The integrability of this model has already been explicitly shown for s 0 = 1/2 in Ref. [31], and for s 0 = 1 in this work. As a first numerical check, we consider the central spin-3/2 model and calculate the spectral statistics, a widespread numerical check for integrability. For an ordered eigenvalue spectrum {E n } with eigenvalue spacing s n = E n+1 − E n we consider the distribution of the level spacing ratio,r n = min(s n , s n+1 )/max(s n , s n+1 ). This level spacing ratio is expected to obey different universal distributions in integrable and chaotic systems [80]. Due to the presence of dark states there will be a large number of states for which s n = 0, which are here excluded from the statistics. Numerical results are presented in Fig. 10. The level spacing ratio agrees with the Poissonian prediction for integrable systems and clearly differs from the GOE prediction expected for chaotic systems.
We support our conjecture through two additional pieces of evidence: perturbative conserved charges for H in the limit of large ω 0 (Appendix C.1), and numerical signatures of integrability in the classical limit of s 0 , s j → ∞ (Appendix C.2). We also present a semi-classical argument for integrability, assuming exactness of the truncated Wigner approximation.
If the central spin Hamiltonian is integrable for all s 0 , then it is likely that the related inhomogeneous Tavis-Cummings model (with off-diagonal disorder), is also integrable. In Appendix C.3, we expand on this additional conjecture.

C.1 Conserved charges far from resonance
A Schrieffer-Wolff transformation to H in the ω 0 → ∞ limit provides an effective Hamiltonian (24) where H(α) is a factorizable Hamiltonian from Eq. (14) and commutes with S z 0 , and so may be treated as a scalar in each S z 0 sector. Thus, in each sector, H eff is Richardson-Gaudin integrable and has an extensive set of conserved charges given by Q j (α) (15). Combined with the fact that both the s 0 = 1/2 and s 0 = 1 cases are known to be integrable, the existence of this integrable limit is suggestive of integrability for any s 0 .
Higher-order corrections to the effective Hamiltonian may also be computed, though it is unclear if they preserve integrability. We discuss them here for completeness and only note some connections with known integrable models.
The next correction to the effective Hamiltonian (at O g 4 ) is given by Evaluating these corrections in a sector with fixed s 0 and m 0 results in linear combinations of four different operators, where all terms of the form G − G + G + G − and G + G − G − G + cancel out. Each of these operators can again be shown to be Richardson-Gaudin integrable-although this does not guarantee that linear combinations will be integrable. G − G + G − G + = (G − G + ) 2 is the square of a factorizable Hamiltonian, as is G + G − G + G − . The terms G ± G ± G ∓ G ∓ are also each integrable. To see this, we use the following generalization of a relation between Richardson-Gaudin charges introduced below Eq. (28), where Q (±) j = Q j ± S z j and k ≥ 1 is an integer. From this expression and its Hermitian conjugate, we see that a complete set of conserved charges for, say, G + G + G − G − is given by 3Q where j ∈ {1, . . . , L}. Similar charges may be constructed for The two quartic terms that drop out in the effective Hamiltonian are the only terms that are not known to be integrable, such that the fact that they cancel out suggests that H eff does not imply that the higher-s 0 central spin XX models are nonintegrable-for instance, the same perturbative expansion holds for s 0 = 1, which is integrable. Of course, conversely, the existence of an integrable limit in a model does not imply that the model is integrable for all parameters. For example, while the Bose-Hubbard model is believed to be nonintegrable, it does possess integrable limits [81,82]. However, for both the central spin-1/2 and the central spin-1 model the exact conserved charges reduce to the conserved charges of the effective Hamiltonian in the limit ω 0 → ∞, motivating our investigation of the effective Hamiltonian for arbitrary central spin values. As such, the integrability of H (1) eff + H (2) eff appears consistent with the conjecture of integrability and the structure observed in central spin models.

C.2 Classical limit
Taking s 0 , s j → ∞ while keeping g j s 0 s j ∼ g cl j /2 and ω 0 s 0 ∼ ω cl 0 finite results in a model which is formally identical to H, with the spinsS µ 0 andS µ j being classical degrees of freedom. If H is integrable for all values of ω 0 , g j , s 0 and s j , it is natural to suspect that this limit model is also integrable. Conversely, integrability of the classical model suggests special structure with s j < ∞. In this section, we numerically search for nonergodicity (a consequence of integrability) in the classical model (107), which can be simulated efficiently. We note that it is sometimes necessary to include additional corrections to preserve integrability when passing from quantum to classical or vice versa [83,84]. However, the numerics in this section suggests that the classical central spin model does not require such additional corrections in order to result in the nonergodic dynamics expected in integrable models.
Equations  where ϵ µνρ is the Levi-Civita tensor and summation over the index ρ is implied.
If the central spin is initially aligned along the z-axis while the environment is in an infinite temperature state, and a quench to resonance ω cl 0 = 0 is performed, then ergodicity implies a late time average value of the central spin polarization given by Here, is the average ofS z 0 (t) over an ensemble of initial states with a fixed central spin state, and an infinite temperature environment. Fig. 11 demonstrates that the late time average value of ⟨S z 0 (t)⟩ is not 1/(L + 1). Instead, the late time value decreases more slowly with L, as numerically demonstrating nonergodicity. The 1/ √ L scaling of the remanent magnetization is also a feature of the spin-1 model (see Fig. 8). That the same phenomenology persists in the classical limit favors the hypothesis of integrability of the classical model.
Integrability of the quantum model at any values of the central spin, including the classical large s 0 , s j limit, also follows from semiclassical considerations. Namely, for this system one can anticipate that the truncated Wigner approximation (TWA) [46,47] accurately describes dynamics in the large L limit for any values of the central and environment spins s 0 , s j . This feature is general for all large L-models with long range interactions, where classical dynamics governed by Eq. (108) emerges as a saddle point within the path integral formulation of the Heisenberg evolution on a Schwinger-Keldysh contour (see for example Refs. [85][86][87]). Within the TWA the values of the spins are encoded in the Wigner function representing the initial state. Because the spin-1/2 or spin-1 systems are integrable, Eqs. (108)-which are expected to describe dynamics of the central spin in the large L limit-must be non-ergodic as well to avoid thermalization. Because these equations are independent of the spin quantum numbers, we can anticipate that integrability holds for all values of s 0 , s j .
These qualitative considerations cannot be viewed as a proof of integrability, as in general the limits L → ∞ and t → ∞ do not commute, and the TWA is expected to be accurate in the limit L → ∞ first. The opposite limit is much harder to analyze analytically within the TWA and needs further study. Nevertheless, putting these subtleties aside, a combination of analytical and numerical evidence we presented in this paper suggests that the model is integrable in the limit s 0 → ∞ for any L and is integrable in the limit L → ∞ for small values of s 0 = 1/2, 1. Because both of these limits are described by the same semiclassical equations of motion it is natural to assume that the model is integrable for any s 0 .

C.3 Inhomogeneous Tavis-Cummings model
Several interesting models may be obtained as limits of H. Assuming the integrability of H for any values of s 0 , s j , ω 0 , and g j , it is natural to suspect that the limit models are also integrable. One such model was the classical central spin model of Appendix C.2. Here, we remark upon another notable large s 0 limit, which results in an inhomogeneous Tavis-Cummings model.
Taking an alternative large-s 0 limit of s 0 → ∞ with g j √ s 0 ∼ g TC j increases the central spin Hilbert space while maintaining its level spacing. The limit model above the ground state may be expressed as an oscillator model, The Hamiltonian H TC is a generalization of the Tavis-Cummings model (which is known to be integrable [19,48]) with inhomogeneous couplings. Due to its connection with the central spin XX model, we conjecture that the inhomogeneous Tavis-Cummings model is also integrable. The Hamiltonian H TC is known to be integrable when additional non-linear terms are introduced [88], but no r-matrix is known for the model without such additional couplings.
The consequence of integrability in the central spin-1 XX model is that the structure of the homogenous limit persists to large inhomogeneity of the couplings. We speculate that this is also the case in the Tavis-Cummings model-the inhomogeneous model continues to exhibit the phenomenology of the homogeneous model, such as a superradiant transition (as expected from mean-field calculations).