Entanglement and matrix elements of observables in interacting integrable systems

Much work has been done in the past decade to understand the far-from-equilibrium dynamics and description after equilibration of isolated nonintegrable (generic) and integrable quantum many-body systems ETH_review ; eisert_friesdorf_review_15 ; polkovnikov_sengupta_review_11 . Despite tremendous progress in recent years essler_fagotti_2016 ; calabrese_cardy_2016 ; cazalilla_chung_2016 ; bernard_doyon_2016 ; caux_2016 ; GGE_review ; ilievski_medenjak_2016 ; langen_gasenzer_2016 ; vasseur_moore_2016 ; deluca_mussardo_2016 , the microscopics of interacting integrable systems are those that remain less understood. On the one hand, interactions are present in those systems as in nonintegrable ones (making their study challenging), and on the other hand, they exhibit extensive numbers of local conserved quantities as noninteracting systems do. The presence of such quantities precludes thermalization in integrable (noninteracting and interacting) quantum systems Kinoshita2006 ; GGE_Rigol_2007 ; wouters_denardis_14 ; pozsgay_mestyan14 ; ilievski_denardis_15 .

Thermalization does occur in generic isolated quantum systems, and is understood to be a consequence of quantum chaos and eigenstate thermalization Deutsch1991 ; Srednicki1994 ; ETHRigol ; ETH_review . Essentially, the matrix elements of observables (few-body operators) $\hat{O}$ in eigenstates of generic quantum Hamiltonians are described by the eigenstate thermalization hypothesis (ETH) ansatz ETH_review ; Srednicki_1999 :

where $\bar{E}\equiv(E_{\alpha}+E_{\beta})/2$, $\omega=E_{\alpha}-E_{\beta}$, and $S(\bar{E})$ is the thermodynamic entropy at energy $\bar{E}$. The functions $\mathcal{O}(\bar{E})$ and $f_{O}(\bar{E},\omega)$ are smooth, and $R_{\alpha\beta}$ is a normally distributed variable with zero mean and unit variance (variance 2) for $\alpha\neq\beta$ ($\alpha=\beta$) in Hamiltonians that exhibit time-reversal symmetry. The smoothness of the diagonal matrix elements allows observables described by Eq. (1) to thermalize (i.e., to be described by traditional ensembles of statistical mechanics) for experimentally relevant initial conditions. The off-diagonal matrix elements control the approach to and fluctuations about equilibrium ETH_review .

The ETH ansatz (1) has been extensively tested in exact diagonalization studies of nonintegrable Hamiltonians ETH_review ; ETHRigol ; Breakdown ; rigol_09b ; Santos2010Localization ; Khatami ; Steinigeweg2013Eigenstate ; Beug1 ; steinigeweg_khodja_14 ; sorg_vidmar_14 ; Kim2014Testing ; Beug2 ; 2DTFIM ; 2DTFIM-2 ; Sagawa2018 ; khaymovich_haque_19 ; Jansen2019Eigenstate ; mierzejewski_vidmar_19 . While the diagonal matrix elements display a shrinking support and an exponentially decaying variance with increasing system size in nonintegrable systems, confirming that the diagonal ETH is valid for them, the same is not true in integrable systems Breakdown ; rigol_09b ; Santos2010Localization ; Beug1 ; GGE_Rigol_2011 ; GGE_review ; mierzejewski_vidmar_19 . In integrable systems, the support of the diagonal matrix elements need not shrink, and their variance is expected to decay as a power law in system size PhysRevLett.105.250401 ; Ikeda ; Alba . The latter is consistent with the fact that the variance of the diagonal matrix elements of intensive translationally invariant observables must vanish at least as a power law in system size because it is bounded from above by the Hilbert-Schmidt norm (which vanishes as a power law in system size) mierzejewski_vidmar_19 .

The off-diagonal matrix elements of observables in the eigenstates of interacting integrable quantum many-body systems have received little attention. While such results have been reported for specific models and system sizes alongside those of quantum chaotic systems rigol_09b ; Santos2010Localization ; Beug2 , there has been no systematic study of their properties. For noninteracting models (or models mappable to them), the existence of an increasingly large (with increasing system size) fraction of vanishing off-diagonal matrix elements precludes the definition of a meaningful function $f_{O}(\bar{E},\omega)$, in contrast to quantum chaotic models Khatami . On the other hand, recent results by two of us (K.M. and M.R.) in the context of periodically driven systems provided strong evidence that one can define (and experimentally measure) a function $|f_{O}(\bar{E},\omega)|^{2}=e^{S(\bar{E})}\overline{|O_{\alpha\beta}|^{2}}$ for interacting integrable models 1907.04261 . Exploring this, along with other properties of the off-diagonal (and diagonal) matrix elements in the spin-1/2 XXZ chain, is one of the two central goals of this work.

The other central goal of this work is to study the structure of highly excited energy eigenstates by means of their bipartite entanglement. Recently, much work has been devoted to understanding entanglement properties of highly-excited eigenstates of many-body Hamiltonians alba09 ; deutsch_10 ; santos_polkovnikov_12 ; hamma_santra_12 ; deutsch13 ; moelter_barthel_14 ; storms14 ; beugeling15 ; yang_chamon_15 ; lai15 ; nandy_sen_16 ; ee_fermion ; ee_chaotic ; dymarsky_lashkari_18 ; riddell_muller_18 ; zhang_vidmar_18 ; liu_chen_18 ; garrison_grover_18 ; nakagawa_watanabe_18 ; vidmar_hackl_18 ; huang_19 ; hackl_vidmar_19 ; murciano_ruggiero_19 ; lu_grover_19 ; sun_nie_19 ; giovenale_pont_19 ; bertini_kos_19 ; huang_gu_19 ; murthy_srednicki_19 ; morampudi_chandran_18 ; modak_nag_19 ; miao_barthel_19 ; bianchi_dona_19 ; faiez_sefranek_19 ; faiez_sefranek_19b ; jafarizadeh_rajabpour_19 . Here we study the average entanglement entropy over all eigenstates of the spin-1/2 XXZ Hamiltonian in the zero-magnetization sector. We argue that this average universally reveals the fundamentally different natures of interacting integrable and quantum chaotic models. While in both nonintegrable and interacting integrable systems the leading term of the average entanglement entropy exhibits a volume-law scaling, we show that the corresponding volume-law coefficient is markedly different between the two. In quantum chaotic systems ee_chaotic it matches the prediction by Page page for random pure states in the Hilbert space, while it is smaller for interacting integrable systems. Remarkably, our results for interacting spin-1/2 integrable systems are consistent with this coefficient being very close to, or the same as, the one for translationally invariant free ee_fermion and (more generally) quadratic vidmar_hackl_18 ; hackl_vidmar_19 fermionic Hamiltonians. This suggests that, entanglement-wise, the overwhelming majority of the eigenstates are very similar between interacting spin-1/2 integrable Hamiltonians and noninteracting fermionic Hamiltonians.

The presentation is organized as follows: In Sec. II, we discuss the specific integrable and nonintegrable models and observables considered, as well as details about the numerical calculations carried out. In Sec. III, we compare the average entanglement entropy of eigenstates of the spin-1/2 XXZ chain with that of eigenstates of noninteracting and nonintegrable models. In Sec. IV, we discuss the distributions and scaling properties of the diagonal matrix elements of two local observables at the center of the spectrum. In Sec. V, we discuss the off-diagonal matrix elements of the same observables: their distributions, scaling properties, and functional dependence of $|f_{O}(\bar{E},\omega)|^{2}$ on $\omega$, for $\bar{E}$ at the center of the spectrum. Lastly, in Sec. VI, we summarize our results.

We study the spin-1/2 XXZ chain with the addition of next-nearest-neighbor interactions, with $L$ sites and periodic boundary conditions. The Hamiltonian is

where $\hat{S}_{i}^{\nu}$ are spin-1/2 operators in the $\nu\in\{x,y,z\}$ directions on site $i$, and $\hat{S}_{i}^{\pm}=\hat{S}_{i}^{x}\pm i\hat{S}_{i}^{y}$ are the corresponding ladder operators. When $\lambda=0$, Hamiltonian (2) is integrable and can be solved exactly using the Bethe ansatz RevModPhys.83.1405 . When $\lambda\neq 0$, Hamiltonian (2) is quantum chaotic PhysRevE.81.036206 . We set $\lambda=0$ and $1$ to compare the integrable and nonintegrable regimes, respectively. Unless otherwise specified, we show results for $\Delta=0.55$ and $\Delta=1.1$ to illustrate that they are qualitatively similar in the (nearest-neighbor) easy-plane ($\Delta<1$) and easy-axis ($\Delta>1$) regimes.

We study the matrix elements of two local operators: The nearest-neighbor $z$-interaction

and the next-nearest-neighbor flip-flop operator

To study the entanglement entropy of energy eigenstates, as well as the matrix elements of $\hat{A}$ and $\hat{B}$ in those eigenstates, it is important to resolve all the symmetries of the Hamiltonian ETH_review ; Santos2010Localization . First, we note that Hamiltonian (2) conserves the total magnetization in the $z$-direction, $M^{z}=\sum_{i}\hat{S}_{i}^{z}$. In this work, we focus on the zero magnetization sector in chains with even numbers of lattice sites. The next important symmetry is translation, which allows one to block-diagonalize the Hamiltonian in different total quasimomentum $k$ sectors. All quasimomentum sectors are used in the average entanglement entropy calculations reported in Sec. III. Within the $M^{z}=0$ and $k=0$ sector, there are two additional symmetries, namely, spin inversion ($Z_{2}$) and space reflection ($P$). In our studies of matrix elements, we focus on the even-$Z_{2}$ even-$P$ sector within the $M^{z}=0$ and $k=0$ sector. We use full exact diagonalization of periodic chains with up to $L=26$ sites. The even-$Z_{2}$ even-$P$ sector within the $M^{z}=0$ and $k=0$ sector of the chain with $L=26$, the largest considered, has $101,340$ states.

In this section, we study the entanglement properties of eigenstates $\{|\alpha\rangle\}$ of Hamiltonian (2) in the zero magnetization sector. We consider a bipartition into a subsystem $A$ and its complement $\bar{A}$ that consist of $L_{\rm A}$ and $L-L_{\rm A}$ consecutive lattice sites, respectively. We calculate the bipartite entanglement entropy of an eigenstate $|\alpha\rangle$ as

where $\hat{\rho}_{A}={\rm Tr}_{\bar{A}}\{|\alpha\rangle\langle\alpha|\}$ is the reduced density matrix of the subsystem $A$. We average $S_{\alpha}$ over all Hamiltonian eigenstates in the zero magnetization sector to obtain the average entanglement entropy $\bar{S}={\cal D}^{-1}\sum_{\alpha}S_{\alpha}$, where ${\cal D}=\binom{L}{L/2}$.

The upper bound for the entanglement entropy of pure states in a given magnetization sector (or, equivalently, in a given particle-number sector when mapping spin-1/2 systems onto hard-core bosons or spinless fermions), for a given $L_{\rm A}/L$, depends both on the magnetization and on $L_{\rm A}/L$. The leading term, which scales with the volume, depends on the magnetization ee_chaotic ; garrison_grover_18 . There is also a subleading, $O(1)$, term that depends on $L_{\rm A}/L$ ee_chaotic .

In the zero magnetization sector, for $L_{A}\leq L/2$, the leading and first subleading terms in the average entanglement entropy of random pure states with normally distributed real coefficients are ee_chaotic

On the r.h.s. of Eq. (III), the first two terms are the upper bound for the entanglement entropy of pure states in the $M^{z}=0$ sector ee_chaotic , while the third term is the generalization of the correction derived by Page page for the $M^{z}=0$ sector of a system with conserved $M^{z}$. Motivated by the numerical results in Ref. ee_chaotic , we think of $\bar{S}_{\rm ran}\left(L_{\rm A},\frac{L_{\rm A}}{L}\right)$ as an upper bound for the average entanglement entropy over all eigenstates of any given physical Hamiltonian, with $\bar{S}_{\rm ran}^{\rm max}\equiv\bar{S}_{\rm ran}\left(\frac{L}{2},\frac{1}{2}\right)$ as the maximum. The dashed line in Fig. 1(a) shows $\bar{S}_{\rm ran}\left(L_{\rm A},\frac{L_{\rm A}}{L}\right)$ in the thermodynamic limit.

On the opposite (low-entropy) side of physical Hamiltonians, one has noninteracting (free) fermions. Translationally invariant free fermionic Hamiltonians exhibit the same leading term of the average entanglement entropy as the XX model, Eq. (2), with $\Delta=\lambda=0$ hackl_vidmar_19 . It was proved in Ref. ee_fermion that the leading (volume) term of the average entanglement entropy over all (i.e., including all particle-number sectors) eigenstates in those models is

with $s_{\rm free}\left(0\right)=1$, and $0<s_{\rm free}\left(\frac{L_{\rm A}}{L}\right)<1$ for $L_{\rm A}/L>0$. In Ref. ee_fermion , $s_{\rm free}\left(\frac{L_{\rm A}}{L}\right)$ was computed numerically [dashed-dotted line in Fig. 1(a)] and, for $L_{\rm A}=L/2$, it was found that $s_{\rm free}\left(\frac{1}{2}\right)=0.5378(1)$. Subsequently, it was conjectured that $s_{\rm free}\left(\frac{L_{\rm A}}{L}\right)$ is universal for all translationally invariant quadratic fermionic models vidmar_hackl_18 ; hackl_vidmar_19 . The horizontal lines in Figs. 1(b) and 1(c) show $s_{\rm free}\left(\frac{1}{2}\right)/\bar{S}_{\rm ran}^{\rm max}$.

In Fig. 1(a), we show the average entanglement entropy over all eigenstates within the half-filled sector of noninteracting fermions, as well as within the zero-magnetization sector of integrable ($\Delta=0.55$ and $\lambda=0$) and nonintegrable ($\Delta=0.55$ and $\lambda=1$) points of Eq. (2), for chains with $L=22$. The results at the nonintegrable point are closest to the thermodynamic limit ones for random pure states. Figure 1(b), for $L_{\rm A}/L=1/2$, shows that the small differences seen in Fig. 1(a) appear to vanish in the thermodynamic limit, in agreement with complementary results reported in Ref. ee_chaotic . For free fermions, on the other hand, the results in Fig. 1(a) are closest to the thermodynamic limit ones obtained in Ref. ee_fermion by averaging over all fillings. Figure 1(b), for $L_{\rm A}/L=1/2$, shows that the differences seen in Fig. 1(a) appear to vanish in the thermodynamic limit as expected (the zero magnetization sector is the one dominant in the thermodynamic limit when averaging over all fillings).

The numerical results at the interacting integrable point ($\Delta=0.55$) in Fig. 1(a) are in between the ones for random pure states and the noninteracting ones. However, the finite-size scaling analysis reported in Fig. 1(c) for $L_{\rm A}/L=1/2$ suggests that the leading term of the average entanglement entropy at $\Delta=0.55$ (easy-plane regime), is very close to, or the same as, the one for noninteracting fermions. As a matter of fact, finite-size scaling analyses in Fig. 1(c) for $\Delta=1.1$, 1.45 (easy-axis regime) and $\Delta=1.0$ (Heisenberg point, the most symmetric point in the spin-1/2 XXZ model) suggest that this is true independently of the value of $\Delta$.

The finite-size scaling analyses in Figs. 1(b) and 1(c) suggest that the qualitatively new effect of interactions in integrable systems is subleading, as they change the first subleading term from $O(1)$ in noninteracting models [the leading correction to $\bar{S}/\bar{S}_{\rm ran}^{\rm max}$ in Fig. 1(b) is $\propto 1/L_{\rm A}$] to $O(\sqrt{L_{\rm A}})$ in interacting integrable models [the leading correction to $\bar{S}/\bar{S}_{\rm ran}^{\rm max}$ in Fig. 1(c) is $\propto 1/\sqrt{L_{\rm A}}$].

In this section, we study expectation values of observables in eigenstates of interacting integrable and nonintegrable Hamiltonians, referred to in what follows as the eigenstate expectation values (EEVs) of the observables, in the even-$Z_{2}$ even-$P$ sector within the $M^{z}=0$ and $k=0$ sector (see Sec. II).

In Fig. 2, we show the EEVs of observables $\hat{A}$ and $\hat{B}$ as functions of the eigenenergies per site ($E_{\alpha}/L$) for different chain sizes at integrable [Figs. 2(a) and 2(c)] and nonintegrable [Figs. 2(b) and 2(d)] points of Hamiltonian (2). At the nonintegrable point, for both observables, one can see in the plots that the support (maximum spread) of the EEVs around each $E_{\alpha}/L$ (away from the edges of the spectrum) shrinks upon increasing the chain size $L$. The scaling of a quantifier of the support, see the insets in Figs. 2(b) and 2(d), indicates that the support vanishes exponentially fast with increasing $L$. This suggests that, in the thermodynamic limit, the EEVs are described by the smooth function $O(E)$, which, in turn, is the thermal expectation value of observable $\hat{O}$ at energy $E$ ETH_review . Hence, one expects all EEVs at the nonintegrable point away from the edges of the spectrum to be thermal in the thermodynamic limit. On the other hand, at integrability, Figs. 2(a) and 2(c) show that the support of the EEVs is wide and does not shrink with increasing system size [see the insets in Figs. 2(a) and 2(c)]. The wide nonshrinking support indicates that at the integrable point ETH is not satisfied, as nonthermal states persist in the thermodynamic limit. The next question to address at the integrable point is how those EEVs are distributed.

In Fig. 3, we show the normalized distribution (color coded) of the EEVs for observables $\hat{A}$ and $\hat{B}$ for two different system sizes ($L=22$ and 26) at the integrable point in Fig. 2, along with the microcanonical averages (solid lines) for the respective observables. The microcanonical averages are calculated using the results from $L=26$ in an energy window $\delta E$ that is small enough to yield a smooth curve independent of $\delta E$. By comparing the results for $L=22$ and $L=26$ for each observable, one can see that, despite the large non-vanishing support, the distribution of EEVs (on the $y$-axis) becomes increasingly peaked (with increasing system size) about the microcanonical average (further evidence for this is reported in Fig. 5). Similarly, on the $x$-axis, the distribution becomes increasingly peaked about the center of the spectrum ($E_{\alpha}/L=0$). The latter occurs because of the known Gaussian behavior of the density of states in local Hamiltonians ETH_review ; Hartmann2004 ; Hartmann2005 .

In Fig. 4, we plot the many-body density of states DOS($E_{\alpha}$) as a function of the energy per site $E_{\alpha}/L$ at integrable and nonintegrable points along with Gaussian functions that have the same mean and variance as the data. The Gaussian functions agree well with the numerical results. This makes it apparent that, even for the small chains that one can solve using full exact diagonalization, the density of states is very close to a Gaussian function (away from the edges of the spectrum). This is true regardless of whether the Hamiltonian is integrable or not. The insets show that the variances in our calculations decay as a power law in $L$, with a power that is close to the expected $L^{-1}$ behavior. This shows that, with increasing system size, the overwhelming majority of eigenstates of local Hamiltonians are at the center of the spectrum (with a vanishing energy per site in our case). In what follows, we focus our scaling analyses on that region of the spectrum ($E_{\alpha}/L\simeq 0$).

To quantify the differences seen in Fig. 2 between the EEVs of observables in integrable and nonintegrable systems, we study the average of the absolute value of the eigenstate-to-eigenstate fluctuations Kim2014Testing ; 2DTFIM ; Jansen2019Eigenstate

where the index $\alpha$ labels the eigenenergies $E_{\alpha}$ (sorted in increasing order), and $\overline{|\dots|}$ denotes an average over the central 20% of eigenstates. To carry out an accurate comparison between our results for $\overline{|\delta O_{\alpha\alpha}|}$ and the ETH ansatz for quantum chaotic systems, a modification needs to be introduced to the ansatz in order to tailor it to our observables of interest. This modification is related to the fact that we focus on intensive operators that are defined via extensive sums, as seen in Eqs. (3) and (4), in the presence of translational invariance. The Hilbert-Schmidt norm of those operators scales as $1/\sqrt{L}$ mierzejewski_vidmar_19 , as opposed to the $O(1)$ Hilbert-Schmidt norm one has in mind when writing Eq. (1). As a result, for the diagonal part of our operators $\hat{A}$ and $\hat{B}$, Eq. (1) needs to be rewritten as:

Since we are focusing on the regime $E_{\alpha}/L\simeq 0$, in which $S(E_{\alpha})\simeq\ln(D)$ ($D$ is the dimension of the Hilbert space of the symmetry sector studied), we expect the average eigenstate-to-eigenstate fluctuations of $\hat{O}$ to be $\propto{(LD)}^{-1/2}$. For integrable systems, on the other hand, the average eigenstate-to-eigenstate fluctuations are expected to be proportional to $1/\sqrt{L}$ Alba .

In Fig. 5, we show the finite size scaling of $\overline{|\delta A_{\alpha\alpha}|}$ [Figs. 5(a) and 5(b)] and $\overline{|\delta B_{\alpha\alpha}|}$ [Figs. 5(c) and 5(d)] at two integrable [Figs. 5(a) and 5(c), $\lambda=0$] and two nonintegrable [Figs. 5(b) and 5(d), $\lambda=1$] points. For the nonintegrable points, we observe a near perfect scaling $\propto 1/\sqrt{LD}$ for both observables. This is in agreement with the ETH ansatz and indicates that the scaling is robust against the parameters of the model and the choice of observables. On the other hand, at integrability, $\overline{|\delta A_{\alpha\alpha}|}$ and $\overline{|\delta B_{\alpha\alpha}|}$ exhibit a much slower decay with increasing system size, and also exhibit very strong finite-size effects. While we expect the decay to be $\propto 1/\sqrt{L}$ Alba , this is not the exponent of the power law we find if we fit the data to $c_{1}L^{-c_{2}}$. Instead, we have fitted the data to the function $\text{fit}(L)=c_{1}/\sqrt{L}+c_{2}/L$ and we find a reasonably good agreement. This suggests that higher powers of $1/\sqrt{L}$ still play an important role in the system sizes accessible to us through full exact diagonalization.

In this section, we study the off-diagonal matrix elements of observables in the eigenstates of Hamiltonian (2). We focus on their distributions, scaling properties, and functional dependence of $|f_{O}(\bar{E}\simeq 0,\omega)|^{2}$ on $\omega$, for eigenstates that are in the even-$Z_{2}$ even-$P$ sector within the $M^{z}=0$ and $k=0$ sector (see Sec. II).

We studied the bipartite von Neumann entanglement entropy and matrix elements of local operators in highly excited eigenstates of interacting integrable (the spin-1/2 XXZ chain) and nonintegrable models.

For the average entanglement entropy over all eigenstates in the zero magnetization sector, we found that the leading term is extensive at interacting integrable points with a coefficient of the volume-law that is smaller (for nonvanishing ratios $L_{\rm A}/L$) than the universal $\ln 2$ coefficient in quantum chaotic models. Finite-size scaling analyses suggested that the coefficient at $L_{\rm A}/L=1/2$, and for arbitrary ratios $L_{\rm A}/L$ (not reported), is (almost) independent of the XXZ chain anisotropy parameter $\Delta$, and that it is very close or equal to that of translationally invariant free fermionic Hamiltonians. Since the average entanglement entropy over all eigenstates is dominated by eigenstates at “infinite temperature”, the presence, or lack thereof, of interactions (and their values) at integrability may play no role in the leading extensive term. What may be essential is that the system is integrable so that it has an underlying quasiparticle description. Hence, we find it plausible that all translationally invariant (two-state per site) integrable models (noninteracting and interacting) have the same average entanglement entropy for any given ratio $L_{\rm A}/L$. If this is the case, then one could think of two universality classes for the average entanglement entropy of all “q-bit” based physical Hamiltonians, (translationally invariant) free fermions characterizing integrable models, and random matrices characterizing nonintegrable ones.

For the diagonal matrix elements of observables at the center of the spectrum and at interacting integrable points, we showed evidence that the support does not vanish with increasing system size and that the average eigenstate-to-eigenstate fluctuation vanishes as a power law in system size. At nonintegrable points, however, both the support and the average eigenstate-to-eigenstate fluctuation vanish exponentially with increasing system size. For the off-diagonal matrix elements with $\bar{E}=(E_{\alpha}+E_{\beta})/2$ at the center of the spectrum, we showed that at interacting integrable points they follow a distribution that is close to (but not quite) log-normal, and that their variance is a well-defined function of $\omega=E_{\alpha}-E_{\beta}$ whose magnitude scales as $1/D$, where $D$ is the Hilbert space dimension. The latter is a known property of the off-diagonal matrix elements of observables in nonintegrable models, which, however, exhibit a Gaussian distribution. We also studied the smooth function $|f_{O}(\bar{E}\simeq 0,\omega)|^{2}$ that characterizes the variance and contrasted its behavior at interacting integrable and nonintegrable points. It was recently argued that this function can be measured in experiments with periodically driven systems, both nonintegrable and interacting integrable ones, by studying how heating rates change when changing the frequency of the drive 1907.04261 . An interesting open question is whether the Bethe ansatz can be used to analytically learn about the smooth function $|f_{O}(\bar{E},\omega)|^{2}$ in interacting integrable systems. This would pave the way for an analytic understanding of the effect of interactions in the matrix elements of observables in many-body quantum systems.