Finding the Dynamics of an Integrable Quantum Many-Body System via Machine Learning

We study the dynamics of the Gaudin magnet ("central-spin model") using machine-learning methods. This model is of practical importance, e.g., for studying non-Markovian decoherence dynamics of a central spin interacting with a large bath of environmental spins and for studies of nonequilibrium superconductivity. The Gaudin magnet is also integrable, admitting many conserved quantities: For $N$ spins, the model Hamiltonian can be written as the sum of $N$ independent commuting operators. Despite this high degree of symmetry, a general closed-form analytic solution for the dynamics of this many-body problem remains elusive. Machine-learning methods may be well suited to exploiting the high degree of symmetry in integrable problems, even when an explicit analytic solution is not obvious. Motivated in part by this intuition, we use a neural-network representation (restricted Boltzmann machine) for each variational eigenstate of the model Hamiltonian. We then obtain accurate representations of the ground state and of the low-lying excited states of the Gaudin-magnet Hamiltonian through a variational Monte Carlo calculation. From the low-lying eigenstates, we find the non-perturbative dynamic transverse spin susceptibility, describing the linear response of a central spin to a time-varying transverse magnetic field in the presence of a spin bath. Having an efficient description of this susceptibility opens the door to improved characterization and quantum control procedures for qubits interacting with an environment of quantum two-level systems. These systems include electron-spin and hole-spin qubits interacting with environmental nuclear spins via hyperfine interactions or qubits with charge or flux degrees of freedom interacting with coherent charge or paramagnetic impurities.


I. INTRODUCTION
Predicting the dynamics of quantum many-body systems is crucial for understanding many important physical phenomena.For example, the dynamics of the Fermi-Hubbard model can advance our understanding of superconductivity and quantum magnetism in correlated materials [1,2].However, brute-force numerical approaches such as exact diagonalization on a classical computer have an exponential cost in time and/or memory and can therefore only be used to simulate small quantum systems.To tackle the problem of quantum many-body simulation, a large-scale fault-tolerant quantum computer could be used to efficiently run quantum simulations with predictable bounded errors [3,4].The hardware challenge behind building a useful quantum computer is, however, significant.General-purpose quantum simulation on a quantum computer may not be feasible until far in the future.
Despite the limitations of classical computers in simulating quantum systems, there are nevertheless many classical algorithms that are effective in special cases.A notable example is the variational Monte Carlo (VMC) method, which finds an approximate ground state upon minimizing the estimated energy with respect to a variational ansatz [5].Recently, Carleo et al. have used a neural-network ansatz in VMC to calculate ground states of many-body Hamiltonians with a better accuracy than other existing methods [6].A number of follow-up works on neural-network quantum states also strongly suggest that neural networks are a promising model for representing a large subset of quantum states and for approximating the ground states of many important many-body Hamiltonians [7][8][9][10][11][12].Using this type of ansatz, it is also possible to leverage well-studied optimization strategies from the deep-learning community.
A number of recent works have extended neural-network methods from the realm of static properties to many-body quantum dynamics [13][14][15][16][17].These have demonstrated success in some test cases (often where an analytic solution is already available), but they are unlikely to be successful for any arbitrary problem [18][19][20].Finding a subclass of problems that have known practical applications, that are nontrivial (where no closed-form analytic solution is available), and that have a good chance of admitting an efficient neural-network solution is an important open question.Here, we apply machine-learning methods to determine the low-energy spectrum and non-perturbative low-frequency dynamic response for one such model: the Gaudin magnet (Figure 1).This model applies directly to the many-body decoherence dynamics of a "centralspin" problem (describing, e.g., an electron spin qubit interacting with a bath of nuclear spins) [21][22][23][24][25][26][27][28][29][30][31].A linear combination of commuting Gaudin magnets can also result in the Bardeen-Cooper-Schrieffer (BCS) pairing Hamiltonian describing non-equilibrium superconductivity, with fermionic creation/annihilation operators represented in terms of Anderson pseudospins [32][33][34][35].A complete solution for the dynamics of the Gaudin magnet thus implies a solution for BCS pairing dynamics (through the composition of a product of commuting unitary time-evolution operators, one for each Gaudin-magnet problem).The Gaudin magnet is an integrable system, with a large number of conserved quantities.This suggests that drastic simplifications are possible.Although these simplifications are difficult to realize through direct analytic study, a well-designed machine-learning procedure may be able to exploit the high degree of symmetry in this problem.
For integrable N -particle systems, the O(e N ) set of linear equations that would normally describe the Schrödinger equation in a complete basis can be recast into a set of O(N ) nonlinear Bethe ansatz equations [36].The Gaudin magnet is one such system, and its integrability has been exploited in a number of works to obtain eigenstates and observable dynamics [25,[37][38][39].Despite the reduction to a much smaller system of Bethe ansatz equations, solving these nonlinear equations is still nontrivial in general.For example, a hybrid method based on a combination of the algebraic Bethe ansatz and direct Monte Carlo sampling proposed by Faribault and Schuricht is still limited to modest system sizes (N ≲ 50) [37].
Knowing that the eigenstates of the Gaudin magnet can be derived from O(N ) nonlinear equations, we view a neural-network quantum state ansatz as a promising candidate for obtaining the low-lying eigenstates and dynamics of this system.Neural networks such as the restricted Boltzmann machine (RBM) [40] can be used to describe quantum states using a number of network parameters that grows only polynomially in the system size N , but they can nevertheless be used to describe complex quantum states through the nonlinearity induced by tracing over hidden layers.We speculate that neuralnetwork methods may be able to calculate the low-lying eigenstates by learning the symmetries of the integrable Gaudin magnet and other integrable models.
In this work, we calculate the low-lying eigenstates of the Gaudin magnet (central-spin model).To do this, we use a penalty-based variational algorithm with an RBM neural-network ansatz.We also compute the nonperturbative transverse spin susceptibility for the central spin, giving the linear response of the central spin to a time-varying magnetic field, accounting for highly accurate representations of the many-body eigenstates of the central spin interacting with environmental spins.We find these approximate eigenstates (and the associated transverse spin susceptibility) numerically in a regime where conventional perturbation theory fails.
Due to the widespread practical importance of the Gaudin magnet, there have been many other proposed strategies to determine the non-perturbative dynamics of this model.The most natural approach is to restrict to uniform coupling coefficients (the box model), in which case the problem can be diagonalized exactly analytically [41,42].The box model can accurately describe dynamics for the original problem (with a broad distribution of coupling coefficients) at short times (where semiclassical methods also apply).The box model is, however, too restrictive to be used to understand long-time quantum dynamics in most practical contexts, except where uniform coupling coefficients can be engineered directly [43].Several authors have considered an extension to a socalled wedding-cake model, composed of layered rings of environmental spins that behave as single larger spins having equal coupling constants [44][45][46].The wedding-cake model allows for efficient exact numerical calculation while sacrificing the original distribution of coupling coefficients.Recent work has analyzed the case where the distribution of piecewise-constant coupling coefficients is chosen to reproduce moments of the original coupling-constant distribution [31].This method accurately reproduces correlation functions at infinite temperature and at full polarization, but it is not clear how this wedding-cake method will perform for low-temperature and non-equilibrium initial states, as considered here.Finally, numerical methods based on the time-dependent density-matrix renormalization group [29,30] and analytic methods based on a partial resummation [28] are both limited to short time scales (set by the largest coupling coefficient), while we will be interested in dynamics at and beyond this time scale in the present work.
The rest of this paper is organized as follows: In Section II, we introduce the Gaudin-magnet model Hamiltonian and describe a particular physically relevant set of coupling constants.In Section III, we describe the variational ansatz and variational algorithms used to find accurate representations of the ground state and several low-lying excited states.Section IV describes the dynamic transverse spin susceptibility for this problem and presents sample calculations for the non-perturbative spectral function and linear response to a low-frequency transverse field.Section V summarizes the potential advantages and shortcomings of this approach, along with conclusions and possible future directions.A technical derivation of the gradient expression used for gradient-based optimization is given in the Appendix.

II. MODEL
The central-spin model (Gaudin magnet) consists of one central spin coupled to N environment spins (Figure 1) [22,47,48].The central spin is coupled to the k th environment spin through a Heisenberg interaction with coupling coefficient A k .The central spin alone is additionally subject to a constant external field B. The Hamiltonian for the Gaudin magnet is thus where we have set ℏ = 1, S α j = 1 2 σ α j (α = x, y, z) is the spin-1/2 operator for spin j (j = 0 for the central spin and j = k = 1, 2, . . ., N for the environment spins).The operator σ α j (α = x, y, z) is the Pauli operator for spin j.In numerical evaluations, we choose the coupling coefficients A k to have an exponential distribution, The total number of environment spins in the model is N while N 0 controls the scale where the coupling coefficients decay exponentially.We will typically be interested in systems where N > N 0 ; in numerical evaluations below, we choose N 0 = (N + 1)/2.An exponentially decaying coupling strength corresponds, e.g., to the distribution of hyperfine couplings expected for an electron spin (central spin S 0 ) interacting with nuclear spins (environment spins S k ) in a two-dimensional quantum dot with parabolic confinement [23].In this example, N 0 is the number of nuclear spins within a quantum-dot Bohr radius, while N → ∞ is the number of nuclear spins in the entire crystal.The low-lying spectrum of this model is shown for N = 5 [N 0 = (N + 1)/2 = 3] in Figure 2 for a range of B. In dynamics calculations below, we set B = 0.35A/N 0 (blue dashed line in Figure 2), chosen to avoid degeneracies up to the fourth excited state.For B ≫ A, we can directly apply perturbation theory to give the eigenstates of H as the simultaneous eigenstates of all S z j and in this case the spectrum is trivial.However, for B ≲ A, this simple perturbation theory does not generally apply and more advanced methods are required to understand the detailed spin dynamics arising from this model.

III. VARIATIONAL ALGORITHMS
In this section, we introduce the RBM variational ansatz used to approximate the eigenstates of the Hamiltonian H.We then describe the variational algorithms used to optimize the RBM network parameters and we explain strategies used to select accurate representations of the states.

A. Variational Ansatz
We represent a quantum state vector with an RBM ansatz [6].An RBM is a two-layer network model (Figure 3), fully defined by a set of parameters W = {a j , b i , w ji }, where here we choose the parameters to be complex-valued.The RBM representation of a quantum state |ψ⟩ in the computational basis (σ z j -eigenbasis) where the "∝" symbol indicates that the RBM representation of a quantum state is generally unnormalized and where σ z j |σ⟩ = σ j |σ⟩ with σ j = ±1.The unnormalized amplitude Ψ W (σ) of a particular spin configuration is given by where a 0 , ..., a N and b 1 , ..., b M are the bias weights for the visible and hidden nodes, respectively.The visible nodes are associated with Ising variables σ j = ±1 and the hidden nodes are assigned values h i = ±1.The parameters w ji assign an independent weight to each link between a visible node j and a hidden node i.For simplicity, we set the number of hidden nodes to be equal to the number of visible nodes, N +1 M = 1.In general, this need not be Visible Layer Hidden Layer A restricted Boltzmann machine (RBM) architecture applied to the Gaudin magnet.The visible layer, shown on the left, is composed of Ising spin variables σ = {σ0, σ1, ..., σN }, describing the eigenstates of Pauli operators σ z j with eigenvalues σj = ±1 at each node j.The hidden layer, on the right, is made up of the Ising variables h1, h2, ..., hM .The parameters a0, ..., aN and b1, ..., bM are the bias weights of the visible and hidden nodes, respectively.The parameters w01, ..., wNM are weights connecting the visible and hidden nodes.
the case and increasing the number of hidden nodes for a fixed number of visible nodes can increase the expressivity of the RBM.Marginalizing over (tracing out) the hiddennode variables h i yields a nonlinear function that only takes the variables σ as its input and this function is fully determined by the network parameters W.

B. Ground State
To find the ground state, we minimize the energy expectation value [6], Calculating the exact expectation value would generally involve a sum over 2 N +1 basis states σ.To avoid an exponential cost in the computation time, we instead evaluate this quantity approximately ⟨...⟩ σ via samples σ obtained from the Metropolis algorithm [49], where the samples are drawn from the probability distribution, We use the stochastic reconfiguration method [50] to minimize the energy subject to variations in the parameters W. Stochastic reconfiguration exploits information encoded in the quantum geometric tensor (the covariance of the logarithmic derivative of Ψ W ) to precondition the gradient used in stochastic gradient descent (SGD).The energy gradient (force vector) is defined as where the logarithmic derivative of the RBM state vector is defined as The partial derivative in Equation ( 7) is taken with respect to the complex conjugate of the variational parameter, since we are looking for the steepest descent of a real-valued function parameterized by complex-valued variables [51].The quantum geometric tensor is defined as At the beginning of an optimization run, a random initial ansatz |W(0)⟩ is generated by setting the real and imaginary parts of all variational parameters independently to values obtained from a symmetric Gaussian distribution with standard deviation σ W = 0.25.At each iteration l = 1, . . ., L, the parameters are updated according to where γ L is the learning rate.Here, S λ = S + λI includes a diagonal shift λ, a regularization that helps with the matrix inversion [6].The specific choice of λ is given in Sec.III D, below, along with an explanation of other hyperparameters.See Figure 4 for a typical successful search for the ground state.The figure shows the estimated energy after each iteration with L = 8000.This constitutes a single optimization run.At the end of each run, an accurate estimate is obtained for the energy of the final state (for l = L) and this is stored.In practice, several runs p are performed (p = 1, 2, . . ., P ) and of those P runs, the run corresponding to the lowest final energy is selected.See Algorithm 1, including the procedure for finding excited states, which we now describe in detail.

C. Excited States
Most of the work on neural-network quantum states has focused on the ground state, but a few works have also found approximate excited states by projecting out the previously determined eigenstates [11,50].In general, exactly projecting out the ground state (or any previously determined lower-energy states [52,53]) may introduce Input : The n lowest-energy eigenstates in the form of RBMs, { W 0 , W 1 , W 2 , ... W n−1 }, and penalty coefficients {β0, β1, β2, ...βn−1}.
Output : An RBM representation for the n th excited state, |W n ⟩. an exponential cost, associated with the cost of exactly calculating an overlap between the current trial state and the previously determined states.Here, as in the approach of Choo et al. [11], we instead estimate the relevant overlaps with approximate Monte-Carlo sampling.However, rather than directly projecting out the previously calculated eigenstates, we adopt a penalty method [53] to calculate the excited states.In this method, the n th excited state is calculated as the ground state of the modified Hamiltonian, Here, the terms β j > 0 are penalty coefficients and W j is the approximate j th excited state, defined in terms of variational parameters W j .The ground state is indexed by j = 0. Explicitly, we replace the energy expectation value with the new loss function, After optimization, Ẽn (W) is the approximate energy of the n th excited state of the original Hamiltonian H.
The individual penalty coefficients β j must be sufficiently large to ensure that the previously calculated eigenstates (j = 0, 1, . . ., n − 1) are all raised in energy (in H n ) above E n .However, very large values for β j will be detrimental to optimization performance, a well-known phenomenon in penalty-based constrained optimization [54].Here, our goal is to find a subset of energy eigenstates describing the low-frequency response of the Gaudin magnet to external perturbations.We therefore calculate eigenvalues only up to a maximum cutoff ω max above the ground-state energy.For simplicity, in numerical evaluations we have therefore set all penalty coefficients to a j-independent value, For the specific spin-dynamics results presented below for the Gaudin magnet, we choose ω max = 0.15 A/N 0 , which is sufficient to characterize the first five levels with the chosen parameters: N = 5, N 0 = 3, B = 0.35A/N 0 , and with the exponential distribution of coupling constants A k given in Equation ( 2).The minimization strategy for excited states differs from that of the ground state, since the loss function is modified.While the definition of the quantum geometric tensor S ik remains unchanged, the gradient vector must be modified to (see the Appendix for details): where the modified energy gradient F n i is evaluated approximately through Metropolis sampling.The samples {σ} are collected based on the distribution π(σ; W) [Equation (6)], whereas {σ j } are sampled accounting for the variational parameters W j , according to π(σ; W j ).See Figure 5 for a typical successful run resulting in the first excited state.In Figure 6 we show a histogram of outcomes for the first five energy levels E 0 , E 1 , . . ., E 4 after 50 optimization runs.

D. Hyperparameters and Postselection
Here we summarize the chosen model hyperparameters.The learning rate was set to γ L = 0.02 for all runs.We perform Metropolis sampling with N σ = 5000 samples to estimate the energy loss function [E(W) for the ground state, Ẽn (W) for the n th excited state], the gradient vector F i , and the quantum geometric tensor S ik .As described following Equation (10), a diagonal shift λ = 0.01 is added to the quantum geometric tensor before inversion.The total number of iterations is fixed to L = 8000.
To further decrease the error of the final approximate ground state, we use the following postselection strategy: 50 independent runs (each consisting of 8000 iterations) were performed with different random seeds, storing the sampled energy and the variational parameters from the last iteration of each run.We then used 10 4 Metropolis samples (twice as many samples as in the optimization runs) to obtain a more accurate estimate of the energy in each of the 50 runs.We then select the state with the lowest energy (out of the 50 states found from the 50 runs) to be the best approximate ground state.The postselection process for each excited state was the same, after replacing E(W) (the estimated energy alone) with Ẽn (W) (estimated energy plus penalty terms).
See Algorithm 1 for a summary of the n th excited state calculation with the number of independent runs P = 50 and the number of iterations L = 8000 for both the ground state and the excited states.With the chosen hyperparameters, we found good approximations for the lowest five eigenstates.By comparing to exact diagonalization, we find the error in each of the energy eigenvalues is less than 10 −3 A/N 0 .For all five eigenstates, the state infidelity was bounded by 1 − |⟨ψ exact |ψ⟩| 2 < 2 × 10 −3 .

E. Time Scaling
The average runtime (averaged over 50 runs) is shown in Figure 7 for a determination of the ground state of the Gaudin magnet for N = 1 to N = 11.The runtime for the RBM-based variational algorithm (green triangles) is compared to a brute-force exact diagonalization (blue triangles).Exact diagonalization was performed using the QR algorithm for Hermitian matrices, implemented in NumPy [55].While the average runtime of brute-force exact diagonalization scales exponentially, the RBM-based variational algorithm with O(N 2 ) variational parameters (based on the RBM architecture shown in Figure 3) has the potential to scale polynomially.Figure 7 suggests that the RBM-based variational algorithm will have a shorter runtime than brute-force exact diagonalization as the system size grows beyond N ≳ 12.For this comparison, we have fixed the number of samples and the number of iterations.We cannot rule out the possibility that, as the system size grows, the number of samples and iterations required to reach a target error may increase exponentially.
As we show in Section IV below, the ground state of this problem can actually be found through a refined exact diagonalization procedure in a reduced subspace with cost O(N ), since the lowest-energy state always lies within a manifold of at most one spin flipped.The comparison here to brute-force exact diagonalization of the entire Hamiltonian is only illustrative and likely not of practical relevance for this particular problem (finding the ground state).However, finding a general excited state (or collection of excited states to describe the system at For each system size, 50 runs were performed for each of the two algorithms to obtain the average runtime.All numerical runs were performed on the Graham cluster located at the University of Waterloo, Waterloo, ON Canada, with a single CPU and 1024 MB of memory for each run.We also fix both the number of samples and the number of iterations to be 1000. finite temperature) will typically have exponential cost for direct numerical diagonalization.

IV. DYNAMICAL RESPONSE
Given an accurate representation for both the low-lying spectrum of a many-body quantum system, as well as the low-lying eigenstates, dynamical response functions can be found.In the presence of a time-varying transverse field B ⊥ (t) = B x (t)x + B y (t)ŷ acting on the central spin of the Gaudin magnet, the total Hamiltonian becomes with time-dependent perturbation Here, we have introduced We assume an initial state ρ(t 0 ) that is stationary with respect to the Gaudin-magnet Hamiltonian H: This is true for a thermal state or any other statisitical mixture of H eigenstates.If the initial state is a mixture of non-degenerate Gaudin-magnet eigenstates, the linear response of Here, S − 0 (t) evolves under the action of the total Hamiltonian H tot (t), while Ŝ− 0 (t) = e iHt S − 0 e −iHt evolves in the interaction picture under H.The average is understood to be taken with respect to the initial state, We have also used the following identities, valid for a mixture of non-degenerate H eigenstates: and These identities follow directly from the fact that the Gaudin-magnet Hamiltonian commutes with the total z-component of spin: while the operators S ± 0 couple sectors of different J z .Rewriting Equation ( 21) in terms of Fourier transform variables and taking t 0 → −∞ gives where the transverse spin susceptiblity is Here, θ(t) is the Heaviside step function and 0 + is a positive infinitesimal.In physical applications, the Gaudinmagnet Hamiltonian will only be an approximate description.Corrections to this description (coupling to a continuum) will lead to a finite decay rate.We account for this effect with the phenomenological replacement We now specialize to the case where the Gaudin magnet is prepared in a non-degenerate ground state, ρ(t 0 ) = |0⟩ ⟨0|, where H |0⟩ = E 0 |0⟩.Expanding in a complete set of energy eigenstates then gives where the j th excitation energy is The excitations that can be generated through centralspin spin flips will result in peaks of the central-spin spectral function, where we have introduced a lineshape function For γ → 0 + , the lineshape function approaches a Dirac delta function, δ γ (ω) → δ(ω).Positive-frequency (ω > 0) contribuions to A 0 (ω) arise from excitations generated through the action of S + 0 on the ground state, while negative-frequency (ω < 0) contributions arise from excitations produced through S − 0 .Without loss of generality, we take B ≥ 0 to define the +z direction of spin.A finite value, B ̸ = 0, favors a ground state with the central spin pointing down, |⇓⟩.When the couplings are negative, A k < 0, the environment spins will then favor an orientation along the central spin in the ground state.This fully polarized state is an exact eigenstate of H, while for A k > 0, the ground state will have mixed character: where For A k < 0, we have This leads to peaks in the spectral function only at positive frequencies ω > 0 for A k < 0, with amplitudes determined by the overlaps ⟨j| ⇑↓↓ A time-dependent transverse field By(t) (subplot a)) drives the central-spin operator ⟨S x 0 (t)⟩ (subplot b)) for a Gaudin magnet with N = 5 environment spins and with exponential distribution of coupling coefficients A k given in Equation ( 2).In subplot b), the red curve shows the approximate result in linear response calculated from an RBM neural-network ansatz and the dotted black curve is the exact result, obtained from exact numerical integration of the Schrödinger equation.The embedded plot in subplot a) shows the spectral weight in the Fourier transform By(ω) relative to the central-spin spectral function A0(ω).
detailed knowledge of the many-body eigenstates |j⟩ as well as the excitation energies ∆ j to accurately determine the spectral function.
For positive coupling, A k > 0, we have instead In contrast with the case of A k < 0, the case of A k > 0 will lead to peaks at both negative and positive frequencies.Since the fully polarized state in Equation ( 39) is an exact eigensate of H with eigenvalue (B + k A k /2)/2, there will be a single peak at ω ≃ B + k A k /2 > 0 controlled by the amplitude β 0 , while there will generally be a family of peaks at low negative frequencies ω ≃ −∆ j , with amplitudes controlled by 8 for an example).In what follows, we focus on this case of A k > 0.
Based on the discussion here, we see that the zerotemperature central-spin spectral function can, in fact, be calculated numerically exactly with only a polynomial cost.For A k < 0, the problem amounts to calculating the excitation energies ∆ j and overlaps ⟨j| ⇑↓↓ • • • ↓⟩ after finding the eigenstates |j⟩ through exact diagonalization in the O(N )-dimensional subspace of one spin flipped.For the opposite case of A k > 0, in addition to calculating the ground state in the subspace of one spin flipped, the associated excitation energies ∆ j and overlaps ⟨j| ⇓↑↑ • • • ↓ k • • • ↑⟩ can then be found through exact diagonalization in the O(N 2 )-dimensional subspace of two spins flipped.In this (zero-temperature) limit, there is likely no practical need for a machine-learning (RBM) approach to dynamics.However, for a finite-temperature thermal state or for initial conditions corresponding to a non-equilibrium configuration described by many spin flips, an efficient alternative to exact diagonalization becomes important.With n flipped spins, the number of states in the subspace ( N n ∼ (N/n) n for N ≫ n ≫ 1) grows exponentially with n.While exact-diagonalization necessarily becomes computationally expensive in this limit, the machine-learning (RBM) approach generalizes directly to different initial conditions and has the potential to scale favorably.
A possible alternative to the RBM method presented here is perturbation theory, but as we show in the following section, perturbation theory fails for this problem in the limit of small-to-moderate B (B ≲ A).

A. Perturbation theory
For a sufficiently large B, it is possible to find perturbative approximations for the excitation energies and amplitudes, by writing H = H 0 + V ff with An expansion in the flip-flop terms V ff then gives: In the second line above, we have used k A k ≃ A for large N and in the third line we have used the specific coupling coefficients given in Equation ( 2) to evaluate the sum, This contribution to the spectral function will dominate (|β 0 | 2 ≃ 1) and perturbation theory is valid for calculating this quantity for any non-negative B whenever N 0 ≫ 1.The parameter |β 0 | 2 is the same quantity that gives rise to the long-time saturation in spin coherence describing non-Markovian partial coherence decay for a central-spin system when the interaction is introduced suddenly [23].If only these highfrequency features were of interest, there would typically be no need for a non-perturbative approach (although higher-order corrections in V ff will broaden this sharp feature [56]).The low-frequency contributions to the susceptibility driving the long-time dynamics have a much more stringent condition for perturbation theory to be valid.In particular, due to the possibility for constructive interference of leading-order contributions in perturbation theory, an accurate description of these terms generally requires [23]: This condition is only satisfied for B ≫ A/2, independent of N 0 .Outside of this regime, this naïve perturbation theory will fail to accurately describe dynamics.Our focus is now on describing the low-frequency contributions with a non-perturbative method.We emphasize that the transverse spin susceptibility χ −+ (ω) gives only the linear response of the central spin to a weak driving field (leading order in B ⊥ ), but an accurate description of the resulting spin dynamics for B ≲ A requires a nonperturbative calculation to all orders in V ff , even in this limit of a weak driving field.

B. Nonperturbative many-body spin susceptibility
In this section, we directly apply the RBM ansatz to the problem of finding the zero-temperature dynamical response.While brute-force exact diagonalization is likely to fail for a more general (finite-temperature or nonequilibrium) initial condition due to the exponential growth of the relevant Hilbert space, the machine-learning procedure presented here is directly generalizable to nontrivial initial conditions and has the potential to remain efficient in this regime.In this paper, we simply present a first example involving a small number (N = 5) of environment spins and we assume a zero-temperature ground state.We leave it to future work to find scaling of the RBM procedure for more general initial conditions.
To find an accurate approximation for the nonperturbative spin susceptibility, we construct RBM representations of the lowest-energy j max + 1 eigenstates of H: where j max is defined (for a fixed maximum frequency ω max ) by The relevant matrix elements required to find the susceptibility are then approximated through optimized RBM representations: where we rewrite each of the quantum averages in terms of a weighted sum: Here, we have introduced the notation where π(σ, W j ) is given by Equation ( 6).On the righthand side of Equation ( 51), we have approximated the average by a sample average over N σj samples σj that are found in practice through Metropolis sampling.
Once we have the subset of approximate eigenstates in the form of an RBM, we use Metropolis sampling to estimate the matrix elements given in Equation ( 49) and (50).We used N σ = 5 × 10 6 samples for the estimate, resulting in a relative statistical error < 10 −2 for the relevant matrix elements.These approximate matrix elements were then combined with the estimated excitation energies ∆ j to construct the central-spin spectral function, Equation (31).See Figure 8 for a comparison of the approximate central-spin spectral function with the result from exact diagonalization.The exact solution accounts for all eigenstates, while the approximate spectral function was calculated using only the lowest five approximate eigenstates.As an illustrative application of the linear-response calculation above, here we consider the real-time dynamics of the central spin in the presence of a slowly-varying transverse field.See Algorithm 2 for the steps taken to perform this computation.The spin susceptibility defined in Equation (28) directly gives this result.In particular, for a total Hamiltonian H total = H +V (t) = H +B y (t)•S y 0 , and for a t = 0 initial ground state of H, the linear response is where Finally, we choose B y (t) to have nonvanishing spectral weight at two of the excitation frequencies (inset of Figure 9a)).For an illustrative example, we make the specific choice where g(x, σ x ) is a normalized Gaussian envelope: Here, B 1 = B 2 = 5 A N0 , t = 200 N0 A , τ 1 = 100 N0 A and τ 2 = 50 N0 A .The functional form of B y (t) is shown in Figure 9a) and the resulting linear response ⟨S x (t)⟩ is shown in Figure 9b).The neural-network (RBM) estimate (solid red curve) accurately reproduces the exact dynamics (black dashed curve) in this linear response regime.The exact result (black dashed line in Figure 9b)) was obtained by integrating the exact Schrödinger equation directly via the Python package Qutip [57], without any linearresponse assumption.
The amplitude of the field B y (t) is deep in the linearresponse regime, so we expect that deviations between the exact and approximate curves shown in Figure 9b) are primarily due to errors in the approximate variational (RBM) eigenstates and in the calculations of matrix elements through Metropolis sampling.There may also be some small correction due to the finite-frequency cutoff ω max = 0.15(A/N 0 ) taken in the approximate RBM calculation.

V. CONCLUSIONS
In this paper, we have applied a neural-network method to calculate approximate low-lying eigenstates and linear response dynamics of an integrable many-body Hamiltonian: the Gaudin magnet.Having an efficient solution for this many-body system could be an important tool in characterizing decoherence sources for qubits that interact with uncontrolled two-level systems.These systems include qubits based on spin, charge, flux, etc., interacting with nuclear spins, coherent charge traps, or paramagnetic impurities.The very same model (the Gaudin magnet) can be used to study the dynamics of nonequilibrium superconductivity since a linear combination of commuting Gaudin-magnet Hamiltonians can be used to represent the Richardson model (BCS Hamiltonian) that drives superconducting pairing dynamics.
This work was motivated in large part by the fact that the Gaudin magnet is a quantum integrable model, admitting a large number of conserved quantities.For such integrable models consisting of N particles, it is well known how to recast the system of O(e N ) linear equations arising from the Schrödinger equation in terms of O(N ) nonlinear equations, via the algebraic Bethe ansatz.It is not generally known how to efficiently solve these nonlinear equations.This understanding has led us to seek a solution to this dynamics problem via a variational neural-network ansatz (the restricted Boltzmann machine).In contrast with brute-force exact diagonalization in the entire Hilbert space, which is guaranteed to have exponential cost, the neural-network approach taken here relies on a limited (polynomial) number of network parameters and approximate sampling from a subset of states.The neural-network approach therefore has the potential to be efficient.The tradeoff is that this variational neural-network approach has no guarantee of convergence or accuracy, but the first results presented here are promising in this respect.
We have presented a systematic procedure that could be used to explore the dynamics of the Gaudin magnet or other interesting models.We have further given some illustrative examples based on very small systems (one central spin and N = 5 environment spins).In this analysis, we have not provided any guarantees of convergence, accuracy, or efficiency (although we have provided direct comparisons with exact results suggesting each of these can be achieved for small systems).In future work it will be important to explore these questions systematically by studying systems of progressively larger size and under conditions that push past the boundaries of what can be done via exact diagonalization.Further improvements can be made by adopting contemporary neural-network quantum state architectures such as autoregressive neural-network models [7,58,59].In this architecture, the neural-network quantum state is normalized by construction and sampling can be done more efficiently.
It remains an interesting open question whether a variational machine-learning approach, as explored here, can easily "learn" the symmetries of a many-body problem and exploit these to find an efficient solution.In this Appendix, we provide a detailed derivation of the gradient expression from Equation (14).The variational gradient used to find excited states is derived by differentiating the Monte-Carlo estimated energy including a penalty term (Equation ( 12)).The derivative is evaluated with respect to the model parameter W. The first contribution in the gradient comes from the estimated energy alone and has been widely used in existing works [6,7,50].The derivation here focuses on the second contribution to the gradient from the penalty term.= 0 [60].Incorporating the normalization terms in the denominator, the first term becomes The second term in Equation (A1) can be expanded as (A4) Summing up the two terms in Equation (A3) and (A4), and then generalizing to multiple penalty terms, we recover the gradient expression in Equation (12).For the normalized neural-network quantum states in, e.g., Ref. [7], the gradient expression would only contain the first term in Equation (A3).

FIG. 1 .
FIG. 1. Diagrammatic representation of the Gaudin magnet, where S k are the individual spin-1/2 operators, A k is the coupling between the k th spin and the central spin, and B is the external field.
FIG. 4. A typical successful run to find the Gaudin-magnet ground state with N = 5, B = 0.35A/N0, N0 = 3, and A k distributed as in Equation (2).The dashed horizontal line indicates the exact energy E1 − E0 corresponding to the first excited state.Note the logarithmic scale.The subplot displays a zoom-in view near convergence using a linear scale for the energy axis, where the solid horizontal line corresponds to the exact ground state energy.

FIG. 5 .FIG. 6 .
FIG. 5. A typical successful run resulting in the first excited state of the Gaudin magnet with N = 5.Horizontal lines indicate the exact excitation energies Ej − E0 for the first eight excited states (j = 1, 2, . . ., 8).The inset displays a zoom-in view near convergence with the excitation energy in a linear scale.Convergence slows at an energy close to E8, which strongly suggests the existence of a local minimum at that particular point of the energy landscape.

FIG. 7 .
FIG.7.Runtime scaling for the ground state of the Gaudinmagnet Hamiltonian.The runtime for finding the ground state via the RBM-based variational algorithm (blue symbols) is compared to brute-force exact diagonalization (green symbols).For each system size, 50 runs were performed for each of the two algorithms to obtain the average runtime.All numerical runs were performed on the Graham cluster located at the University of Waterloo, Waterloo, ON Canada, with a single CPU and 1024 MB of memory for each run.We also fix both the number of samples and the number of iterations to be 1000.
Appendix A: Derivation of the gradient