"Bethe-Ansatz-free" eigenstates of spin-1/2 Richardson-Gaudin integrable models

TL;DR

This paper introduces a Bethe-ansatz-free method for constructing eigenstates of the most general spin-1/2 Richardson-Gaudin models, including fully anisotropic XYZ couplings, simplifying analysis without relying on traditional Bethe ansatz techniques.

Contribution

It provides a novel, simpler approach to eigenstate construction that does not require Bethe roots or pseudo-vacuums, applicable to a wide class of integrable models with anisotropic interactions.

Findings

Eigenstates constructed without Bethe ansatz

Method applicable to fully anisotropic XYZ models

Eigenvalues obtained from quadratic equations

Abstract

In this work we construct the eigenstates of the most general spin-1/2 Richardson-Gaudin model integrable in an external magnetic field. This includes the possibility for fully anisotropic XYZ coupling such that the $S^x_iS^x_j$, $S^y_iS^y_j$ and $S^z_iS^z_j$ terms all have distinct coupling strengths. While insuring that integrability is maintained in the presence of an external field excludes the elliptic XYZ model which is only integrable at zero field, this work still covers a wide class of fully anisotropic (XYZ) models associated with non skew-symmetric r-matrices. The eigenstates, as constructed here, do not require any usable Bethe ansatz and therefore: no proper pseudo-vacuum, Bethe roots, or generalised spin raising (Gaudin) operators have to be defined. Indeed, the eigenstates are generically built only through the conserved charges which define the model of interest and the specification of the set of eigenvalues defining the particular eigenstate. Since these eigenvalues are, in general, solutions to a simple set of quadratic equations, the proposed approach is simpler to implement than any Bethe ansatz and, moreover, it remains completely identical independently of the symmetries of the model. Indeed, the construction removes any distinction between XYZ models and XXZ/XXX models and, generically, that between models with or without U(1) so that any difficulties associated with the use of a Bethe ansatz in any of these cases are avoided.

Equations

$${\hat{R}}_i={\vec{B}}_i\cdot {\vec{\sigma }}_i+\sum _{j\ne i}^N\left(X_{ij}{\sigma }_i^x{\sigma }_j^x+Y_{ij}{\sigma }_i^y{\sigma }_j^y+Z_{ij}{\sigma }_i^z{\sigma }_j^z\right).$$

$$B_i^x=\frac{\gamma }{F_x({ϵ}_i)},$$

$$B_i^y=\frac{\lambda }{F_y({ϵ}_i)},$$

$$B_i^z=1,$$

$$F_x({ϵ}_j)$$

$$F_y({ϵ}_j)$$

$${\hat{R}}_i^2=\sum _{j\ne i}^N{\Gamma }_{ij}{\hat{R}}_j+K_i,$$

$${\Gamma }_{ij}=2Z_{ij}.$$

$$r_i^2=\sum _{j\ne i}^N{\Gamma }_{ij}{r}_j+K_i.$$

$$\hat{Q}({\hat{R}}_1\dots {\hat{R}}_N)=\sum _{n_1,n_2\dots n_N=0}^1{B}_{n_1\dots n_N}{\hat{R}}_1^{n_1}{\hat{R}}_2^{n_2}\dots {\hat{R}}_N^{n_N};$$

$$\displaystyle\hat{P}(r_{1}\dots r_{N})\equiv\mathrm{det}\hat{J}\ \ \ \mathrm{with}\left\{\begin{array}[]{l}\hat{J}_{ii}=r_{i}+\hat{R}_{i}\\ \hat{J}_{ij}=-\Gamma_{ij}\ \ \ \forall\ i\neq j.\end{array}\right.$$

$$\displaystyle N(r_{1}\dots r_{N})\equiv\mathrm{det}N\ \ \ \mathrm{with}\left\{\begin{array}[]{l}N_{ii}=2r_{i}\\ N_{ij}=-\Gamma_{ij}\ \ \ \forall\ i\neq j.\end{array}\right.$$

$$|{\psi }_n⟩⟨{\psi }_n|\propto \hat{P}(r_1^n\dots r_N^n)$$

$$\displaystyle=\mathrm{det}\left(\begin{array}[]{ccccc}r_{1}^{n}+\hat{R}_{1}&-\Gamma_{12}&-\Gamma_{13}&\dots&-\Gamma_{1N}\\ -\Gamma_{21}&r_{2}^{n}+\hat{R}_{2}&-\Gamma_{23}&\dots&-\Gamma_{2N}\\ -\Gamma_{31}&-\Gamma_{32}&r_{3}^{n}+\hat{R}_{3}&\dots&-\Gamma_{3N}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ -\Gamma_{N1}&-\Gamma_{N2}&-\Gamma_{N3}&\dots&r_{N}^{n}+\hat{R}_{N}\\ \end{array}\right)$$

$$|{\psi }_n⟩⟨{\psi }_n|=\frac{\hat{P}(r_1^n\dots r_N^n)}{N(r_1^n\dots r_N^n)}.$$

$$|{\psi }_n⟩\propto \hat{P}(r_1^n\dots r_N^n)|\Omega ⟩,$$

$$\hat{P}(r_1\dots r_N)=\sum _{m=1}^{2^N}{C}_m(r_1\dots r_N)|{\psi }_m⟩⟨{\psi }_m|,$$

$$\displaystyle C_{m}(r_{1}\dots r_{N})=\mathrm{det}\left(\begin{array}[]{ccccc}r_{1}+r_{1}^{m}&-\Gamma_{12}&-\Gamma_{13}&\dots&-\Gamma_{1N}\\ -\Gamma_{21}&r_{2}+r_{2}^{m}&-\Gamma_{23}&\dots&-\Gamma_{2N}\\ -\Gamma_{31}&-\Gamma_{32}&r_{3}+r_{3}^{m}&\dots&-\Gamma_{3N}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ -\Gamma_{N1}&-\Gamma_{N2}&-\Gamma_{N3}&\dots&r_{N}+r_{N}^{m}\\ \end{array}\right).$$

$$(r_i^n{)}^2-(r_i^m{)}^2=\sum _{j\ne i}^N{\Gamma }_{ij}\left(r_j^n-r_j^m\right),$$

$$\sum _{j=1}^N{J}_{ij}(r_j^n-r_j^m)=(r_i^n+r_i^m)(r_i^n-r_i^m)-\sum _{j\ne i}^N{\Gamma }_{ij}\left(r_j^n-r_j^m\right)=0.$$

$$|{\psi }_n⟩⟨{\psi }_n|=\frac{\hat{P}(r_1^n\dots r_N^n)}{N(r_1^n\dots r_N^n)}.$$

$${\int }_{-a}^{a}d{r}_1\cdots {\int }_{-a}^{a}d{r}_N\left(\prod _{i=1}^N{r}_i\right)\hat{P}(r_1\dots r_N)={\left(\frac{2a^3}{3}\right)}^{N}1,$$

$$\hat{P}(r_1\dots r_N)\hat{P}(r_1^n\dots r_N^n)$$

$$\displaystyle=\mathrm{det}\left(\begin{array}[]{ccccc}(r_{1}+r_{1}^{n})&-\Gamma_{12}&-\Gamma_{13}&\dots&-\Gamma_{1N}\\ -\Gamma_{21}&(r_{2}+r_{2}^{n})&-\Gamma_{23}&\dots&-\Gamma_{2N}\\ -\Gamma_{31}&-\Gamma_{32}&(r_{3}+r_{3}^{n})&\dots&-\Gamma_{3N}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ -\Gamma_{N1}&-\Gamma_{N2}&-\Gamma_{N3}&\dots&(r_{N}+r_{N}^{n})\\ \end{array}\right)\left|\psi_{n}\right>\left<\psi_{n}\right|.$$

$$\hat{P}(r_1\dots r_N)=\sum _{\sigma \in S_n}\mathrm{sgn}(\sigma )\prod _{i=1}^n{\hat{J}}_{i,{\sigma }_i},$$

$$\hat{P}(r_1\dots r_N)=\sum _{\pi \in \Pi }{K}_{\pi }\prod _{k=1}^r(r_{{\pi }_k}+R_{{\pi }_k})$$

$$\displaystyle K_{\pi}=\mathrm{det}\left(\begin{array}[]{ccccc}0&-\Gamma_{\bar{\pi}_{1}\bar{\pi}_{2}}&-\Gamma_{\bar{\pi}_{1}\bar{\pi}_{3}}&\dots&-\Gamma_{\bar{\pi}_{1}\bar{\pi}_{M}}\\ -\Gamma_{\bar{\pi}_{2}\bar{\pi}_{1}}&0&-\Gamma_{\bar{\pi}_{2}\bar{\pi}_{3}}&\dots&-\Gamma_{\bar{\pi}_{2}\bar{\pi}_{M}}\\ -\Gamma_{\bar{\pi}_{3}\bar{\pi}_{1}}&-\Gamma_{\bar{\pi}_{3}\bar{\pi}_{2}}&0&\dots&-\Gamma_{\bar{\pi}_{3}\bar{\pi}_{M}}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ -\Gamma_{\bar{\pi}_{M}\bar{\pi}_{1}}&-\Gamma_{\bar{\pi}_{M}\bar{\pi}_{2}}&-\Gamma_{\bar{\pi}_{M}\bar{\pi}_{3}}&\dots&0\\ \end{array}\right),$$

$$K_{\pi }=\left(\prod _{k=1}^M-2g{F}_x({ϵ}_{{\bar{\pi }}_k})F_y({ϵ}_{{\bar{\pi }}_k})\right)$$

$$\displaystyle\times\mathrm{det}\left(\begin{array}[]{ccccc}0&\frac{1}{\epsilon_{\bar{\pi}_{1}}-\epsilon_{\bar{\pi}_{2}}}&\frac{1}{\epsilon_{\bar{\pi}_{1}}-\epsilon_{\bar{\pi}_{3}}}&\dots&\frac{1}{\epsilon_{\bar{\pi}_{1}}-\epsilon_{\bar{\pi}_{M}}}\\ \frac{1}{\epsilon_{\bar{\pi}_{2}}-\epsilon_{\bar{\pi}_{1}}}&0&\frac{1}{\epsilon_{\bar{\pi}_{2}}-\epsilon_{\bar{\pi}_{3}}}&\dots&\frac{1}{\epsilon_{\bar{\pi}_{2}}-\epsilon_{\bar{\pi}_{M}}}\\ \frac{1}{\epsilon_{\bar{\pi}_{3}}-\epsilon_{\bar{\pi}_{1}}}&\frac{1}{\epsilon_{\bar{\pi}_{3}}-\epsilon_{\bar{\pi}_{2}}}&0&\dots&\frac{1}{\epsilon_{\bar{\pi}_{3}}-\epsilon_{\bar{\pi}_{M}}}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ \frac{1}{\epsilon_{\bar{\pi}_{M}}-\epsilon_{\bar{\pi}_{1}}}&\frac{1}{\epsilon_{\bar{\pi}_{M}}-\epsilon_{\bar{\pi}_{2}}}&\frac{1}{\epsilon_{\bar{\pi}_{M}}-\epsilon_{\bar{\pi}_{3}}}&\dots&0\\ \end{array}\right),$$

$$\displaystyle\lim_{\Delta\to 0}\mathrm{det}A=\mathrm{det}\left(\begin{array}[]{cccccc}0&-\frac{1}{\Delta}&\frac{1}{\epsilon_{2}-\epsilon_{3}}&\frac{1}{\epsilon_{2}-\epsilon_{4}}&\dots&\frac{1}{\epsilon_{2}-\epsilon_{{M}}}\\ \frac{1}{\Delta}&0&\frac{1}{\epsilon_{2}-\epsilon_{3}}&\frac{1}{\epsilon_{2}-\epsilon_{4}}&\dots&\frac{1}{\epsilon_{2}-\epsilon_{{M}}}\\ \frac{1}{\epsilon_{3}-\epsilon_{2}}&\frac{1}{\epsilon_{3}-\epsilon_{2}}&0&\frac{1}{\epsilon_{3}-\epsilon_{4}}&\dots&\frac{1}{\epsilon_{3}-\epsilon_{{M}}}\\ \frac{1}{\epsilon_{4}-\epsilon_{2}}&\frac{1}{\epsilon_{4}-\epsilon_{2}}&\frac{1}{\epsilon_{4}-\epsilon_{3}}&0&\dots&\frac{1}{\epsilon_{4}-\epsilon_{{M}}}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ \frac{1}{\epsilon_{{M}}-\epsilon_{2}}&\frac{1}{\epsilon_{{M}}-\epsilon_{2}}&\frac{1}{\epsilon_{{M}}-\epsilon_{3}}&\frac{1}{\epsilon_{{M}}-\epsilon_{4}}&\dots&0\\ \end{array}\right)$$

$$=$$

$$\underset{\Delta \to 0}{\lim }\det A=\frac{\det A_{\hat{1},\hat{2}}^{\hat{1},\hat{2}}}{{\Delta }^2}.$$

$$\det A=\sum _{j>1}\frac{\det A_{\hat{1},\hat{j}}^{\hat{1},\hat{j}}}{{\left({ϵ}_1-{ϵ}_j\right)}^2},$$

$$K_{\pi }=\left(\prod _{k=1}^{N-r}2g{F}_x({ϵ}_{{\bar{\pi }}_k})F_y({ϵ}_{{\bar{\pi }}_k})\right)\sum _{\sigma \in {\Sigma }_{\bar{\pi }}}\prod _{k=1}^{(N-r)/2}\frac{1}{({ϵ}_{i_k}-{ϵ}_{j_k}{)}^2}$$

Full text

The eigenstates, as constructed here, do not require any usable Bethe ansatz and therefore: no proper pseudo-vacuum, Bethe roots, or generalised spin raising (Gaudin) operators have to be defined. Indeed, the eigenstates are generically built only through the conserved charges which define the model of interest and the specification of the set of eigenvalues defining the particular eigenstate. Since these eigenvalues are, in general, solutions to a simple set of quadratic equations, the proposed approach is simpler to implement than any Bethe ansatz and, moreover, it remains completely identical independently of the symmetries of the model. Indeed, the construction removes any distinction between XYZ models and XXZ/XXX models and, generically, that between models with or without U(1) so that any difficulties associated with the use of a Bethe ansatz in any of these cases are avoided.

1 Introduction

Quantum integrable models of the Richardson-Gaudin (RG) [1, 2, 3, 4] class have already found applications in a remarkably vast array of domains: from superconductivity in metallic grains [5, 6, 7], to nuclear pairing [8, 9, 10] through applications in the behaviour of spin qubits coupled to a spin bath [11, 12, 13, 14, 15, 16]. Their value to problems of physical interest cannot be denied and, as such, the capacity to compute their physical properties as efficiently as possible remains an important question.

In recent years, this class of integrable models has been considerably expanded by relaxing the usual assumption of anti-symmetric couplings. Either through non-skew-symmetric r-matrices [17, 18, 19, 20, 21, 22, 23, 24, 25] or, as in Gaudin’s original work [1, 2] by directly building the commuting conserved charges [26, 27]. It was demonstrated that a wide class of such models remain integrable in the presence of an arbitrary external magnetic field, not only in the XXX and XXZ models, but also for certain fully anisotropic XYZ RG models built out of spins-1/2. While similar constructions are possible for higher spins $s>1/2$ (see [20] for example), they must contain additional local terms which takes them outside of the class of RG models studied here which is assumed to only contain $S^{x}_{i}S^{x}_{j}$, $S^{y}_{i}S^{y}_{j}$ and $S^{z}_{i}S^{z}_{j}$ coupling terms between spins $i$ and $j$.

It was recently demonstrated that, for these spin-1/2 systems in an external field, the conserved charges systematically obey quadratic relations so that the eigenvalues characterising their eigenstates can always be found as the solutions to a small set of quadratic equations [28]. Expectation values of local spin operators, in any of these eigenstates, can then also be computed easily by making use of the Hellmann-Feynman theorem giving access to some of the physics of these systems [29]. All of these results can be obtained without any knowledge of whether, or how, one can construct the eigenstates corresponding to these eigenvalues.

When the system of interest retains a U(1) symmetry (associated to its rotational invariance in a given plane we will call the XY plane), they allow for the use of the Algebraic Bethe Ansatz (ABA) to construct their eigenstates [30, 31, 32, 33]. In doing so, one defines the eigenvalues in terms of a set of Bethe roots $\{\lambda_{1}\dots\lambda_{M}\}$. These Bethe roots, when “on-shell”, i.e. solutions of a system of $M<N$ coupled non-linear algebraic equations dubbed the Bethe equations, provide the various eigenvalues of the conserved charges associated with the corresponding eigenstate. This set of Bethe roots then also provide, through the repeated action $M$ times of a generalised creation operator $B(u)$ a representation of the eigenstate itself as $\prod_{i=1}^{M}B(\lambda_{i})\left|\Omega\right>$, where the specific state $\left|\Omega\right>$ can be called the pseudo-vacuum.

Certain models which break $U(1)$ symmetry either through an external magnetic field [24, 26, 27], or through other integrable “boundaries” [23, 34] which can lead to additional terms in the conserved charges, have also be shown to be solvable in a remarkably similar way despite lacking a “proper” pseudo-vacuum. A construction similar to the ABA has indeed been shown to lead to eigenstates of the form $\prod_{i=1}^{N}B(\lambda_{i})\left|\Omega\right>$ (now containing a number of Bethe roots equal to $N$, the system size) acting on the same reference-state as the one used for a normal ABA implementation. This Modified Algebraic Bethe Ansatz (MABA) [34, 35, 36, 37, 38] therefore also provides a compact representation of the eigenstates of these models. However, for XYZ models, which break $U(1)$ symmetry by being fully anisotropic in their couplings, only the elliptic Gaudin models has, to the best of our knowledge, an explicit construction of its eigenstates and it requires a much more complex ABA [39]. The XYZ RG models in an external magnetic field which are discussed in this work do not naturally give access to an ABA or MABA approach for finding their eigenstates. Nonetheless, explicit expressions for the eigenstates will be built, only by exploiting the quadratic relations between the spin-1/2 conserved charges and, therefore, without having to define any Bethe ansatz procedure. Since the approach completely bypasses any reference to Bethe roots, to the need to specify a pseudo-vacuum state or to the necessity to define proper $B(u)$ operators, it provides what can be dubbed a “Bethe-ansatz-free” solution of the generic (XYZ) spin-$\frac{1}{2}$ Richardson-Gaudin integrable models in an arbitrary magnetic field.

The resulting construction gives explicit expression for the projectors on any eigenstate in terms of only the conserved charges of the model and the set of eigenvalues defining the eigenstate of interest. Since finding these eigenvalues has been reduced to finding solutions of a system of quadratic equations, the resulting formal definition of eigenstates could become of great use numerically, since it only requires the easy to compute (compared to Bethe roots) eigenvalues.

The next section will first describe the models this work is concerned with as well as the relevant quadratic relations between conserved charges and eigenvalues. In the following one, the main result of this work is given, namely an explicit expression for the eigenstate projectors as operator-determinant, after which the validity of the construction is proven and a simple resolution of the identity operator is proposed. Finally, we present an explicit expansion of the determinant which highlights some of its underlying mathematical structure while concluding remarks are given in the very last section.

2 The models

The quantum integrable models of interest in this work are all characterised by the existence of $N$ operators $\hat{R}_{i}$ commuting with each other. Here $N$ is the number of spins (or pseudo-spins) in the system and each of them will be considered as a spin-1/2 representation of an SU(2) algebra, i.e. each of them (labelled by an index $i=1,2\dots N$) are defined by the usual $2\times 2$ Pauli matrices $\sigma^{x}_{i},\sigma^{y}_{i},\sigma^{z}_{i}$.

We impose that the conserved charges have only diagonal-couplings ($x$ direction coupled with $x$ direction, etc) and allow for an external magnetic field, i.e. the conserved charges are chosen to have the form:

[TABLE]

One can then insure integrability by simply making sure that the parameters are chosen in such a way that $\left[\hat{R}_{i},\hat{R}_{j}\right]=0\ \ \forall\ \ i,j$. This imposes restrictions on the coupling constants $X_{ij},Y_{ij},Z_{ij}$ and the field components $B^{x}_{i},B^{y}_{i},B^{z}_{i}$. The solutions to these constraints have been completely parametrised in [29] where it was demonstrated that the most general set of parameters compatible with the requirements of integrability can be written as:

[TABLE]

with the functions $F$ given by:

[TABLE]

for arbitrary $g,\gamma,\lambda,\alpha_{x},\alpha_{y},\beta_{x},\beta_{y}$ and $\epsilon_{i}$. Here, the choice of $B^{z}_{i}=1$ remains generic in that conserved charges can always be rescaled by an arbitrary constant. Since $\frac{\hat{R}_{i}}{B^{z}_{i}}$ still all commute together, any model with $N$ arbitrary non-zero values $B^{z}_{i}$ can always be rescaled to an equivalent $B^{z}_{i}=1$ model. The only case which is therefore excluded is that of models which are integrable exclusively at strictly zero field, namely the elliptic Gaudin models. Indeed, although zero-field limits of the models considered here can be properly defined through an appropriate rescaling of the parameters [29], in the fully anisotropic XYZ case such a limit does not correspond to the usual elliptic Gaudin models. The elliptic models being antisymmetric in their coupling constants while the zero-field limit of the XYZ models of interest here can only lead to non-skew-symmetric XYZ models.

For all of the models parametrised by (2), it was shown [28] that integrability is sufficient to establish that the resulting $N$ conserved charges obey, at the operator level, the following quadratic relations:

[TABLE]

where $\Gamma_{ij}=-\frac{2Z_{ij}Y_{ij}}{X_{ji}}$ and $K_{i}=\sum_{\alpha}B_{\alpha}^{2}+\sum_{j\neq i}\left(X_{ij}^{2}+Y_{ij}^{2}+Z_{ij}^{2}\right)$. Since the general parametrisation, once rescaled to $B^{z}_{i}=1$, systematically leads to $X_{ij}=-Y_{ji}$ one can then always use these adequately rescaled operators and therefore set:

[TABLE]

Since they all commute with one another, the conserved charges share a common eigenbasis in which they are all diagonal and therefore the same quadratic relations are also obeyed by the $N$ eigenvalues $(r_{1}\dots r_{N})$ which specify any given eigenstate. The $n^{\mathrm{th}}$ eigenstate: $\left|\psi_{n}\right>$ such that $\hat{R}_{i}\left|\psi_{n}\right>=r_{i}^{n}\left|\psi_{n}\right>\ \forall\ i=1\dots N$ is then, in principle, fully characterised by $(r_{1}^{n}\dots r_{N}^{n})$, i.e. the $n^{\mathrm{th}}$ solution to the system of quadratic equations:

[TABLE]

In the next section, we demonstrate that the projector $P_{n}\equiv\left|\psi_{n}\right>\left<\psi_{n}\right|$, on any given eigenstate $\left|\psi_{n}\right>$ can be defined in a compact way, exclusively in terms of the conserved charges $\hat{R}_{i}$ and the corresponding set of eigenvalues $r^{n}_{i}$ of that particular eigenstate. The resulting representation of the eigenstates is then obtainable exclusively and directly from the solutions of the quadratic Bethe equations for the eigenvalues (6) whose simplicity provides an important numerical advantage in stability and computational speed [40, 41, 42, 43, 44] over the (M)ABA approaches based on Bethe roots. Even in the XXZ cases (including those without U(1)-symmetry), for which eigenstates could be written in terms of Bethe roots, the proposed representation could lead to major simplification in numerically approaching the physics of these systems.

3 Eigenstate projectors

A direct consequence of the quadratic relations (4) between operators is that an arbitrary polynomial (or even a formal power series) in these conserved charges: $\hat{Q}(\hat{R}_{1}\dots\hat{R}_{N})\equiv\displaystyle\sum_{n_{1},n_{2}\dots n_{N}=0}^{\infty}C_{n_{1}\dots n_{N}}\ \hat{R}_{1}^{n_{1}}\hat{R}_{2}^{n_{2}}\dots\hat{R}_{N}^{n_{N}}$, can always be reduced to:

[TABLE]

a polynomial which is at most linear in each of the conserved charges. As a simple example for $N=2$, one has $\hat{R}_{1}^{3}\hat{R}_{2}=\hat{R}_{1}(\Gamma_{12}\hat{R}_{2}+K_{1})\hat{R}_{2}=\Gamma_{12}\hat{R}_{1}\hat{R}^{2}_{2}+K_{1}\hat{R}_{1}\hat{R}_{2}=\Gamma_{12}\hat{R}_{1}(\Gamma_{21}\hat{R}_{2}+K_{2})+K_{1}\hat{R}_{1}\hat{R}_{2}=(\Gamma_{12}\Gamma_{21}+K_{1})\hat{R}_{1}\hat{R}_{2}+K_{2}\Gamma_{12}\hat{R}_{1}$. Through these successive reduction of quadratic terms to linear ones, the same procedure could be carried out for arbitrary orders and for arbitrary number of conserved charges $N$, systematically leading to polynomials at most linear in each $\hat{R}_{i}$. Such a generic linear polynomial is then, ultimately, defined by a set of $2^{N}$ coefficients $B_{n_{1}\dots n_{N}}$, a number which corresponds to the dimension of the full Hilbert space.

Among this large class of operators parametrised by $2^{N}$ coefficients, we define a subclass, compactly parametrised by only $N$ real variables $(r_{1}\dots r_{N})$ as the following determinant of an “operator-valued” $N$ by $N$ matrix:

[TABLE]

Keeping in mind that every conserved charges commute with one another, the determinant representation of $\hat{P}$, despite having operator entries, has no ambiguity in the order of the products involved. We define as well the following scalar-valued polynomial:

[TABLE]

In the next section, we will prove our main result which is stated here: the projector on the normalised eigenstate $\left|\psi_{n}\right>$ such that $\hat{R}_{i}\left|\psi_{n}\right>=r_{i}^{n}\left|\psi_{n}\right>$, is proportional to the operator $\hat{P}$ evaluated at $(r^{n}_{1}\dots r^{n}_{N})$:

[TABLE]

and can be properly normalised by the scalar determinant defined in (13) so that:

[TABLE]

Consequently, an (unnormalised) representation of the eigenstate itself can be obtained from the action of this projector on an arbitrary “vacuum” state $\left|\Omega\right>$ as:

[TABLE]

with the only requirement that $\left|\Omega\right>$ has a non-zero overlap with the target eigenstate $\left|\psi_{n}\right>$. This is not particularly limiting since, in a generic model, eigenstates will have a non-zero overlap with any canonical spin state $\prod_{i=1}^{r}S^{+}_{i_{r}}\left|\downarrow\downarrow\dots\downarrow\right>$ which could then all serve the role of vacuum. The exception is for models which exclusively contains a z-component of the magnetic field. In these cases, the XXZ (and XXX) models retain the U(1) symmetry and therefore each eigenstate have a fixed z-axis total magnetisation. This would require $\left|\Omega\right>$ to have a component in the appropriate magnetisation sector. Consequently, any $\left|\Omega\right>$ which is spread over every sector would provide a proper “vacuum” for every eigenstate. For XYZ models in a z-only magnetic field, the eigenstates actually have fixed magnetisation parity, meaning they either belong to the sub-Hilbert-space which contains only odd OR even numbers of up-pointing spins. Consequently any $\left|\Omega\right>$ which contains at least one odd and one even magnetisation term (such as $\left|\downarrow\downarrow\dots\downarrow\right>+\left|\uparrow\downarrow\dots\downarrow\right>$ for example) would allow one to build any eigenvector.

While the projectors themselves (21) are indeed the ones built using the normalised eigenstate, the resulting vector (22) would, however, have a normalisation which depends on the particular choice of $\left|\Omega\right>$.

4 Proof of the projector’s representation

Considering that the family of operators $\hat{P}(r_{1}\dots r_{N})$ as defined by eq. (20) are built exclusively out of the conserved charges and scalar parameters, $\hat{P}(r_{1}\dots r_{N})$ is always a diagonal operator in the orthonormalised eigenbasis built from the $\left|\psi_{n}\right>$ eigenstates. Each one of the matrix elements in the determinant (20) is therefore also a diagonal matrix in this basis. The orthogonality of the basis states consequently leads to the following explicit expansion of the determinant in that basis:

[TABLE]

where

[TABLE]

In order to prove that $P(r_{1}\dots r_{N})$ becomes proportional to the desired projector when evaluated at $(r_{1}^{n}\dots r_{N}^{n})$, it is sufficient to prove that, at this point, every coefficient $C_{m\neq n}(r^{n}_{1}\dots r^{n}_{N})$ cancels out.

Considering that, being eigenvalues $(r_{1}^{n}\dots r_{N}^{n})$ and $(r_{1}^{m}\dots r_{N}^{m})$ are “on-shell”, i.e. they form a solution to the quadratic Bethe equations (6), one has:

[TABLE]

it then becomes easy to see that, when the matrix $J$ given in eq. (30) is evaluated at $(r_{1}\dots r_{N})=(r_{1}^{n}\dots r_{N}^{n})$, its columns are not independent since the following linear combination of its columns:

[TABLE]

Consequently, the determinant, and therefore $C_{m}(r^{n}_{1}\dots r^{n}_{N})=0$ for every $n\neq m$. The only exception is found when $n=m$ since every coefficient in the linear combination is then explicitly $(r_{j}^{n}-r_{j}^{n})=0$. In this case, the non-zero value of $C_{n}(r^{n}_{1}\dots r^{n}_{N})$ is obviously given by the determinant polynomial $N(r_{1}^{n}\dots r_{N}^{n})$ defined in (13). This therefore proves the proposition: the operator $\hat{P}(r_{1}\dots r_{N})$ introduced in eq. (10), when evaluated at a point $(r_{1}\dots r_{N})$ which corresponds to a solution of the quadratic Bethe equations (6), become the projector on the normalised eigenstate corresponding to these eigenvalues:

[TABLE]

This provides an operator-based concept of “on-shell” and “off-shell”, where the generic $N$-parameter operator $\hat{P}(r_{1}\dots r_{N})$, becomes the projector on an eigenstate of the integrable model whenever its $N$ free parameters are on-shell, i.e. chosen to be a solution of (quadratic) Bethe equations. This is in complete analogy with the usual ABA for which a generic Bethe state $\prod_{i=1}^{M}B(\lambda_{i})\left|\Omega\right>$ becomes an eigenstate whenever the Bethe roots $\lambda_{i}$ are on-shell, i.e.: solution to the Bethe equations. However, the construction proposed in this work avoids the need to define the relation between Bethe roots and eigenvalues of the conserved charges, the need to define the right $B(u)$ operators and the necessity to find a proper pseudo-vacuum $\left|\Omega\right>$ on which to create the excitations. In this sense, it gets completely rid of the possible difficulties of applying a Bethe ansatz and provides, nonetheless, the exact solution to the eigenproblem of spin-1/2 Richardson-Gaudin integrable systems even in the cases where Bethe roots-based approaches typically cannot be used. The proposed solution remains valid in every case without U(1) symmetry (apart from the zero-field elliptic model which does not belong to the class of models for which the proposed quadratic equations are respected). Among these, the XXZ (and XXX) systems in an arbitrary external field could also be approached using the (M)ABA, but the XYZ case in an arbitrary field did not, to the best of our knowledge, have an explicit solution for its eigenstates yet.

4.1 Off-shell resolution of the identity and product of operators

On top of providing, when on-shell, a representation for the eigenstate projectors, the generic (off-shell) family of operators $\hat{P}(r_{1}\dots r_{N})$ also allows one to build simple resolutions of the identity. Indeed, as it was defined in (10), $\hat{P}(r_{1}\dots r_{N})$ is not only a polynomial in the conserved charges $\hat{R}_{i}$, but also a polynomial, of maximal order 1, in each of the scalar variables $r_{i}$. Moreover, the (operator) coefficient in front of the $\left(\displaystyle\prod_{i=1}^{N}r_{i}\right)$ term is given by $\mathbb{1}$ since this term, as seen by expanding the determinant, comes exclusively from the product of the $N$ diagonal elements. Consequently, integrating $\Psi(r_{1}\dots r_{N})\hat{P}(r_{1}\dots r_{N})$ on an N-dimensional hypercube centered around the origin will provide, for any function $\Psi(r_{1}\dots r_{N})$ antisymmetric in every variable, a resolution of the identity. A simple example would be to choose $\Psi(r_{1}\dots r_{N})=\prod_{i=1}^{N}r_{i}$ and integrate, on the real axis, on the symmetric intervals $r_{i}\in\left[-a,a\right]$. This leads to

[TABLE]

since every term in the expansion of $\hat{P}(r_{1}\dots r_{N})$, except for the leading term $\left(\prod_{i}r_{i}\right)$, is missing at least one of the $r_{i}$ variables. Apart from the leading term, every term in $\Psi(r_{1}\dots r_{N})\hat{P}(r_{1}\dots r_{N})$ is odd in the missing variables and their integral over the symmetric intervals will therefore cancel out.

Another important quantity which is simple to evaluate is the product of an “on-shell” and “off-shell” operator, analogous to Slavnov’s determinants [31, 32, 49, 50]. In fact, considering that we have an explicit decomposition of the “off-shell” operator $\hat{P}(r_{1}\dots r_{N})$ given in (23) we directly know that its product with an “on-shell” projector is given by the corresponding coefficient (30), namely:

[TABLE]

This result is highly reminiscent of the eigenvalue-based determinants known for XXX [27, 51] and XXZ [42, 50] models, being again defined simply as determinants in which eigenvalues $r^{n}_{i}$ and their “off-shell” equivalent $r_{i}$ only appear as their sum $(r_{i}+r^{n}_{i})$ and do so exclusively in the diagonal elements of an $N$ by $N$ matrix.

5 Explicit expansion of $\hat{P}(r_{1}\dots r_{N})$

In this last section, we study the explicit expansion of the family (10) of operators $\hat{P}(r_{1}\dots r_{n})$ in a way which provides further insight into their underlying structure. Considering the general expression of a determinant as a sum over permutations:

[TABLE]

the first term present is given by the “correctly” ordered permutation $\sigma_{1}=1,\sigma_{2}=2\dots\sigma_{N}=N$ and leads, in the determinant (10), to the term $\prod_{i=1}^{N}(r_{i}+R_{i})$. Permutations which exchange only two of these indices $(\sigma_{k^{\prime}}=k,\sigma_{k}=k^{\prime})$ would then lead to the replacement of two diagonal elements by the corresponding off-diagonal elements, while changing the sign of the permutation. These then give a series of terms $\displaystyle-\sum_{k\neq k^{\prime}}^{N}\left[\Gamma_{kk^{\prime}}\Gamma_{k^{\prime}k}\prod_{i\neq\{k,k^{\prime}\}}^{N}(r_{i}+R_{i})\right]$, etc. Since the determinant is a sum over products of $N$ matrix elements $\hat{J}_{ij}$ where the $i$ indices need to cover every line and the $j$ indices have to cover every column, one can classify the resulting terms according to the number of diagonal matrix elements (when $\sigma_{k}=k$) which are used to build each individual term. One can then write the determinant as a sum over the subsets $\pi$ of every possible cardinality $r=0,1\dots N$, that one can build out of $\{1,2,3\dots N\}$. A given set $\pi=\{\pi_{1},\pi_{2}\dots\pi_{r}\}$ defines the diagonal matrix elements $\hat{J}_{\pi{{}_{k}}\pi{{}_{k}}}$ present and the complement $\bar{\pi}=\{\bar{\pi}_{1}\dots\bar{\pi}_{N-r}\}$ (such that $\pi\cup\bar{\pi}=\{1,2,3\dots N\}$) defines the indices which will be used for off-diagonal elements. Defining $\Pi$ as the power set of $\{1,2,3\dots N\}$ (containing every one of its subsets of every cardinality), one can write:

[TABLE]

with $K_{\pi}$ given by the $M\times M$ determinant (with $M\equiv N-r$):

[TABLE]

built from the $M$ indices in the complement $\{\bar{\pi}_{1}\dots\bar{\pi}_{M}\}$. Remarkably, the structure of $\Gamma_{i,j}$ imposed by integrability is such that, for $M$ odd, the resulting coefficient is strictly zero. Indeed, using the general parametrisation (2) and the fact that $\Gamma_{ij}=2Z(\epsilon_{i},\epsilon_{j})$, one can write the determinant as:

[TABLE]

showing it is proportional to that of a square real skew-symmetric matrix $A$ of dimension $N-r$. From Jacobi’s theorem, we then know that the resulting determinant is zero for any odd $N-r$ dimension. For even dimensions, “Cayley’s theorem for pfaffians” states that the determinant of a skew-symmetric matrix is given by the square of its pfaffian. However, in the case at hand, the integrability of the models, and the constraints it imposes on the matrix elements $\Gamma_{ij}$, the non-zero determinants found for even $N-r$ can be further simplified into a single sum over partitions into pairs. Redefining, for compactness, and without loss of generality, the indices $\bar{\pi}_{i}\to i$, we are interested in the determinant of the matrix $A_{i,j}=\frac{1-\delta_{ij}}{\epsilon_{i}-\epsilon_{j}}$.

By the very definition of a determinant, one has to encounter twice the index 1 in each term of its expansion, one matrix element having to come from the first line and one from the first column. Looking at a given specific pairing, say $(\epsilon_{1},\epsilon_{2})$, the expansion of the determinant can therefore contain: a) a term quadratic in $\frac{1}{(\epsilon_{1}-\epsilon_{2})}$ coming from $A_{12}A_{21}$, b) terms linear in $\frac{1}{(\epsilon_{1}-\epsilon_{2})}$ which comes from either $A_{12}$ or $A_{21}$ while the other instance of index 1 comes from a different line or column than $2$ and finally c) terms which do not contain $\frac{1}{(\epsilon_{1}-\epsilon_{2})}$, i.e. terms which do not contain $A_{12}$ or $A_{21}$.

Such a classification of terms allows one to study the analytic structure of the determinant by looking at $\epsilon_{1}=\epsilon_{2}+\Delta$ in the limit $\Delta\to 0$ in order to prove that the determinant has no simple pole at $\epsilon_{1}\to\epsilon_{2}$ and therefore that terms linear in $\frac{1}{(\epsilon_{1}-\epsilon_{2})}$ all sum up to zero.

Expanding the determinant in its first two rows, it can be written as:

[TABLE]

where $A^{\hat{1},\hat{2}}_{\hat{j},\hat{k}}$ is the matrix $A$ from which line 1 and 2 as well as columns $j$ and $k$ have been removed while $\epsilon^{jk}$ is the sign of the corresponding partition of the indices.

First, one easily sees that the two by two determinants \mathrm{det}\left(\begin{array}[]{cc}\frac{1}{\epsilon_{2}-\epsilon_{j}}&\frac{1}{\epsilon_{2}-\epsilon_{k}}\\ \frac{1}{\epsilon_{2}-\epsilon_{j}}&\frac{1}{\epsilon_{2}-\epsilon_{k}}\end{array}\right) are all zero. Moreover, since a line permutation and the resulting minus sign makes \mathrm{det}\left(\begin{array}[]{cc}0&\frac{1}{\epsilon_{2}-\epsilon_{j}}\\ \frac{1}{\Delta}&\frac{1}{\epsilon_{2}-\epsilon_{j}}\end{array}\right)=\mathrm{det}\left(\begin{array}[]{cc}-\frac{1}{\Delta}&\frac{1}{\epsilon_{2}-\epsilon_{j}}\\ 0&\frac{1}{\epsilon_{2}-\epsilon_{j}}\end{array}\right), the change in the sign of the partition: $\epsilon^{1k}=-\epsilon^{2k}$, makes both single sums over $j$ also trivially cancel term by term. All in all, one simply has:

[TABLE]

Consequently, the determinant has no simple pole at $\epsilon_{1}\to\epsilon_{2}$; only a double pole $\frac{1}{(\epsilon_{1}-\epsilon_{2})^{2}}$. Since the exact same reasoning can be carried out for any $\frac{1}{\epsilon_{1}-\epsilon_{j}}$ ($j=2,3,4\dots N$), the resulting determinant can only depend on $\epsilon_{1}$ through double poles in all of its possible pairings. It can therefore be written as:

[TABLE]

since $\epsilon_{1}$ has to appear only twice in each term of the determinant (matrix elements coming from line and column 1), and the possible terms $\frac{1}{\left(\epsilon_{1}-\epsilon_{j}\right)^{2}}$ exhaust those two appearances of $\epsilon_{1}$.

This leads to a recursive construction of the determinant, since $A^{\hat{1},\hat{j}}_{\hat{1},\hat{j}}$ has the exact same structure as the original matrix $M$ albeit with a dimension lowered by two and the absence of the indices $1$ and $j$. The exact same reasoning could then be applied to the possible pairings of the next unpaired index and the procedure can be kept going until $(N-r)/2$ pairs of indices have been made. The resulting determinant is then given by the sum (without sign alternances) over all the possible pairings. Since $(-1)^{N-r}=1$, one therefore has:

[TABLE]

with the sum carried over $\Sigma_{\bar{\pi}}$, i.e.: the set of possible partitions in pairs $\sigma=\{\{i_{1},j_{1}\},\{i_{2},j_{2}\}\dots\{i_{(N-r)/2},j_{(N-r)/2}\}\}$ of the indices contained in $\bar{\pi}$.

The generic determinant operator $\hat{P}(r_{1}\dots r_{N})$, which become a projector when on-shell, has therefore been shown to have an explicit expansion in terms of only odd (for $N$ odd) or even (for $N$ even) products of $(r_{i}+\hat{R}_{i})$. Moreover, the coefficient in front of each of these products is given by the sum, over possible pairings of the complementary indices, of the products of factors $\frac{1}{(\epsilon_{i}-\epsilon_{j})^{2}}$, so that:

[TABLE]

where, for $N$ odd (even), $\tilde{\Pi}$ is the set of all subsets of $\{1,2\dots N\}$ having exclusively odd (even) cardinality, while $K_{\pi}$ is given by eq. (77). This result is to be contrasted with the square of the Pfaffian one would get from the determinant of a generic skew-symmetric matrix which would then be expanded a a double sum over these partitions into pairs.

As a simple explicit example for $N=4$ spins, the expansion of $\hat{P}(r_{1},r_{2},r_{3},r_{4})$ only contains $r=4,r=2$ and $r=0$ terms and is explicitly given by:

[TABLE]

While these expanded representations are certainly less compact than the determinant representation (10) they reveal explicitly how the integrability of the system, and the resulting relations between $\Gamma$s come into play to in defining some of the fundamental underlying structure of the (eigen-)states constructed in this work.

6 Conclusions

In this article, we have showed how one can systematically build the eigenstates of spin-1/2 Richardson-Gaudin models in a magnetic field by only specifying the set of conserved charges and the eigenvalues associated with the given eigenstate. To do so, we defined a class of operators, parametrised by $N$ variables $(r_{1}\dots r_{N})$, which, when evaluated “on-shell”, i.e. at $(r_{1}\dots r_{N})$ solution to quadratic Bethe equations, explicitly become the projector on the corresponding normalised eigenstate.

This result completes the eigenvalue-based approach to these models by providing generic expressions for the eigenstates without having to define a Bethe Ansatz. One can hope that this construction could be used to provide expressions, written directly in terms of eigenvalues, for form factors of local operators which form the basic building blocks required to access the dynamical properties of these systems.

The authors would like to thank Pieter W. Claeys for all of his valuable comments and suggestions on this work.