A digital quantum simulation of the Agassi model

A digital quantum simulation of the Agassi model from nuclear physics is proposed and analyzed. The proposal is worked out for the case with four different sites. Numerical simulations and analytical estimations are presented to illustrate the feasibility of this proposal with current technology. The proposed approach is fully scalable to a larger number of sites. The use of a quantum correlation function as a probe to explore the quantum phases by quantum simulating the time dynamics, with no need of computing the ground state, is also studied. Evidence is given showing that the amplitude of the time dynamics of a correlation function in this quantum simulation is linked to the different quantum phases of the system. This approach establishes an avenue for the digital quantum simulation of useful models in nuclear physics.


Introduction
During the past few decades the possibility of using controllable quantum systems to simulate other quantum systems has been explored extensively [1]. Different quantum platforms have been proposed to reproduce quantum models experimentally, including superconducting circuits, ion traps, cold atoms, quantum dots, as well as quantum photonics [1]. One of the emerging fields proposed for quantum simulations is the analysis of nuclear physics models. In particular, a cloud quantum computing of an atomic nucleus [2], quantum simulations of Schwinger-model dynamics [3,4,5,6], and quantum simulations of quantum field theories with trapped ions and superconducting circuits [7,8,9,10] have been proposed and sometimes experimentally realized. For a thor-ough review of this research field with updated references see Ref. [11]. However, a paradigmatic quantum nuclear system such as the Agassi model [12] has not been analyzed in the context of quantum simulations. Its relevance in Nuclear Physics, but also in a wide variety of fields, including many-body quantum systems and quantum phase transitions, as well as the difficulty to numerically compute the dynamics and static properties of large quantum systems, motivates the quantum simulation of the Agassi model.
The importance of quantum information science (QIS) in Nuclear Physics cannot be underestimated in view of the report of the U.S. Department of Energy (DOE) [13] where in its Research Opportunity II establishes "A broad theory program should be supported, which can, e.g., develop methods to address problems in NP using digital quantum computers and quantum simulators, utilize QIS concepts to better understand nuclear phenomena (such as the nuclear many-body problem and hadronization), and develop new QIS applications of importance to nuclear physics". The present work represents a first step in the direction of this general goal. The use of quantum computing for solving present Nuclear Physics problems has been reviewed in [14]. The Agassi model [12] is a simple but far from trivial quantum model that includes a combination of long range monopole-monopole and short range pairing interactions. It was first proposed in nuclear physics since it is an exactly solvable model that provides a schematic version of the pairing-plusquadrupole model that has been extensively used in nuclear structure [15]. From the quantum phase transition view, this model presents a rich quantumphase diagram for the ground state, containing several phase transition lines [16,17,18], and has been widely studied in a variety of fields. Apart from the symmetric phase, the model has two brokensymmetry phases: one superconducting, linked to the pairing interaction, and another parity-broken phase linked to the monopole-monopole interaction. The phase diagram of the Agassi model has been studied within a mean-field formalism [16,17,18]. As known, this kind of formalism is valid for the thermodynamic or large-N limit, where N is the number of sites. However, for mesoscopic systems, where finite-N effects are important, the corresponding phases and transitions are blurred and more detailed studies are needed for a clear understanding. In addition, beyond-mean field methods to calculate finite-N effects are difficult to apply for moderately small-N. For this purpose, quantum platforms could be used to mimic the Agassi model. On the other hand, tools from quantum information, as the quantum discord, have been recently employed to explore the phases in this model to gain insight about its structure [19].
In this paper, we propose and analyze the digital quantum simulation of the Agassi model [12]. Although we propose a fully digital scheme, for some useful comparisons we refer to trapped-ion platforms [20,21]. A quantum simulation of the Agassi model may enable one to carry out a full-fledged analysis of this model for a mesoscopic number of sites, in such situations where all classical methods will fail. For instance, apart from the mean field calculations, no finite-N corrections have been calculated for the model in its simplest version, even for the first correction term. In addition, the extended Agassi model, Ref. [17], includes extra terms producing up to 5 different phases with three control parameters. With our approach, the extension to possible scenarios with inhomogeneous couplings where mean field methods will fail is direct, allowing one for the scalable quantum simulation of nuclear physics models inaccessible to classical supercomputers.
In this work, we also study how to employ quantum correlation functions as a probe to explore the quantum phases in the system via a quantum simulation of the time dynamics, without needing to compute the ground state. Indeed, we give evidence, analysing the time dynamics of a correlation function, that its amplitude can be linked to the different quantum phases of the model. Thus, a measure of this time dynamics, that can be done routinely with present technology, will provide the system phase.

The Agassi model
The Agassi model [12] consists of N interacting fermions which occupy two levels, each of degeneracy Ω, where Ω is even, and j = Ω/2. Note that in the following, we will consider N = 2Ω. The lower level σ = −1 has negative parity, and the upper level σ = 1 has positive parity. The magnetic quantum number takes the values m = ±1, . . . , ±j (note that m = 0 is excluded). Thus, a single-particle state is labeled by (σ = ±, m). The model is an extension of the Lipkin-Meshkov-Glick [22] model introduced by D. Agassi as a toy model to test many-body theories and to explore the interplay between particle-hole and superfluid correlations. However, the appearance of the model in the literature is scarce. Davis and Heiss [16] derived the phase diagram of the model and the different collective excitations in the existing phases, using Hartree-Fock-Bogoliuvov (HFB), particle-hole RPA and QRPA approximations. The Agassi model was also used to test some cumbersome numerical methods such as the merging of Coupled Cluster with the symmetry restored HFB theory [23]. In [17,18], the authors extended the model by the introduction of new interaction terms that give rise to an extremely rich phase diagram. In [24], the model was used as a test-bed for a number conserving particle-hole RPA theory. Finally, in [19], the authors use the Agassi model to study the so called two-orbital quantum discord.
The Agassi Hamiltonian is

sumed. The operators in H are
where c † σ,m (c σ,m ) are fermion creation (annhilation) operators in the state |σ, m .
The Agassi model (1) has a phase diagram with three phases, namely, spherical, deformed, and superfluid, in other words, a symmetric phase and two broken-symmetry phases (deformed HF and superfluid BCS). It is customary to divide the Hamiltonian by ε and to define the scaled parameters χ and Σ as g = Σ 2j−1 and V = χ 2j−1 to correctly scale the parameters with the system size. The phase diagram of the Agassi model is sketched in Fig. 1.

Quantum simulation of the Agassi model
To simulate a quantum model, a mapping between the original Hamiltonian and one suited for the digital simulation is needed. Here an Agassi Hamiltonian with j = 1 is considered. This contains four different sites, to analyze a case that may be experimentally realized with current technology. To simplify the notation we relabel the original fermion operators as (7) and the corresponding relationships for the creation operators. The mapping is carried out through the Jordan-Wigner image of the above fermions and is written as, and the corresponding Hermitian conjugate one for the creation operator. σ x i , σ y i , σ z i are Pauli matrices at position i and σ ± i = 1 2 (σ x i ± iσ y i ). The symbol ⊗ stands for the tensor product. We consider in this work N = 4, then the model space is of dimension 2 4 = 16 and, therefore, each operator is given by a 16 × 16 matrix. It is worth to mention that the Jordan-Wigner transformation produces a non-local Hamiltonian.
The spin image of the building block operators is Finally, one can to write down the Agassi Hamiltonian (1) for the case of j = 1 as, where Note that H 1 and H 2 only depend on σ z and, therefore, any state with well defined σ z components will be its eigenstate. H 3 depends on g and V and it vanishes for g = −V . Moreover, one should consider that, 3 The term H 3 can be further decomposed in terms of tensor products of Pauli matrices, where the symbols ⊗ have been taken out to simplify the notation. It is worth noting that for this simple case, j = 1, the ion-mapped Hamiltonian (14) depends on just one effective control parameter, g + V (see Eq. (19)), and not on g and V separately, as in the thermodynamic limit of the model [16,17,18]. Therefore, it is only possible to distinguish for j = 1 between a symmetric phase (SP) that is obtained for g + V < 1 and a broken-symmetry phase (BSP) emerging for g + V > 1. In this simplest case, the phase diagram is of dimension 1 as shown in Fig. 5 (upper coloured  panel). The critical point in the transitional path between these two phases is g + V = 1.

Theoretical model for the implementation
In order to carry out a quantum simulation with the Agassi model, we propose to employ a digital protocol, via a Lie-Trotter-Suzuki decomposition [1]. The protocol will rely on expressing the quantum evolution operator U (t) = exp(−iHt) for the Hamiltonian H in Eq. (14) by means of a Trotterized dynamics, in terms of H 1,2,3 of Eqs. (15), (16), and (17), where the error produced will depend on the commutator [(H 1 + H 2 ), H 3 ] and scale as 1/n T , where n T denotes the number of Trotter steps.
Once the dynamics has been decomposed into the previous blocks, each of these can be implemented efficiently with trapped-ion systems. The operator exp(−iH 1 t) consists of single-qubit gates which are customary with trapped ions, to fidelities often above 99.99% [25]. The operator exp(−iH 2 t) is composed of two two-qubit gates which can be carried out via Mølmer-Sørensen gates with fidelities above 99.9% in some experimental setups [26], in addition to singlequbit gates to rotate the basis from x to z. And finally, exp(−iH 3 t) consists of the exponential of sum of tensor products of four Pauli matrices, which can be carried out efficiently with trapped ions [27,28]. Namely, each exponential of a tensor product of four Pauli operators can be implemented via two Mølmer-Sørensen gates and a local gate, together with the necessary single qubit gates to rotate the bases of the Pauli operators in the tensor product to those needed. Given that all the 4-body terms in Eq. (19) commute, they may be carried out sequentially without digital error, namely, with one Trotter step for the whole H 3 term.
The scaling in our protocol is efficient, given that the number of necessary elementary gates in trapped ions, i.e., single and two-qubit gates, is polynomial in the number of interacting fermions, N , in the Agassi model. On the other hand, with a classical computer the scaling would be inefficient given that the Hilbert space dimension would grow exponentially in N . Of course, under highly symmetric configurations one may obtain a solution, in some cases, in terms of polynomial resources. However, in general terms, with a generalized Agassi model with, e.g., inhomogeneous couplings, this will not be possible and a quantum simulator such as the one proposed here will provide an exponential gain in resources with respect to a classical computer.

Numerical simulations
Note that for all the calculations presented in this section a certain initial state is considered, in our case, | ↓ 1 ⊗ ↓ 2 ⊗ ↑ 3 ⊗ ↑ 4 which corresponds to the state with the minimum value of the angular momentum projection, J 0 = −1 (see definition of J 0 (10)).
We plot in Fig. 2 our numerical results for the digital decomposition. In Figs. 2a we show the fidelity | φ(t)|φ(t) T | 2 as a function of (g + V )t with n T = 10, where |φ(t) and |φ(t) T denote the exact numerical state and the one obtained via the Trotterized digital dynamics, respectively. Panel c) is just a zoom of a) for small t. In Figs. 2b and 2d we depict the fidelity | φ(t)|φ(t) T | 2 as a function of n T for t = 5 and t = 1, respectively, where t denotes the total simulated time interval. The red dots in panels b) and d) correspond to the red dots in panels a) and c), respectively. This figure makes clear that the Trotter dynamics match very efficiently the exact calculation in a large time interval even for small number of Trotter steps, n T . Calculations with larger systems [29] points into the same conclusions raised in this work, however digital quantum simulations for j >> 1 are beyond current technology.
Also, we can observe in Fig. 2 that the digital error remains negligible for a sizeable time evolution with a nontrivial dynamics and a sufficiently large n T . With respect to the total gate error in a plausible implementation with trapped ions, one can estimate its magnitude via adding the single and two-qubit gate errors times the corresponding number of gates. In our specific 4-qubit proposal, there are 52 single-qubit gates and 50 two-qubits gates. If one assumes experimentally achieved values of 0.0001 for the single-qubit gate error [25] and 0.001 for the two-qubit one [26], an estimate for the total gate error E G will be, assuming n T = 5, E G ≃ 5 × (52 × 0.0001 + 50 × 0.001) ≃ 0.28. Thus, with a conservative gate counting, we estimate that the achieved fidelity may be above 70%. Moreover, the number of gates is such that one may perform the experiment well before the decoherence time, in less than ten milliseconds. Therefore, our proposal may be carried out in trapped-ion setups with current technology, for a proof-of-principle model with j = 1, i.e., N = 4. We point out that in this heuristic analysis we have assumed a wellcontrolled experiment with uncorrelated errors. In addition, the protocol is efficiently scalable to many fermionic modes, namely, N ≫ 1, once the single and two-qubit gate fidelities, as well as coherence times, are improved. This is due to the fact that the number of terms in our digital decomposition scales polynomially in N , as opposed to classical supercomputers, for which the scaling would be exponential in the general case. Moreover, the nuclear physics models we consider will always have a polynomial number of terms when expressed as sum of many-body tensor products of Pauli matrices. This is due to the fact that the H 3 Hamiltonian will contain at most products of four c fermionic operators, and this implies that the number of σ + and σ − operators will be at most of four per term, independently of the number N of modes. Regarding the connection to usual observables in nuclear physics, in a quantum simulation experiment such as this, one may compute the quantum state via quantum tomography [30], for systems up to 8 qubits, and the Hamiltonian spectra via quantum phase estimation algorithm [31], which is polynomial (i.e., efficient) in the size of the system. We point out that quantum tomography would only be useful for a quantum experiment with few qubits, as the one explicitly described here. For scaling up the experiment to many qubits, we propose to employ instead two-point correlation functions as shown in Fig. 4, which can be measured directly in trapped ions via resonance fluorescence. In nuclear physics, sometimes observables evaluated in different times are desirable. In this sense, we point out that it is possible to carry out this kind of measurement in a digital quantum simulator, as was proposed, e.g., in Ref. [32,33].
In Fig. 3 we plot the survival probability | φ(t)|φ(0) | 2 as a function of (g + V )t to show that the dynamics of the system is not trivial and significantly changes in the time interval considered in Fig. 2. Finally, in Fig. 4 we depict the correlation function σ z (1, 2) ≡ σ z 1 σ z 2 − σ z 1 σ z 2 , showing that the time dynamics alone can be used as a probe to explore the different quantum phases of the system via this correlation function. As mentioned above, the critical point in the Agassi model for j = 1 is given by g + V = 1. Fig. 4 shows three calculations for this correlation function: a) g + V < 1 (symmet-ric phase), b) g + V = 1 (phase transition point), and c) g + V > 1 (non-symmetric phase). One can clearly see that at one side of the phase transition the correlation amplitude maximum is smaller than one (symmetric phase, Fig. 4a with g + V < 1), it is already one at the transition point (Fig. 4b, with g + V = 1) continue being 1 at the other side (broken phase, Fig. 4c, with g + V > 1 ), and extra oscillations appear which amplitudes depend on the control parameter value (also visible in Fig. 4c). A quite similar behaviour is obtained for other initial states and correlation functions, such as σ z (1, 3) or σ z (3, 4), obtaining also a maximum amplitude for g + V = 1. This is more clearly shown in Fig. 5, in which the amplitude of the oscillation is plotted as a function of the control parameter g = V , where the critical point is at g = V = 0.5. The figure shows that this amplitude reaches the value 1 at the critical point and keeps this value in the broken-symmetry phase (g + V ≥ 1). In the upper part of the figure the phase diagram of the model is sketched, separating the symmetric phase (SP) and the broken-symmetry phase (BSP). This is an evidence that one does not need to compute the ground state to distinguish the different quantum phases in the system. The most direct time dynamics for a typical initial state allows one to obtain signatures of these quantum phases in the amplitude of the time dynamics of the correlation function. This type of procedure resembles a dynamical quantum phase transition (see Ref. [34]), which is a type of quantum phase transition in the time domain. In Ref. [35], the author studied the Rabi model through a quench, noticing that its evolution provides hints on the phase of the ground state of the system. The analysis of the time evolution of the correlation function σ z (1, 2), among other possible functions, provides information on the whole Hilbert space of the system, including its ground state. We plan to extend the present formalism to larger systems and to other observables, such as the Loschmidt echo, to explore how robust are the obtained results. Even very small systems, such as the one we are considering with j = 1, can present some precursors of quantum phase transitions as it was explained in Ref. [36]. Although the quantum phase transition is strictly defined in the thermodynamic limit, already for really small sizes its effect is noticeable in several observables that can act as order parameters.

Conclusions
We have proposed and analyzed the quantum simulation of the Agassi model. Our numerical simulations and analytical estimations show that this protocol is feasible with current technology, for instance, using trapped ions. The proposal has been exemplified with four sites to be implemented with four trapped ions, while it is scalable to many sites with polynomial resources. We also give evidence that the time dynamics of a quantum correlation function for typical initial states can serve as a probe to explore the different quantum phases of the model, with no need of computing specifically the ground state. Indeed, the different phases of the system can be matched to the time dynamics of the amplitudes of the correlation function. With recent advances in trapped-ion quantum platforms approaching a few tens of ions in a quantum processor [37, 38], we are already going through the crossover for outperforming the fastest classical supercomputers for useful scientific problems. Our approach is a step in this direction, for the efficient quantum simulation of the Agassi model and related nuclear physics systems with digital quantum platforms. An appeal of trapped-ion quantum platforms is the all-to-all connectivity that enables one to implement the Nbody tensor products of Pauli matrices with just two Mølmer-Sørensen gates. However, to be able to carry out a full-fledged quantum simulation of the Agassi model with trapped ions, two-qubit gate fidelities will still need to improve. Even though the scaling of our quantum algorithm is polynomial in the number of quantum particles, beyond a few hundred spins one will need to employ a full fledged error corrected quantum computer, with the consequent overhead in resources that will be needed.