Digital quantum simulation of an extended Agassi model: Using machine learning to disentangle its phase-diagram

PHYSICAL REVIEW C 106, 064322 (2022) Digital quantum simulation of an extended Agassi model: Using machine learning to disentangle its phase-diagram Álvaro Sáiz ,1José-Enrique García-Ramos ,2,3José Miguel Arias ,2,4Lucas Lamata ,2,4and Pedro Pérez-Fernández 1,2 1Departamento de Física Aplicada III, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, 41092 Sevilla, Spain 2Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Fuentenueva s/n, 18071 Granada, Spain 3Departamento de Ciencias Integradas y Centro de Estudios Avanzados en Física, Matemática y Computación, Universidad de Huelva, 21071 Huelva, Spain 4Departamento de Física Atómica, Molecular y Nuclear, Facultad de Física, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain (Received 14 June 2022; accepted 29 November 2022; published 21 December 2022) A digital quantum simulation for the extended Agassi model is proposed using a quantum platform with eight trapped ions. The extended Agassi model is an analytically solvable model including both short range pairing and long range monopole-monopole interactions with applications in nuclear physics and in other manybody systems. In addition, it owns a rich phase diagram with different phases and the corresponding phase transition surfaces. The aim of this work is twofold: on one hand, to propose a quantum simulation of the model at the present limits of the trapped ions facilities and, on the other hand, to show how to use a machine learning algorithm on top of the quantum simulation to accurately determine the phase of the system. Concerning the quantum simulation, this proposal is scalable with polynomial resources to larger Agassi systems. Digital quantum simulations of nuclear physics models assisted by machine learning may enable one to outperform the fastest classical computers in determining fundamental aspects of nuclear matter. DOI: 10.1103/PhysRevC.106.064322 I. INTRODUCTION Realistic physical systems, in general, and nuclear physics problems, in particular, are usually too complex to be solved exactly. Thus, many-body approximations have to be used to get reasonably approximate results. Frequently, the complexity of the problem is such that even the use of standard many-body approximations faces insuperable difficulties. To tackle this issue, algebraic models, either for fermions or bosons, having only few degrees of freedom have been developed. These can be solved analytically in some particular scenarios and numerically in general cases, even for very large systems. These facts make them very useful to check approximations and to provide solid references to test more complex and fundamental calculations. Among these types of algebraic models, with relevance in nuclear physics, one can mention the Elliot SU(3) model [1], which played a major role in the explanation of rotational structure emergence from single particle degrees of freedom in nuclei, the interacting boson model (IBM) for the study of low-lying collective states in medium and heavy nuclei [2], or the Lipkin-Meshkov- Glick (LMG) model [3] introduced first to check approximate many-body techniques in nuclear models with many nucleons interacting via a monopole-monopole force. Those are examples of simple models that have been able to permeate other fields of physics, e.g., condensed matter, quantum optics, or molecular physics, due to its simplicity and ability to capture the key aspects of a large variety of physical systems. The Agassi model was proposed as a simple and exactly solvable version of the far more complex pairing-plus- quadrupole model of nuclear physics [4]. It includes a combination of a two-level pairing model and a monopolemonopole interaction such as the one of the LMG model. The extended Agassi model was introduced in Ref. [5] and its main difference with respect to the original model is the addition of a more general pairing interaction. The model presents a rich quantum phase diagram and, specifically, gives rise to quantum phase transitions of different nature. This makes the study of the extended Agassi model of much interest in the field of quantum phase transitions (QPT’s), in particular, to test new techniques able to identify the phase of the system through the study of appropriate observables. Quantum simulations provide a fast-developing and powerful tool to provide realizations of countless quantum systems of great interest, such as spin models, quantum chemistry, as well as QPT’s [6]. In the near future, quantum simulators should be able to outperform classical computers and solve previously intractable problems. Consequently, many experimental setups are being proposed to validate the feasibility to deal with different physical models [7]. The usefulness of the Agassi model and its relevance in a variety of fields, including that of QPT’s, motivates its quantum simulation. Furthermore, tools from quantum information have been recently used, as the quantum discord, to explore the phases in this model to gain insight about its structure [8]. The quantum simulation of the original Agassi model has been already proposed for a system of trapped ions with 2469-9985/2022/106(6)/064322(12) 064322-1 ©2022 American Physical Society ÁLVARO SÁIZ et al. PHYSICAL REVIEW C 106, 064322 (2022) four sites in [9], and the quantum phase diagram of the extended Agassi model has been extensively studied through a mean-field formalism [10]. In the present work, the quantum simulation of the extended Agassi model with eight sites in a system of trapped ions is worked out, as a natural followup to these previous studies. Note that the analysis is performed in a classical computer and the exact results are compared with the ones from a Trotter expansion which is used to simulate the quantum computer. In the case of the standard Agassi model with four sites, its quantum simulation only presents two phases: a symmetric and a broken symmetry phase. However, the extended Agassi model with eight sites presents a richer phase diagram with five different quantum phases. This makes the model specially issued to be studied for two reasons: (i) It is at the limit of what the current resources in facilities with trapped ions could do and (ii) its rich phase diagram allows evaluating and validating the use of machine learning methods to identify the phase of the system. The phase/shape of a system is a property of its ground state. However, in Ref. [9] it was proved that it can also be determined through the time evolution study of selected operator matrix elements acting on a non eigenstate of the system. However, it is not always obvious how to determine the phase of the system through this procedure. We can find an alternative solution to this challenge using another technique that has seen a fast growth in the last decade, namely, machine learning. Thanks to its great versatility, it has been tested with large success in many and diverse fields of science, such as medicine [11], biology [12] and astrophysics [13]. Machine learning algorithms are specially suitable to identify features that are not immediately apparent and they work as ideal classifiers [14]. In this paper, two different machine learning approaches are presented to aid in the determination of the system phase: One simple approach using a multilayer perceptron (MLP) and another more complex and precise approach, using a convolutional neural network (CNN). Both are able to accurately predict the quantum phases of the system just taking as an input the quantum simulation of the time evolution of an appropriate operator. The paper is organized as follows. In Sec. II, the extended Agassi model Hamiltonian is presented and then it is written in terms of the Pauli matrices using the Jordan-Wigner spin mapping approach. In Sec. III, the feasibility of a digital quantum simulation of the model is studied through its Trotterized dynamics that can be implemented efficiently with trappedion systems. In Sec. IV, a brief study of the QPT’s of the model is presented. In Sec. V, the deep learning techniques, used for the determination of the quantum phase are presented and their results are reported. Finally, in Sec. VI, the conclusions of the paper are summarized. II. SPIN MAPPING OF THE EXTENDED AGASSI MODEL The Agassi model is a two-level system in which the upper level is labeled as σ=1 and the lower level as σ=−1. Each level has degeneracy, being an even number [4]. One can then introduce jsuch that =2j. Furthermore, in order to label the states within a given level, it is introduced a magnetic quantum number m=±1,±2,...,±/2. In Ref. [5], an extended Agassi model was considered, which includes a more general pairing interaction with the additional term −2hA† 0A0[see Eq. (1)]. The Hamiltonian for the extended Agassi model is written as H=εJ0−g σ,σ=−1,1 A† σAσ −V 2[(J+)2+(J−)2]−2hA† 0A0.(1) The operators of the Hamiltonian are all defined in terms of fermion creation and annihilation operators, J+= j  m=− j c† 1,mc−1,m=(J−)†,(2) J0=1 2 j  m=− j (c† 1,mc1,m−c† −1,mc−1,m),(3) A† 1= j  m=1 c† 1,mc† 1,−m=(A1)†,(4) A† −1= j  m=1 c† −1,mc† −1,−m=(A−1)†,(5) A† 0= j  m=1 (c† −1,mc† 1,−m−c† −1,−mc† 1,m)=(A0)†,(6) where c† ±1,mand c±1,mare the fermion creation and annihilation operators, respectively. Using the Jordan-Wigner spin mapping approach [15,16], the Hamiltonian of the extended Agassi model can be mapped into Pauli matrices, which can later be experimentally implemented into a digital quantum simulator, as it will be discussed in the following section. In the following, a system with eight sites (j=2) will be considered. Hence, to simplify the notation, the fermions are relabeled as c1,2→c1,c−1,2→c5, c1,1→c2,c−1,1→c6, c1,−1→c3,c−1,−1→c7, c1,−2→c4,c−1,−2→c8, (7) and the corresponding relabeling for the creation operators holds. The spin image of the above fermions is c† i=I1⊗... ⊗Ii−1⊗σ+ i⊗σz i+1⊗... ⊗σz N, ci=I1⊗... ⊗Ii−1⊗σ− i⊗σz i+1⊗... ⊗σz N, (8) where σ±are written in terms of the Pauli matrices in the standard way σ+=σx+iσy 2,σ −=σx−iσy 2,(9) σx,y,zare, respectively, σx=01 10 ,σy=0−i i0,σz=10 0−1,(10) 064322-2 DIGITAL QUANTUM SIMULATION OF AN EXTENDED … PHYSICAL REVIEW C 106, 064322 (2022) ⊗stands for the Kronecker product and Ifor the identity operator. The Hamiltonian can be written in terms of strings of tensor products of Pauli matrices. Indeed, there are terms which only depend on the Pauli matrix σz, and there are others that need to be written in terms of the Pauli matrices σxand σy, so that they can be mapped into the corresponding quantum gates. Once σ+and σ−are decomposed in terms of σxand σy, each of the tensor product strings will give rise to 8 different terms. The total number of terms scales polynomially with the number of sites. III. FEASIBILITY OF THE IMPLEMENTATION Given a certain time-independent Hamiltonian, the evolution operator is written as U(t)=e−itH ,(11) considering ¯h=1. For the case of a Hamiltonian written as the sum of k different noncommuting parts, the evolution operator can be approached using the so-called Trotter expansion as U(t,n)T= k e−itHk/nn ≈U(t).(12) The study of the many Hamiltonian terms arising from the Jordan-Wigner mapping shows that many of them commute with each other. Hence, this can be exploited by adding these terms together in a total of eight noncommuting groups [k=8inEq.(12)], saving computational resources. Each of these groups can then be implemented with trappedion systems using a combination of single and multiqubit (Mølmer-Sørensen) gates and simulation of the terms within each single group can be carried out sequentially without digital error [17–19]. The accuracy of this approximation can be measured by its fidelity, defined as F(t,nT)=|φ|U(t,nT)† TU(t)|φ|2.(13) This observable depends on time, t, and on the number of Trotter steps nT. The considered initial state to which we apply the evolution will be along the rest of the paper |↓ 1↓2↓3↓4↑5↑6↑7↑8, where |↓ iand |↑ istand for the eigenvectors of σz iwith eigenvalues −1 and 1, respectively. Hence, the trial state is not a Hamiltonian eigenstate, except for V=g=h=0. Figure 1shows the fidelity as a function of both tand nT. In a previously studied case with four sites (j=1), it was seen that already nT=5 Trotter steps provided enough accuracy [9]. Here, the fidelity drops faster as time increases and for nT=5 the approximation only remains accurate at small times (t<1) [panel (a)], but this drop is greatly reduced for nT=15 [panel (c)]. This is expected, since eight sites are considered and, therefore, there exist many more terms that need to be considered for the Trotterized Hamiltonian, which in turn means that a significantly higher amount of Trotter steps are required to attain the same accuracy in a given time range. In panels (b) and (d) the value of the fidelity is depicted FIG. 1. Fidelity as a function of tfor (a) nT=5and(c)nT=15 andasafunctionofnTfor (b) t=0.5and(d)t=5. The initial state in all cases is |↓ 1↓2↓3↓4↑5↑6↑7↑8and the Hamiltonian parameters are ε=1, g=0.25, V=0.25, and h=0.25. The highlighted points correspond to nT=5,15 and t=0.5,5, respectively. as a function of the Trotter steps, nT, for two different values of time Figure 2shows as an example the survival probability, |φ(t)|φ(0)|2, computed for both the exact and the Trotterized evolution operator with different values of V=0.25 (top row) and V=1.5 (bottom row) and nT=5,15 [panels (a) and (b)] and nT=5,30 [panels (c) and (d)], keeping constant g=0.5 and h=1.5. For a small value of V=0.25, a good fidelity can be achieved with only nT=5 Trotter steps. However, this is not true in general. There are situations where an accurate study require a larger number of Trotter steps as it is seen in Fig. 2, panels (c) and (d). The implementation of the simulation is efficient, that is, it scales with polynomial resources. However, considering the need of two Mølmer-Sørensen gates for each tensorial product and two single-qubit gates for each σxand σyto be changed into a σzand back (assuming σzas the basis), around two thousand quantum gates would be needed to implement this model for j=2. It is important to note that the aggregated error from all the gates will then increase with the number of Trotter steps needed for a precise approximation. Nevertheless, the estimations seem to agree that this proposal should be achievable in the near future as a good progress is being made in scaling trapped ion systems to the hundreds or thousands of qubits [20,21]. Some approaches with promising results are the implementation of high fidelity fast gates [22,23] and two-dimensional ion traps, since the fidelities in each dimension are independent from one another [24]. The use of a higher-order Trotter expansion would also help 064322-3 ÁLVARO SÁIZ et al. PHYSICAL REVIEW C 106, 064322 (2022) FIG. 2. Survival probability |φ(t)|φ(0)|2as a function of tusing the exact and the Trotterized evolution operator. In all cases, the initial state is |↓ 1↓2↓3↓4↑5↑6↑7↑8and ε=1, g=0.5, and h=1.5. Top row:V=0.25 and (a) nT=5and(b)nT=15. Bottom row: V=1.5and(c)nT=5and(d)nT=30. reducing the required number of Trotter steps for an accurate approximation [15]. IV. EXPLORING PHASE TRANSITIONS The original Agassi model presents a rich quantum phase diagram with a symmetric, a parity broken phase and a superconducting one [25]. By adding the additional pairing interaction term of the extended Agassi model, the quantum phase diagram becomes even richer. In order to explore the phase diagram and the appearance of QPT’s it is customary to rescale the parameters of the Hamiltonian (1) by introducing the new parameters χ,, and , V=εχ 2j−1,g=ε 2j−1,h=ε 2j−1.(14) The phase diagram of the extended Agassi model presents five different quantum phases, as shown in Fig. 3.They appear in the following areas of the phase diagram (see Ref. [10] for a detailed study of the quantum phases of the system): (i) Symmetric or spherical solution: χ<1,<1,<1. (ii) Hartree-Fock (HF) deformed solution: χ>1,<χ,< 1+χ2 2χ. (iii) Bardeen-Cooper-Schrieffer (BCS) deformed solution: χ<, >1,< 1+2 2. FIG. 3. Graphical representation of the phase diagram of the extended Agassi Hamiltonian. Red vertical planes represent second order QPT surfaces. The green surface (=1forχ<1and <1, =1+χ2 2χfor χ>,and=1+2 2for χ<) and the blue vertical one (χ=and < 1+2 2) correspond to first order critical surfaces. Red sphere, blue oval, black oval, black thick oval, and crossed green ovals correspond to the symmetric solution, the HF deformed solution, the BCS deformed solution, the closed valley solution, and the HF-BCS deformed solution, respectively. Figure adapted from [5]. (iv) Combined HF-BCS deformed solution: >1forχ<1 and <1; >1+χ2 2χfor χ> >1; > 1+2 2for >χ>1. (v) Closed valley solution: χ=,  < 1+2 2. It is important to note that this phase diagram has been obtained using a mean-field approach, which is only valid for the large Nlimit, where Nis the number of sites, and it allows to calculate the value of the order parameters as a function of the control parameters. Here, N=8 and therefore, finite-N effects become relevant, causing the QPT’s to be diffuse and harder to be distinguished. In previous studies, time evolution of a correlation function has been successfully used as a tool to explore the QPT’s of the standard Agassi model with j=1[9]. In the rest of this work, we will use precisely this observable to study the phase diagram of the extended Agassi model with j=2. A correlation function between two given sites can be defined as Cα,β (i,j)=σα i⊗σβ j−σα iσβ j,(15) where α, β =x,y,or z. In our case, we will focus on a correlation function for sites 1 and 2 and α=β=z, Cz(1,2) ≡σz 1⊗σz 2−σz 1σz 2.(16) Figure 4(a) shows the systematics of the correlation function Cz(1,2) as we move through a line of fixed χ=1.5 and =0 from the HF phase (<χ) to the BCS phase (>χ) as a function of time. The correlation function oscillates with time and, despite the fact that the transition is not apparent, if we move away from the critical surface, e.g., at slices A 064322-4 DIGITAL QUANTUM SIMULATION OF AN EXTENDED … PHYSICAL REVIEW C 106, 064322 (2022) FIG. 4. (a) Time evolution of the correlation function Cz(1,2) from =0.5to=2.5, with χ=1.5, =0, and an initial state |↓ 1↓2↓3↓4↑5↑6↑7↑8. The phase transition point =χ=1.5 is marked by a dashed red line and the points χ=1.5, =0.7 (slice B) and χ=1.5, =2.3 (slice A) with blue dashed lines. (c) Maximum value of the amplitude of the oscillations of Cz(1,2) as a function of  value. (b) Time evolution of Cz(1,2) for =2.3 (slice A) and (d) for =0.7 (slice B). and B, we can see a clear difference in the amplitude of these oscillations. Hence, it increases as one moves towards the BCS phase, as shown in panel (c) of Fig. 4. In panels (b) and (d), the time evolution of the correlation function for the points A and B are depicted. These two panels clearly show how the amplitude and the period of the oscillation greatly change depending on the value of the control parameter. In the Appendix, other examples are gathered showing the variety of behaviors that one can get depending on the position in the phase diagram. Other correlation functions between different sites, applied to the same and/or different initial states, yield similar results (not shown in this work), except when they vanish, which proves that correlation functions are good observables to obtain information about the phase of the system without the need of computing the ground state. Even in the case where differences between quantum phases can be appreciated in the time evolution, it is hard, and sometimes impossible, to visually identify the quantum phase of the system, specially when looking at a single point in the phase space. Hence, two machine learning models will be trained for quantum phase classification using the time dynamics of the correlation function Cz(1,2). The benefit of using machine learning methods for classification is that, given the data to classify, the method can determine itself the optimal procedure and, moreover, it will be able to accurately label not previously considered cases. V. MACHINE LEARNING METHODS FOR CLASSIFICATION A. Multilayer perceptron A multilayer perceptron (MLP) is a supervised machine learning algorithm that learns a nonlinear function which, given a set of features as input X=x1,x2,...,xm, is able to approximately predict a target output y. The algorithm learns this function via gradient descent, modifying its weights by minimizing a loss function based on how accurate the predictions yare in relation to the real target values y, which are known in the training process, as with any supervised method [26]. Once optimized, it can reliably predict the target y≈y for inputs for which the real target is not already known. The implementation of this model is very simple thanks to the existing standardized machine learning software packages, e.g., SCIKIT-LEARN [27], but first it is needed to define a set of training and test data for the model to learn with. For this purpose, we compute the exact as well as the Trotterized (using a number of Trotter steps nT=6) time evolution of the correlation operator Cz(1,2) for a large set of points throughout the phase diagram. We create a mesh of points with the parameters χ,, and in the interval [0,2] with 21 steps each for a total of 9261 points. For each of them, we compute the time dynamics for 100 time steps. This information together with the knowledge of the phase diagram (see Sec. IV) feeds the model as part of the training process. 064322-5 ÁLVARO SÁIZ et al. PHYSICAL REVIEW C 106, 064322 (2022) Once the dataset is ready, it is randomly split into two sets. One set, the training set, will be used to train the model, while the other, the test set, will be used to test the accuracy of the model. It is compulsory for the model to be tested with data not used during the training process to be sure that the algorithm is not simply remembering the already used data. Therefore, it is important to have enough data points for testing and to reach an acceptable metric, while keeping as many points as possible for training. In most cases, to use from 10 to 20 % of data for testing is enough. In our case, we will use 10% of data points for testing. SCIKIT-LEARN provides us with several optimization algorithms to perform training on the data. In this particular case, the default Adam solver gives optimal results. Adam solver is an optimization algorithm that is based on the stochastic gradient descent [28]. It combines the advantages of two other extensions of the stochastic gradient descent known as adaptive gradient algorithm (AdaGrad) and root mean square propagation (RMSProp). For a general overview on the theoretical and practical aspects of the learning algorithms based on stochastic gradient descent, see, e.g., Ref. [29]. The model actually computes the probability for each quantum phase and classifies the data based on which one is the most probable. The classification with the MLP correctly predicts the quantum phase in the 96.7% of the points in the test data when trained and tested with the exact time evolution. When the MLP is trained and tested with the Trotterized time evolution, it has an accuracy of 93.4% (nT=6). These results are very valuable, since it proves that the quantum phase can be determined using the correlation function without the need of computing the ground state, but it also shows that they can be identified despite using the Trotter approximation with a relative small number of steps, for which the agreement with the exact results is limited. Note that when reducing the number of Trotter steps, the aggregated error of the quantum gates is reduced too. In Fig. 5, the results for the MLP method are presented through different trajectories in the parameter space of the model. Generally, the incorrect predictions occur near the critical points, which is expected, since the QPT’s are not clearly defined because of finite-Neffects (N=8 in our case). Note that at the critical points (denoted by a dashed line), two or three quantum phases can coexist, therefore, it is considered that the model correctly classifies the quantum phase of the system when it predicts either of the two phases, except in the case of the closed valley solution, where that phase is considered the only valid answer. This is done to ensure that the model learns to identify this specific phase instead of choosing any of the other two phases with which it coexists. Perhaps the most surprising case is, in fact, the one regarding this solution. In the results portrayed in Fig. 5(d), the model is able to predict the quantum phase of the system as it crosses the closed valley solution at χ==1.5 when trained with the exact solution. As mentioned, the model is not completely accurate in some cases [as can be seen from Fig. 5(c)-Trotter], but the QPT’s are still fairly clear. In next section we will present a more accurate machine learning model. B. Convolutional neural network Convolutional Neural Networks (CNN) are types of deep learning algorithms that are well suited for image recognition. Their efficacy for extracting features and learning to classify images has been proven in multiple instances [30–32]. In our case, instead of feeding two-dimensional (2D) arrays of pixel values, one-dimensional (1D) arrays containing the time dynamics of the correlation function will be used. When working with 2D CNN, the data have actually the shape (x,y,z), being zthe value contained in the pixel, while xand yare its positions. 1D CNN work effectively in the same way, just with data of shape (x,y), being ythe value of the variable of interest at time x. The general structure of a CNN includes more than one convolution layer, so that one can extract higher level features. They also include additional layers, e.g., the pooling and the fully connected layers, each with a different purpose. For the present model the used layers are (i) Convolution layer: The layer responsible of performing the convolution operation. (ii) Activation layer: The layer that applies the activation function together with the filter of the convolution layer. (iii) Pooling layer: The pooling layer performs a dimension reduction of the data, collapsing data by connecting clusters of neurons to a single neuron each. (iv) Dropout layer: This optional layer temporarily deactivates, or drops out, randomly selected training parameters from the previous layer that has trainable parameters. Its goal is to avoid the “overfitting”. (v) Fully connected layer: Also know as dense layers, they connect every neuron of the input to every neuron of the output. (vi) Softmax layer: This layer is a fully connected or dense layer that applies a specific kind of activation function, called a softmax function, which is a normalized exponential function. A fast and efficient way of implementing such a model is through the use of the Keras library [33], a deep learning application programming interface built on top of the machine learning platform Tensorflow [34]. Keras provides the ability to apply all the necessary layers sequentially, giving the model the final structure shown in Table I. Four cycles are applied, consisting on convolution, activation, pooling and dropout, allowing for a hierarchical decomposition of the input and the extraction of higher level features, as previously mentioned. Similarly, five additional fully connected (dense) layers are applied before the final output layer. The model includes a layer named “Flatten”. This layer simply reshapes the output of the previous layer into a single dimension array so that it can be passed as the input to the dense layer. Table Ishows, alongside the layers, the shape of the output and the number of trainable parameters of each layer. Note that only the convolution and dense layers have trainable parameters and the rest of the layers simply apply a chosen transformation to the data which is never changed throughout the training process. The convolution layers are applied with 064322-6 DIGITAL QUANTUM SIMULATION OF AN EXTENDED … PHYSICAL REVIEW C 106, 064322 (2022) FIG. 5. Quantum phase prediction of the system via the MLP method applied to theCz(1,2) correlation function. The different panels show the probability for the system of being in a given phase for specific values of the Hamiltonian parameters. Left panels correspond to the results using the exact time evolution, while right ones using the Trotter approach with nT=6. In all cases, the initial state |↓ 1↓2↓3↓4↑5↑6↑7↑8is used. (a) =0.5and=0.5movingfromχ=0 (symmetric phase) to χ=2 (HF phase); (b) χ=0.5and=0.5 moving from =0 (symmetric phase) to =2 (BCS phase); (c) χ=0.5and=0.5movingfrom=0 (symmetric phase) to =2 (combined HF-BCS phase); (d) χ=1.5and=0.5movingfrom=0 (HF phase) to =2 (BCS); (e) χ=1.5and=0.5movingfrom=0 (HF phase) to =2 (combined HF-BCS phase); (f) χ=0.5and=1.5movingfrom=0 (BCS phase) to =2 (combined HF-BCS phase). The dashed black line in each panel marks the theoretical critical point. padding to match the first dimension of the output with the size of the corresponding input. The second dimension will be equal to the number of applied filters, 32, 64, 128, and 256 for the first, second, third, and fourth convolution layers, respectively. For convolution layers, the trainable parameters define each of the filters. Since the filters have a size 3, each one will have three parameters or weights for each of the input channels, plus one bias weight. This number, multiplied by the number of applied filters, gives the total number of trainable parameters shown on the right column of Table I. In the case of the dense layers, each connection has its own weight and all input neurons are connected to all output ones, so the 064322-7 ÁLVARO SÁIZ et al. PHYSICAL REVIEW C 106, 064322 (2022) TABLE I. Number of parameters and shape of the output of the used convolution layers. Layer Output shape # of parameters Convolution 1D (1) 100 ×32 128 Leaky ReLU (1) 100 ×32 0 Avg. pooling 1D (1) 34 ×32 0 Spatial dropout (1) 34 ×32 0 Convolution 1D (2) 34 ×64 6208 Leaky ReLU (2) 34 ×64 0 Avg. pooling 1D (2) 12 ×64 0 Spatial dropout (2) 12 ×64 0 Convolution 1D (3) 12 ×128 24704 Leaky ReLU (3) 12 ×128 0 Avg. pooling 1D (3) 4 ×128 0 Spatial dropout (3) 4 ×128 0 Convolution 1D (4) 4 ×256 98560 Leaky ReLU (4) 4 ×256 0 Avg. pooling 1D (4) 2 ×256 0 Spatial dropout (4) 2 ×256 0 Flatten 512 0 Dense (1) 512 262656 Leaky ReLU (1) 512 0 Dropout (1) 512 0 Dense (2) 512 262656 Leaky ReLU (2) 512 0 Dropout (2) 512 0 Dense (3) 512 262656 Leaky ReLU (3) 512 0 Dropout (3) 512 0 Dense (4) 512 262656 Leaky ReLU (4) 512 0 Dropout (4) 512 0 Dense (5) 512 262656 Leaky ReLU (5) 512 0 Dropout (5) 512 0 Softmax 5 2562 Total parameters 1445445 total number of parameters will be the size of the input times the size of the output, plus one bias weight for each output neuron. Once the layer structure is defined, the model should be trained in a similar manner to the MLP case. The data are divided in three sets, namely, the test, the training and the validation set. The training one is used to get the trainable parameters, the validation one to check the consistency of the results and eventually to modify the network layer structure and hyperparameters (such as the number of neurons or filters), finally, the test dataset is used to report the prediction accuracy of the process. When working with neural networks, a common problem is the overfitting of the model, specially when the available data for training is limited. The model may “memorize” the training data and achieve good results for that input specifically, but becomes unable to generalize and predict correctly for data outside of the training set. In Fig. 7, training and validation accuracy are shown across the epochs of the training process. In the training case, the accuracy in predicting the phase is computed using the training data, while in the validation accuracy a separate set of data is used, i.e., the validation set. Note that neither of these sets contain the test dataset. One can see that both training and validation accuracies steadily increase at a similar pace, which suggests that the model is not overfitted. In the case of overfitting, the training accuracy would increase faster than the validation accuracy [14] or the latter would stagnate at an earlier epoch, which is not the case. The total number of training epochs considered in this work was 200, as further cycles do not significantly increase the accuracy, assuring that the model has converged. In our case, we get an accuracy in the prediction of the phase, using the test dataset, of a 98.7% (corresponding to the exact evolution). Note that this reported accuracy do not use data neither from the validation dataset nor from the validation one. The increase in accuracy from 96.7% using the MLP to 98.7% using the CNN model is small, but significant. In particular, now the model has a better performance in the most complicated regions, i.e., near the QPT’s. This is clearly seen in Fig. 6, especially for the cases b), c), and d) where the QPT becomes much more clear. The main reason is that not only the model is more precise, but also that inaccuracies are mainly limited to points with values of χand/or equal to or near zero, with only a few cases occurring at the critical points. The reason for these inaccuracies is that the amplitude of the correlation function is very small when both χand are small, and as such, it presents very little variation as one moves across . Hence, this inaccuracy is mostly present between the symmetric and the combined HF-BCS phases. In both cases, the time evolution is quite similar but the model seems to slightly favor the combined HF-BCS phase as the amplitude is generally lower in this case. It is worth to mention once again, that the results are robust to the error introduced from the Trotter approximation even in the case of a relatively small number of Trotter steps. As a matter of fact, the model was trained and tested using the Trotterized time dynamics of the correlation function Cz(1,2) computed for nT=6 and nT=20 with accuracies of 99.2% and 98.9%, respectively, to be compared with 98.7% of the exact evolution. As a conclusion of this section, we can say that the accuracy in determining the phase is larger for the CNN than for the MLP method, which is expected, specially in regions nearby the QPT’s. Moreover, a really appealing fact is that the accuracy of the procedure is even larger when using the Trotter approximation with a small value of nT, which has clear practical advantages. A possible explanation is that the larger oscillations observed in the approximate evolution, compared with the exact one, somehow helps the CNN to better predict the most problematic situations (see Appendix). The use of a relatively small number of Trotter steps reduce the number of needed quantum gates and, therefore, will reduce the experimental error. Because the above reasons, the practical implementation of this algorithm in a real quantum computer seems to be feasible and potentially valid to determine the phase of the system. As a matter of fact, the impact of noise on the learning accuracy, assuming it is small, it is not expected to affect the phase prediction. 064322-8 DIGITAL QUANTUM SIMULATION OF AN EXTENDED … PHYSICAL REVIEW C 106, 064322 (2022) FIG. 6. Same caption as Fig 5but using the CNN. C. Further exploration Other correlation functions, as well as different initial states were studied. There is a lot of combinations of different possible correlation functions and different initial states to evolve. On the one hand, it is obvious that equivalent correlation functions will give the same results, in this case Cz(1,2) =Cz(1,3) =Cz(2,4) =Cz(3,4) =Cz(5,6) = Cz(5,8) =Cz(6,8) =Cz(7,8) and their conjugates for the initial state |↓ 1↓2↓3↓4↑5↑6↑7↑8. On the other hand, the rest of the studied cases give comparable, although less accurate results. The CNN trained and tested on Cz(1,4) for the initial state |↓ 1↓2↓3↓4↑5↑6↑7↑8obtained an accuracy of 94.2%, while using the correlation function Cz(1,2) for a different initial state |↑ 1↑2↑3↓4↑5↓6↓7↓8achieved only an accuracy of 86.5%. Initial states that are eigenstates of one or more of the terms of the Hamiltonian are unsuitable as the information of the time evolution is partially (or totally) lost. Similarly, some correlation functions may vanish, for instance, σxor σy when using initial states in the zbasis. Even when computing the time dynamics of the correlation function Cx(1,2) for an initial state |↓ x 1↓x 2↓x 3↓x 4↑x 5↑x 6↑x 7↑x 8(where |↓ x iand |↑ x i stand for the eigenvectors of σx iwith eigenvalues −1 and 1, respectively) which does not vanish, the CNN was unable of providing valuable results, i.e., it fails in finding the phases of the system. 064322-9