The extended Lipkin model: proposal for implementation in a quantum platform and machine learning analysis of its phase diagram
S. Baid A. Sáiz L. Lamata P. Pérez-Fernández A. M. Romero A. Ríos J. M. Arias J. E. García-Ramos





  ^1Departamento de Física Atómica, Molecular y Nuclear, Facultad de Física, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain

  ^2Departamento de Física Aplicada III, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, E-41092 Sevilla, Spain.
  
  ^3 Departament de Física Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), c. Martí i Franqués, 1, 08028 Barcelona, Spain

  ^4 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), c. Martí i Franqués, 1, 08028 Barcelona, Spain

  ^5Departamento de  Ciencias Integradas y Centro de Estudios Avanzados en Física, Matemática y Computación, Universidad de Huelva, 21071 Huelva, Spain.

  ^6Instituto Carlos I de Física Teórica y Computacional,  Universidad de Granada, Fuentenueva s/n, 18071 Granada, Spain

  ^aNew address: Fujitsu Research of Europe, antonio.marquezromero@fujitsu.com  





  Background
In recent years, the implementation of Nuclear Physics models in quantum computers has emerged as a promising and novel area of research. Simultaneously, the study of quantum shape phase transitions in nuclear models has gained significant attention. Specifically, the phase diagram of the Interacting Boson Approximation (IBA) has been extensively explored, particularly in connection with 
large-particle-number-limit considerations. Interestingly, the Extended Lipkin Model (ELM) serves as a valuable alternative for mimicking the IBA phase diagram and holds the advantage of being more straightforward to implement within a quantum computing platform.



  Purpose
We explore the ELM and provide: i) calculations of the ground state energy using a variational quantum eigensolver; ii) a comprehensive formulation for implementing the dynamics of the ELM within a quantum computing platform, enabling the experimental exploration of the IBA phase diagram for systems with a small number of particles; and iii) a determination of the phase diagram of the model using different Machine Learning (ML) methods. We note that in the ELM, unlike the usual Lipkin model, both first- and second-order quantum shape phase transitions take place depending on the model parameters.



  Method
 Initially, we employ the Adaptive Derivative-Assembled Pseudo-Trotter ansatz Variational Quantum Eigensolver (ADAPT-VQE) to calculate the ground-state energy of the ELM.  Next, we introduce the essential formulation and procedures required to implement this model effectively in a quantum computing environment. Finally, we use ML techniques to identify the different phases and critical points of the ELM. 



  Results
We successfully reproduce the ground-state energy of the ELM across the complete phase space of the model using the ADAPT-VQE algorithm. We provide the necessary framework for implementing the ELM in a quantum computing platform, ensuring that the model can be executed with controlled errors. Finally, we obtain meaningful ML predictions that allow us to determine the phase diagram of the model. 



  Conclusions
Our findings offer compelling evidence that the implementation of a nuclear model like the ELM in a quantum computing environment is not only feasible but can also be achieved with manageable error rates. This realization opens the door to detailed experimental investigations of the phase diagram of the ELM (and indirectly of the IBA) in a quantum computer, further advancing our understanding of quantum shape phase transitions and nuclear structure.







The extended Lipkin model: proposal for implementation in a quantum platform and machine learning analysis of its phase diagram
    S. Baid^1, A. Sáiz^1,2, L. Lamata^1, P. Pérez-Fernández^2, A.M. Romero^3,4,a
, A. Rios^3,4, J.M. Arias^1 and J.E. García-Ramos^5,6
    April 29, 2024
======================================================================================================================================





§ INTRODUCTION


A phase is a specific state of matter characterized by a well-defined structure and physical properties. Phase transitions correspond to abrupt changes in some physical properties of a system as it passes from one phase to the other <cit.>. Phase transitions are associated with sudden shifts in the internal energy of the system that occur at specific values of a control parameter, that is, a physical magnitude that can be changed externally to the system. The point of change is called critical point and, in addition to changes in energy, qualitative changes in other physical quantities, called order parameters, also appear. Order parameters behave differently depending on the phase, so they can potentially be used to distinguish among different phases. Phase transitions are well defined and characterized in macroscopic systems <cit.>.

When one looks into mesoscopic quantum systems, a natural question is whether one can find similar phase transitions to the macroscopic ones or not. One may also wonder whether the same concepts introduced to study classical phase transitions apply at the quantum level. 

The study of Quantum Phase Transitions (QPTs) is a popular and widely discussed topic that encompasses various branches of quantum many-body physics <cit.>. In recent years, there has been a growing interest in investigating low-dimensional lattice models due to advances in quantum computing <cit.>. In addition, concepts from quantum information theory have been used to characterize quantum critical phenomena <cit.>. On a different note, there has been a resurgence in the examination of structural changes in finite-size systems, where early signs of transitions can be observed <cit.>.

In the realm of atomic nuclei, the Interacting Boson Approximation (IBA) <cit.> provides a straightforward, yet detailed, framework to study first- and second-order phase transitions. In the IBA, there are three limiting situations, called dynamical symmetries, in which exact analytical solutions to the complex nuclear many-body problem exist. These are called U(5), SU(3) and O(6) dynamical symmetries. Remarkably, each of these symmetries corresponds to a specific shape in terms of the usual nuclear deformation parameters: β, that measures the deviation from sphericity, and γ, that measures the deviation from axial symmetry. U(5) corresponds to spherical shapes, SU(3) to axial prolate deformed shapes and O(6) to γ-unstable deformed shapes <cit.>. Thus, the phase space of the IBA is usually represented as a triangle such that each dynamical symmetry is located at one vertex <cit.>. We show this so-called Casten triangle in Fig. <ref>. In this figure, the critical line delineating the phase transition from spherical to deformed shapes is depicted by a solid arc connecting the e(5) critical point (representing a second-order phase transition on the U(5)-O(6) line) and the x(5) critical point (representing a first-order phase transition on the U(5)-SU(3) line). The dashed line (anti-spinodal) and dotted line (spinodal) demarcate the regions where shape coexistence occurs within the model.

At the mean-field limit, there are two phases in the model: a spherical phase and a deformed one. The transitions between these phases are of first or second order depending on the path followed in the model parameter space. For instance, the transition from the spherical U(5) dynamical symmetry to the deformed γ-unstable O(6) dynamical symmetry is exactly solvable. This corresponds to the boson pairing model initially solved by Richardson <cit.>, and later applied to the IBA <cit.>. These methods provide access to solutions for a very large number of bosons along this line.  Graphically, this line corresponds to one side of the Casten triangle →. This is the only integrable line in the IBA phase diagram and there, a second order QPT exists. 
The interior region of the Casten triangle, as well as the other two sides, remain beyond the scope of large-scale numerical studies that are limited by the capabilities of standard IBA codes, which can handle systems with approximately 100 bosons <cit.>. In all this area a critical line of first order QPTs appears. 
Other phase diagrams corresponding to algebraic models  have been studied, e.g., the Agassi model <cit.>. The model has been implemented in a quantum platform and the phases of the system was determined resorting to the analysis of the time evolution of a properly chosen operator <cit.>. In a further analysis, the phase diagram has been determined using Machine Learning (ML) methods <cit.>.

Strictly speaking, QPTs only occur in macroscopic systems. One may consider different avenues to connect QPTs to the physics of few-body systems. First, we may want to gain access to numerical simulations that can address low as well as large particle numbers (large-N limit). A second complementary approach would be to encode the model in a quantum platform, where large-N scalings may be tamed down in the long term. These two research avenues are particularly challenging for the IBA, but other models may not suffer from the same limitations. The purpose of this paper is precisely to overcome the limitations of the two analyses suggested above with a model that is similar, but not as challenging, as the IBA. 





We use here the extended Lipkin model (ELM) proposed in Ref. <cit.>. In the large-N limit, the ELM mimics the IBA-shape phase diagram. Moreover, the ELM is simpler than the IBA, and it can be formulated and solved in terms of spin-1/2 Pauli operators. Consequently, the ELM can naturally be implemented in a quantum platform with standard qubit registers. Furthermore, the numerical solution to the ELM is 
accessible classically for a wide range of particle numbers. In this sense, the ELM is a perfect test-bed for QPTs employing both classical and quantum techniques. In this paper, we exploit a variety of ML and quantum computing methods to assess the different phases of the model at finite N and to demonstrate the usefulness of the model for a wide range of applications. We ultimately hope to prove that this model has the potential to become a new standard  for QPT studies in quantum devices. 




The paper is organized as follows. In Sect. <ref>, we provide a short revision of the usual Lipkin model, the IBA and the ELM. We specifically discuss the relation between the IBA and the ELM phase diagrams. In Sect. <ref>, we present an ADAPT-VQE implementation for the ELM, providing benchmarks with ground-state energies for a wide range of values of the control parameters.  
We introduce a proposal for the implementation of the ELM in a quantum platform, together with an evaluation of the possible errors,  in Sect. <ref>. 
ML techniques are applied to investigate the phase diagram of the ELM in Sect. <ref>.  Finally, Sect. <ref> provides the summary and conclusions.  



§ THE MODELS





 §.§ The Lipkin model


The Lipkin-Meshkov-Glick model (LMG) <cit.> (Lipkin model in short) describes a set of N particles interacting with a long-range interaction. Each particle is allowed to be in a two-level system, upper and lower, separated by an energy gap. If each particle is considered to have spin 1/2, the complete set can be represented by a collective spin S⃗. Defined as S⃗ = ∑_i=1^Ns⃗_i where s⃗_i is the spin of each individual particle. The total collective spin S=N/2 is the relevant sector for this work and it includes the ground state of the system.

The original Lipkin Hamiltonian <cit.> can be written[The original Lipkin Hamiltonian is really H_L = ε J_z - V/2 (J_+^2 + J_-^2) - W/2 (J_+J_- +J_-J_+).
Taking W=V, up to a global scale factor, it can be written as,
H_L = (1-λ) J_z - 4 λ/N J_x^2 that is, apart from a constant factor, Eq. (<ref>). The factor 1/N is included to scale properly the relative strengths of the one- and two-body interactions.] as

    H_L= (1-λ) (S+S_z) - λ/N (S_+ + S_-)^2 = (1-λ) (S+S_z) - 4 λ/N S_x^2  ,


where S_ν = ∑_i=1^N s_i,ν for ν = x,y,z and S_± = S_x ± i S_y. Here, λ is a control parameter, and the Hamiltonian can be diagonalized in the collective basis |S M_S⟩, where S=N/2 is fixed and M_S runs from -N/2 to +N/2 in unit steps and the dimension of the Hilbert space is =N+1. 
In this way, the complete energy spectrum can be obtained as a function of λ. 

Although the Lipkin model has a very simple structure, it displays interesting symmetry properties and phase transitions. 
Its phase diagram has been extensively studied <cit.>. It is known that, in the form given by Eq. (<ref>), the LMG presents two phases in the large-N limit. These symmetric and non-symmetric phases are separated by a second-order shape phase transition at the critical point, λ_c=1/5.
We note that the Lipkin model has become a standard benchmark and test-bed for implementation of quantum algorithms focused on nuclear structure <cit.>. 



 §.§ The interacting boson approximation





The interacting boson approximation (IBA) <cit.> is a successful and widely employed nuclear structure model valid to describe low-lying levels of medium and heavy nuclei. The model is formulated in a second quantization formalism that includes a scalar boson (s boson, L=0) and five quadrupole bosons (d_μ boson, L=2). The IBA has u(6) as dynamical algebra and three specific situations, called dynamical symmetries, can be distinguished for which analytical solutions to the complex nuclear many-body problem are obtained. These dynamical symmetries are named: U(5), SU(3) and O(6). A simplified form of the IBA Hamiltonian is the so-called Consistent-Q Hamiltonian (CQH) <cit.>,

    H_CQH = ϵ  n_d - κ  Q^(χ)· Q^(χ) ,

where n_d is the d-boson number operator, 
Q^(χ) = s^†d̃ + d^† s + χ (d^†d̃)^(2) is the quadrupole operator and ϵ, κ and χ are the model parameters. The phase space of the IBA is usually represented  in the Casten triangle (Fig. <ref>). 

Although the IBA is formulated within the context of the second quantization formalism, it is possible to develop a geometric representation of the model by utilizing the concept of coherent states and the intrinsic state formalism <cit.>. This alternative approach unveils geometric shapes reminiscent of the Bohr model <cit.>, which are associated with the three dynamical symmetries of the IBA. Consequently, within this framework, the U(5) limit corresponds to spherical shapes, while the SU(3) 
limit leads to prolate 
axially symmetric shapes. Additionally, the O(6) limit gives rise to deformed γ-unstable shapes. 


Since the model presents several phases, studies of these and of the corresponding quantum phase transitions are of great interest.
Sometimes, to study phase transitions, the CQH in Eq. (<ref>) is rewritten, taking out a global scale factor, as

    H_CQH = (1-ρ)   n_d - ρ/N  Q^(χ)· Q^(χ),

where N is the total boson number of the system. The specific values of ρ and χ that produce the three IBA dynamical symmetries are given in Fig. <ref>. Excellent reviews presenting the IBA in relation with its phase diagram and the corresponding shape phase transitions are Refs. <cit.>.

However, it is important to recognize that QPTs manifest themselves genuinely in macroscopic systems. Consequently, to study them without using mean-field methods, one needs to consider a significantly large number of bosons, which is not always computationally feasible. 
A possible solution is to look for a simpler model with a similar phase diagram, as the one presented in the next subsection.





 §.§ The extended Lipkin model


The Extended Lipkin Model (ELM) <cit.> is a simple two-level model designed to mimic the IBA phase diagram. Although the Lipkin model was originally formulated in terms of fermions or spin operators, as discussed above, it can be alternatively formulated in terms of scalar bosons using the Schwinger representation <cit.> (see Fig. <ref>). Using this boson formulation, the ELM was introduced as a two-level  model constructed on two types of L=0 bosons, s and t, with Hamiltonian

    H= (1-λ)   n_t-λ/N  Q^(α)· Q^(α),

where the operators n_t, which  corresponds to the number of t-bosons, and Q^α are defined as

    n_t= t^† t,  Q^(α)=(s^† t + t^† s) + α  (t^† t) .

λ and α are two independent control parameters that play the role of ρ and χ in the IBA Hamiltonian, Eq. (<ref>). The total number of bosons N=n_s + n_t is a conserved quantity.

We have deliberately written the ELM two-level boson Hamiltonian in the form of Eq. (<ref>) so that it resembles the QCH of Eq. (<ref>). This way, it becomes clear that s and t play the role of the s and d bosons of the IBA respectively. Clearly, the difference resides in the quadrupole character of the d boson that leads to the u(6) dynamical algebra of the IBA, while the ELM is described by a u(2) dynamical one. Despite this important difference, the ELM captures the main features of the phase diagram of the IBA <cit.>. 

Using mean-field techniques similar to those used in the IBA <cit.>, one can obtain the energy per particle surface for the ELM of Eq. (<ref>) as a function of the control parameters, λ and α, and a variational parameter, β,

    E(λ,α;β)/N = β^2/(1+β^2)^2 [ β^2 (1-(1+α^2) λ) - β 4 αλ + (1-5 λ) ].

This allows us to obtain the classical limit of the Hamiltonian and to produce a phase diagram of the ELM model.



The connection between the Hamiltonian in Eq. (<ref>) and its form in terms of spin operators can be obtained by making the inverse Schwinger transformation, defining the collective spin S⃗ as:

    S_+= t^† s,    S_-= s^† t,    S_z=1/2(t^† t- s^† s).

Using this fermionic formulation, given a number of spin-1/2 particles, 
the sector containing the ground state is given by a collective spin S=N/2 and the model Hamiltonian is written as,

    H_EL= H_L - λ/N[α^2 (S + S_z)^2 - 2 α(S_x (S+S_z) + (S+S_z) S_x )  ],

in such a way that for α=0 the usual Lipkin Hamiltonian, H_L in Eq. (<ref>), is recovered. This Hamiltonian can also be easily diagonalized in a basis given by |S M_s⟩ with S=N/2.











 §.§ The phase diagram of the IBA and the ELM


We employ Eq. (<ref>) to determine the phase diagram of the ELM Hamiltonian in 
Eq. (<ref>) (or Eq. (<ref>) in the fermionic formulation). 

The corresponding ELM phase diagram mimics the IBA phase diagram depicted in Fig. <ref> for the IBA Hamiltonian in Eq. (<ref>). As a matter of fact, the connection between the parameters of both models is rather simple,

    ρ→λ,     χ→ -√(7/2)α .


In addition to the critical line displayed by the solid curve,

    λ_c=1/5+α^2 ,

we show the spinodal (dotted) λ_s and antispinodal (dashed) λ_as lines, given by the expressions

    λ_s =α ^2-√(α ^4+10 α ^2+16)+6/α ^2+10 ,     λ_as=1/5 .


We provide a schematic representation of the phase diagram in Fig. <ref>.
The spinodal and antispinodal lines in this diagram
define the region where the energy surface presents two minima: a spherical one and a deformed one, indicating a coexistence of both phases. In the symmetric phase, there is a single energy minimum at the origin (β=0). When the spinodal line is reached, a second, non-symmetric minimum develops for  β≠ 0. At the critical line, the symmetric and the non-symmetric minima are degenerate. Then, at the anti-spinodal line the symmetric minimum at β=0 disappears to become an inflection point. From there on, the β=0  point becomes a maximum and there are only non-symmetric deformed minima. 


The energy functional given by Eq. (<ref>) is a case of the cusp catastrophe whose germ is x^4 <cit.>. 

For α=0, which corresponds to the usual Lipkin model, there is an isolated point of second-order phase transition, with a β^4 behavior at the origin. In that case, the spinodal, anti-spinodal and critical points coincide at the critical value, λ_c=1/5. For α≠ 0, the phase transition changes its character to first order with the critical point given by Eq. (<ref>). Here, both minima (one spherical and one deformed) are degenerate and there is a coexistence region around the critical point, λ_c, between λ_as and λ_s. Moreover, for λ>λ_as, two non-symmetric minima coexist. 





 §.§ ELM numerical calculations


Having established the similarity of the IBA and the ELM models at the mean field level, we now proceed to present numerical ELM calculations for a finite particle number, N. We aim to show that precursors of the phase transitions are clearly observed even for small systems. All calculations are done for N=6.



We select three distinct paths that represent various regions of the phase diagram.
Path 1 follows the α=0 line, corresponding to the standard Lipkin model, which exhibits a second-order phase transition at the mean-field level. Path 2 traverses the α=1/√(2) line, associated with a first-order phase transition. Finally, path 3 keeps λ=1 while varying α as the control parameter.
Figures <ref>(a), <ref>(b), and <ref>(c) illustrate the correlation energies along these paths in the ELM phase diagram.

In the large-N limit, the ground state energy per particle displays a change in slope at the critical point, λ_c of Eq. (<ref>). This transition is gradual for α=0 (path 1) due to the second-order nature of the phase transition. For α≠ 0, the transition is more abrupt, particularly for α=1/√(2) (path 2), as a result of a first-order phase transition. Despite the finite size effects in our simulations for a small system with N=6, Fig. <ref> demonstrates precursors of a QPT in the ground state energy around λ_c for paths 1 and 2, depicted in panels (a) and (b). Notably, the slope change near λ_c is smooth for path 1 (panel a), indicative of a second-order phase transition, and more abrupt for path 2 (panel b), suggesting a first-order phase transition. In contrast, no indication of a QPT is observed for path 3 (panel c). It is worth noting that the apparent crossings observed for some excited states in panel (c) are not genuine, but rather avoided crossings, as the mechanism of level repulsion is at play in this scenario.









In addition to the energy of the system, there are other physical properties that reflect the change of phase. 
For the ELM, the expectation value of (S_z+S) in the ground state can be considered as an order parameter that, in the large-N limit, is zero in the symmetric phase, and different from zero in the broken one. Panels a), b) and c) in Fig. <ref> represent this order parameter, normalized to N, for N=6 and the same three paths discussed in the preceding figure. 
For finite particle number, N, the sudden change of behaviour in this order parameter is not exactly at the mean-field critical value of the control parameter.

Also for this observable, precursors for the phase transition are clearly observed in the ground state for paths 1 and 2 in the vicinity of λ_c, Eq. (<ref>). 
In contrast, we do not see change in behavior along path 3, where the system is always in the broken phase. In addition, in Fig. <ref> the change in the order parameter at the critical point is more abrupt for path 2 than for path 1. This is due to the fact that, in the large-N limit, the transition in path 1 is of second order (continuous change in the order parameter at the critical point), while in path 2, it is of first order (discontinuous change in the order parameter at the critical point). 





§ ADAPT-VQE ALGORITHM FOR THE ELM


We have so far discussed the phase diagram of the ELM. We have illustrated the richness of the phase diagram of this model, which is reflected in the properties of even relatively small systems. We now turn to discuss how simulations for such small systems may potentially be performed in quantum devices. 




 §.§ The ADAPT-VQE algorithm 
 

The ADAPT-VQE algorithm <cit.> is a method within the class of Variational Quantum Eigensolver (VQE) approaches <cit.>. VQEs aim at solving variational problems, such as the ground state of quantum many-body systems, by 
parametrizing their 
wave functions in terms of quantum circuits, and minimizing the corresponding energies employing classical algorithms. The hope is that such hybrid methods may provide a quantum advantage over classical simulations <cit.>

ADAPT-VQE  was originally devised to calculate bond dissociation curves in molecules <cit.>, but has since been employed to study other physical systems like the Lipkin model <cit.> or nuclei <cit.>. 
The methods is iterative and starts from an initial reference trial state, |ref⟩. This state is updated by means of encoded operators in a previously-defined operator pool. Thus, at iteration n, the resulting quantum state |n⟩ is given by a product

    |n⟩ = e^iθ_n A_n|n-1⟩ = ∏_k=1^n e^iθ_k A_k|ref⟩.

Here, θ_k with k=1, ⋯, n are n variational parameters
associated to the n operators, A_k. 

The optimal values of these parameters are obtained by minimization of the energy surface resulting from the expectation value of this state with the Hamiltonian. 

At iteration n-1, the algorithm must decide which operator to choose at iteration n. 


The gradient at θ_k=0 can be written in a compact expression as <cit.> 

    ∂ E^(n)/∂θ_k|_θ_k=0 = i ⟨ n | [H, A_k] | n ⟩.

The algorithm selects the operator A_n that produces the largest gradient, and thus is expected to have the optimal descent in the energy surface. Compared to approaches without adaptive gradients, this should allow for a faster approach to the real energy minimum, potentially avoiding plateaus <cit.>. We note that the operators A_k
are selected from a common pool, and may be repeated along the iterative minimization. 


Within ADAPT, quantum circuits can be used to construct the states and measure the commutators. The variational parameters are optimized at each iteration minimizing the total energy surface, on a classical computer, aided by measurements in the quantum computer <cit.>. 
For the purpose of evaluating its performance, we do not operate explicitly with quantum circuits, but rather  simulate the evolution of the algorithm on a classical computer. 


We intend to show that the formalism developed in Ref. <cit.> for the standard Lipkin model is also valid for the ELM. We note this is indeed the case, even though the ELM has a more complex phase diagram including first- and second-order phase transitions.
In order to apply the method, we need to express the collective spins appearing in the ELM Hamiltonian, Eqs. (<ref>) and (<ref>), in terms of individual spins,

    S_+ = ∑_k a^†_k,+a_k,-≡∑_k σ_+^k,  
       S_-= ∑_k a^†_k,-a_k,+≡∑_k σ_-^k, 
       S_z = 1/2∑_k (a^†_k,+a_k,+ - a^†_k,-a_k,-) ≡1/2∑_k σ_z^k.

where σ_±^k ≡ (σ_x^k ± iσ_y^k)/2.




Here, σ_x^k, σ_y^k, and σ_z^k are the Pauli matrices for the k-th site. a^†_k,± and a_k,± are the creation and annihilation operators for a particle in the k-th position with spin projection ± corresponding to the upper or lower levels, respectively. Because the Lipkin model only allows vertical excitations, one can treat the upper and lower states within the same site as independent qubit projections. 
The model can be interpreted in the context of quantum information processing <cit.> without explicitly resorting to the Jordan-Wigner transformation <cit.>, which reduces the size of the Hilbert space by a factor of 2. In this approach, the standard Lipkin model Hamiltonain is written as,

    H_L= (1-λ) ( S+1/2∑_i σ_z^i)- λ/2N( ∑_i<j (σ_x^i σ_x^j - σ_y^i σ_y^j)+∑_i<j (σ_x^i σ_x^j+σ_y^i σ_y^j) ),

which leads to the ELM Hamiltonian

    H_EL   = H_L - λ/N[α^2 (S^2+1/4∑_i<jσ_z^i σ_z^j+S ∑_iσ_z^i) 
          - α( 2S∑_i σ_x^i+1/2∑_i<jσ_x^i σ_z^j+1/2∑_i<jσ_z^i σ_x^j) ] .










In order to keep the circuits as shallow as possible, we restrict the operator pool to one- and two-body operators that act on particles in at most two pairs of levels. The specific form of the one-body operators that we include in the pool 
are:

    G_+^k    =   σ_+^k + σ_-^k = σ_x^k,
    
    
    G_-^k    =    -i(σ_+^k - σ_-^k) = σ_y^k,
    
    
    G_0^k   =   σ_z^k,

while for the two-body operators, with j<k, our pool contains the operator structures:

    T_+^jk   =   σ_+^j σ_+^k + σ_-^j σ_-^k = 1/2 (σ_x^j σ_x^k - σ_y^j σ_y^k),
    
    
    T_-^jk   =    -i (σ_+^j σ_+^k - σ_-^j σ_-^k) = 1/2 (σ_x^j σ_y^k + σ_y^j σ_x^k),
    
    
    U_+^jk   =   σ_+^j σ_-^k + σ_-^j σ_+^k = 1/2 (σ_x^j σ_x^k + σ_y^j σ_y^k),
    
    
    U_-^jk   =    -i (σ_+^j σ_-^k - σ_-^j σ_+^k) = 1/2 (σ_y^j σ_x^k - σ_x^j σ_y^k),
    
    
    V_+^jk   =    (σ_+^j + σ_-^j) σ_z^k = σ_x^j σ_z^k,
    
    
    V_-^jk   =    -i (σ_+^j - σ_-^j) σ_z^k = σ_y^j σ_z^k,
    
    
    V_0^jk   =   σ_z^j σ_z^k.

Note that this is the same pool of operators used for the standard Lipkin model in Ref. <cit.>.





 §.§ Results for the ground state energy


In this subsection, we present a comparison of the exact ELM ground state energy values with the approximate results obtained with ADAPT-VQE. 
We employ the same three paths discussed in Section <ref>, as they provide a representative sample of the physics at play. 

In addition, we analyze the number of iterations needed for the relative error, (E_ADAPT-VQE^gs-E_exact^gs)/E_exact^gs, in the ADAPT-VQE calculation to be less than 10^-5, and how this error changes with the number of iterations.
For simplicity, we take as a reference state the lowest energy Slater determinant of the system. We run the ADAPT-VQE simulation for 400 iterations. Our classical minimization employs the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm 


First, we compare in Table <ref>

the exact results for the ground state energy per particle for N=6 to those obtained with the ADAPT-VQE simulation. The results are presented as a function of the control parameters (λ, α). The same comparison is also shown in Fig. <ref>, where the full orange dots correspond to the ADAPT-VQE results, while the full black line corresponds to the exact ones. Panel a) corresponds to path 1 (α=0), panel b) to path 2 (α=1/√(2)), and panel c) to path 3 (λ=1). It is clear

that the ADAPT-VQE results 
agree closely with the exact results in all the analyzed cases.
The bold figures in Table <ref> highlight the difference between the exact and the ADAPT-VQE results. 
We find a level of agremeent that in all but one case reaches the fifth significant digit. 

This is remarkable, in the sense that we are exploring a relevant fraction of the model parameter phase space. Note that in path 1 a second order QPT is present in the large-N limit, whereas in path 2 a first order QPT shows up.  
In these two cases, we find some minor differences at the end of our ADAPT-VQE simulation.


In order to complete the information for the ADAPT-VQE calculations, we show in Fig. <ref> the convergence of ADAPT-VQE calculations for the ground state energy of the ELM. Panel a) illustrates the number of iterations required for convergence, set conventionally at a relative difference with the exact value below 10^-5, against the system size N for selected points in path 1, α=0. These points include λ=0.1 (blue line-full dots), λ=0.2 (green line-crosses) and λ=0.9 (red line-full squares). Panel b) displays a similar analysis for path 2, α=1/√(2).

Equivalent figures for path 3, λ=1, in which the system is always in the non-symmetric phase, are similar to the λ=0.9 results shown in panels a) and b), independently of the α value.  
From Fig. <ref>, it is clear that the convergence in points corresponding to the non-symmetric phase, λ=0.9, is reached very fast (red lines) independently of the value of α.  
We hypothesize that this is due to the fact that, for all particle numbers, the system has a well-defined deformed minimum. This is true except close to the critical line, which is around the point λ=0.2 (green lines). Here, there are two competing minima and the convergence of a variational method is expected to be slow. This is indeed confirmed by the corresponding green lines in panels a) and b). In any case, the increase in the number of iterations grows approximately linearly with the size of the system.

In Fig. <ref>, we present the relative errors obtained in the ADAPT-VQE calculations with respect to the exact values of the ELM ground state energies for a system size N=6 versus the number of iterations needed. Three selected points in different paths going from symmetric (λ=0.1) to non-symmetric, well-deformed (λ=0.9) shapes are presented. Panel a) corresponds to path 1 (α=0), panel b) corresponds to an intermediate value of α=0.4, and panel c) corresponds to path 2 (α=1/√(2)). 
Again, we find that the faster convergence occurs for λ=0.9 (red lines), where the system has a well-defined, deformed minimum and ADAPT-VQE finds it easily. The most difficult convergence occurs at the region around the critical line, λ=0.2, (green lines) where the competition between a spherical and a deformed minimum 
hampers the convergence of the variational method. The symmetric cases, λ=0.1 (blue lines) are intermediate situations. Regardless of the path, a moderate number of iterations (between 5 and 60) is needed to reach convergence in small systems.

We conclude that the ADAPT-VQE with the pool of operators given in (<ref>-<ref>), as proposed for the standard Lipkin model in Ref. <cit.>, gives an excellent approximation for the ground state energy of the complete parameter space (λ,α) of the ELM.





§ PROPOSAL FOR THE IMPLEMENTATION OF THE ELM MODEL IN A DIGITAL-ANALOG QUANTUM PLATFORM


We have so far discussed the ELM and its potential implementation within a VQE. Our findings indicate that ADAPT-VQE can reach the energy minimum even in situations where the system is close to a phase transition in the large-N limit. 
We can alternatively exploit the distinct dynamical features in different regions of the parameter space to identify the different phases of the system <cit.>. This, however, requires access to the temporal evolution of the many-body model at hand.

In this section, we intend to determine the necessary resources for the implementation of the time evolution operator of the ELM Hamiltonian in a quantum simulator. To this end, let us first rewrite the Hamiltonian of Eq. (<ref>) as


    H_EL = g_0 + [1-(1+α^2)λ] S_z + 2αλ S_x - α^2λ/N S_z^2 - 4λ/N S_x^2 + 2αλ/N (S_xS_z + S_zS_x).

Considering S_xS_z + S_zS_x = (S_x + S_z)^2 - S_x^2 - S_z^2 one gets 

    H_EL = g_0 + g_z S_z + g_x S_x + g_zz S_z^2 + g_xx S_x^2 + g_xz(S_x+S_z)^2,

where the factors of the collective spin operators have been absorbed into the corresponding g_j coefficients.

Equation (<ref>) makes apparent what the necessary blocks for the simulation are. The evolution of the operators S_z and S_x can be implemented in the quantum circuit using high-fidelity Single-Qubit Rotation (SQR) gates. Meanwhile, under a purely digital quantum simulation, entangling terms would require a number of lower-fidelity two-qubit (CNOT) gates on top of the SQR gates, which increases with the system size <cit.>. Luckily, in trapped-ion systems, one can efficiently simulate S_x^2 by means of a single Mølmer-Sørensen (MS) gate affecting all qubits simultaneously, which evolves analogously to S_x^2 <cit.>. Since the error of the two-qubit gates dominates the infidelity of the whole system, this approach offers a significant advantage with regards to purely digital quantum simulations by greatly reducing the required number of two-qubit gates, which is only one, independently of the system size. Note, though, that the fidelity of MS gates in trapped-ions systems still decreases as the number of qubits increases, even if the number of MS gates does not.

Considering the convenience of these MS gates, we propose an implementation based on the recently developed Digital-Analog Quantum Computation (DAQC) paradigm <cit.>. By means of employing analog blocks as the entanglement sources of the system, alongside high-fidelity SQR gates, one can achieve universal quantum computation, exploiting both the robustness of Analog Quantum Simulations and the versatility of Digital Quantum Computers. Fig. <ref> shows a scheme of the proposed quantum circuit.

First, the evolution of a single body term G_ν = g_ν∑_i^N σ_ν^i, where σ_ν^i is the Pauli operator ν∈x,y,z acting upon qubit i, can easily be implemented with a SQR gate R_ν^i(θ) = exp[iθ/2σ_ν^i], with θ = g_ν t, applied to every qubit. So, the terms S_x and S_z only require N single-qubit gates each.



Second, as previously mentioned, making use of the properties of trapped-ion systems, the term S_x^2 can be directly implemented as an analog block via a single MS gate, which is independent of N. The MS gate consists of a bichromatic laser field acting on the qubits, and the parameter can be controlled through the phase and frequency of this laser field. Lastly, the terms S_z^2 and (S_x + S_z)^2 can be implemented in a similar manner, by applying SQRs R_y(π/2) and R_y(π/4), respectively, to all qubits before and after the MS gate. Effectively, this performs a local rotation from the z or x+z axes to the x axis, which can then be evolved analogically with a MS gate, then rotated back to its original axis. This means that all three two-body terms can be implemented using only three two-qubit MS gates (one for each term), plus 4N single-qubit gates to rotate S_z^2 and (S_x + S_z)^2 into the x direction and then back (2N gates each).

Therefore, for each single-body term that is implemented, we need N gates, for a total of 2N (there are two single body terms). For each two-body term that we need to rotate to the x axis, we also need 2N single-qubit gates. Because we have three two-body terms, but one of them already points along the x axis, we only need 4N single-qubit gates. The total number of gates for the circuit would then be 6N single-qubit gates plus three MS gates. However, since not all the terms in the Hamiltonian commute with each other, applying these gates sequentially results in an evolution that differs from that of the target Extended Lipkin Hamiltonian. A common way to deal with this error is to employ the first-order Lie-Trotter-Suzuki formula <cit.>:

    U(t) ≈( ∏_k e^-i H_k t / n_T)^n_T = U_T(t,n_T),

where n_T are the so-called Trotter steps and H_k are the k non-commuting terms that conform the Hamiltonian, i.e, H = ∑_k H_k. For this particular case, the terms are chosen as follows:

    H_1 =     g_zS_z + g_zzS_z^2, 
    
        H_2 =     g_xS_x + g_xxS_x^2, 
    
        H_3 =     g_xz (S_x + S_z)^2.

One can approximate the exact evolution by performing the evolution in time t in many small steps  of time Δ t, called Trotter steps. 
The approximation becomes exact in the limit n_T →∞ (Δ t → 0), which means that it can be made arbitrarily precise by increasing the number of Trotter steps. Note, though, that, in order to do so, one must increase the depth of the circuit, aggregating the infidelities of the quantum gates in the circuit for each run. Therefore, it is important to study, case by case, the experimental error introduced with each Trotter step versus the digital one due to the Trotter approximation. 

Fidelities of up to 99.9999% have been reported for single-qubit gates <cit.>, while fidelities of up to 99.9% have been achieved for MS gates in trapped-ion systems <cit.>.
Since the infidelities of individual gates are small enough and given the resource quantification we just presented, it is reasonable to assume that the total fidelity of the system will approximately be

    f = [ (f_s)^6N(f_MS)^3 ] ^n_T,

where f_s and f_MS are the fidelities of single-qubit and MS gates, respectively.


One can define the fidelity of the Trotter approximation for a particular initial state |ϕ(0)⟩ as

    F = |⟨ϕ (0) | U_T(t,n_T) U(t) | ϕ(0) ⟩|^2.

The values of the fidelity of the Trotter evolution are shown in Fig. <ref> for different time values and Trotter steps in a system of size N=6. In panel a), we find that, for a low number of Trotter steps, the fidelity decreases at a very fast rate with time. This suggests that the digital error will dominate over the gate error, except for short or very short times. In contrast, a high number of Trotter steps will allow longer evolution times without incurring strong infidelities due to digital errors, but with a high gate infidelity that will dominate the total error. Similarly, in panel b) of Fig. <ref>, we show how the fidelity improves with n_T for different values of time. For a given amount of evolution time t, there is a value of n_T beyond which there is no substantial improvement of the fidelity.

This leads to a compromise in fidelity, as decreasing the errors from one source implies increasing the errors from the other. This often requires finding a balance between the precision of the Trotter expansion and the precision of the physical implementation, but, as it will be shown in the following section, for certain tasks, such as finding the phase diagram of a quantum system, ML approaches can properly work with a reduced number of Trotter steps, saving computational resources and thus minimizing the total gate infidelity.



§ UNVEILING THE ELM PHASE DIAGRAM THROUGH A MACHINE LEARNING ANALYSIS




We now discuss specifically an ML implementation to identify the different phases of the ELM.
In Section <ref>, we have discussed the use of the expectation value of ⟨ S_z ⟩ + N/2 for the ground state as an order parameter. In Fig. <ref>, one can clearly observe the precursors of the QPT for a relatively small system of size N=6, far from the large-N limit. In this section, we will present a different way of identifying the critical points of the QPT. We will employ ML techniques that can enhance our ability to precisely classify the quantum phase of the system. First, we will exploit our previous theoretical knowledge of the phase diagram using supervised ML algorithms; then, we will propose strategies to gather new information in systems where the phase diagram is unknown with the aid of unsupervised ML methods.



 §.§ Supervised quantum phase classification


A common way of identifying a QPT is by measuring an order parameter, which 
for the ELM Hamiltonian can be (⟨ S_z ⟩ + N/2)/N (see  Fig. <ref>). Nevertheless, in practice, this implies  a priori knowledge of the ground state of the system, which may be a big challenge in certain many-body systems e.g., in the Heisenberg or the Ising models. Therefore, other less costly methods may be of interest, especially if it is not necessary to obtain the ground state of the system. 
In Refs. <cit.>, some of us proposed to use the time evolution of an appropriately chosen operator as a proxy of the phase of the system. In particular, the following time-dependent function was used 

    C_ν(i,j,t) = ⟨Φ_0 |U(t)^†σ_ν^i σ_ν^j U(t)| Φ_0⟩ - ⟨Φ_0| U(t)^†σ_ν^i U(t)| Φ_0⟩⟨Φ_0|U(t)^†σ_ν^j U(t)|Φ_0 ⟩,

where σ_ν^i corresponds to the Pauli matrix on the axis ν = { x,y,z } over the i-th qubit and U(t) to the evolution operator. Note that this function depends on the considered state, |Φ_0⟩, that should not be an eigenstate of the system. 
We will also employ  the time evolution of the average value of S_z,

    S_z(t) = ⟨Φ_0 | U(t)^† S_z U(t) | Φ_0 ⟩,

which we intuitively expect to be sensitive to the phase of the system. 
However, inferring the phase is not always a trivial task. 
Among other things, one needs to find a rule that connects the time evolution of these observables with the phase of the system.



A very appealing way to establish the relationship between the phase and the time evolution is to use supervised ML based on our knowledge of the phase diagram (see Fig. <ref>). For example, one could use labeled data to train a Neural Network as a classifier, which would then be able to recognize the phase in which the system is in, based only on the measurement of its order parameter <cit.>. 

Furthermore, ML allows for the use of less intuitive or visually apparent methods. This expands the set of tools available to study QPTs, using different operators as probes of the QPT or even initial states that are not the ground state, which is our objective.

There is no single ML algorithm that consistently performs classification tasks better than any other algorithm <cit.>. In other words, multiple different approaches can be equally viable. In our specific case, we choose to build and train a Convolutional Neural Network (CNN) <cit.>. This choice is motivated by the well-known success of CNNs in recognizing shapes and patterns as image classifiers and other visual recognition tasks <cit.>, as well as by its application in previous similar works <cit.>. The CNN consists of a standard architecture: 3 Convolution-Relu-Pooling layers with 32 kernels of size 3, plus a set of 4 fully-connected (dense) layers consisting of 512 neurons each. We also introduce dropout layers as a measure against data overfitting. 

To train the CNN, a dataset is needed. Since the ELM is exactly solvable, the dataset is created by classically computing the evolution of the expected value of the desired operators, of Eqs. (<ref>) or (<ref>),

for a system with N=6 particles.






These dynamical evolutions are computed for values of the control parameters λ∈ [0,1] and α∈ [0,1], in 0.01 steps, yielding a grid of 10201 data points. Each evolution is computed for a conveniently adimensionalized time t ∈ [0,4] with steps of size 0.04. 
Then, we label our simulations as either 0, if they correspond to the symmetric phase, or 1, if they correspond to the broken-symmetry phase. 
We remove the simulations of the 
points corresponding to the three paths mentioned in Sect. <ref> from the dataset. These, alongside a randomly selected 20% of the remaining points, are kept apart for testing purposes. The remaining data is then used for training and validation of the training. The cross-validation strategy employs a randomly selected 20% of the non-test data at each iteration (epoch) of the training. This validation data is used to confirm that there is no overfitting, sudden forgetting or other anomalies during the training process, as well as to optimize the hyperparameters of the CNN, such as the number of layers and the number of neurons per layer. The test dataset is used exclusively to compute the results reported here.  

First, we use S_z(t) in Eq. (<ref>) to train the CNN with the labeled data. The output of the network is not a discrete classification; instead, it provides the probability, P_BS, 
for the corresponding point in the phase diagram to belong to the broken symmetry (BS) phase. Ideally, this probability should be close to 0 in the symmetric phase, close to 1 in the broken phase, and in between those values when close to a critical point. 

Once the CNN has been trained, the performance is analyzed with the 
dataset that was not considered in the training. 

Figure <ref> depicts the probability P_BS for the testing dataset using the evolution of the operator S_z (t), a system size N=6, an initial state |Φ_0 ⟩ = 1/√(3)(|m_z=-N/2⟩ + |m_z=0⟩ + |m_z=N/2⟩ ) and control parameters α∈ [0,1/√(2)] and λ∈ [0,1]. We show in panels a), b), and c) the results for the three relevant paths, considering the exact evolution for training and predicting. Other than the inaccuracy at the extremes of paths 1 and 3, for values λ=1 and α=0, respectively, the rest of the predictions are precise and mostly with complete certainty. Paths 1 and 2 show a very accurate prediction of the critical point (marked with the red dotted line), as P_BS suddenly jumps from 0 to 1 precisely in this area, despite working with a relatively small system size. Panels d), e) and f) (including the insets which provide a zoomed view) give a more general view of the quantum phase diagram classification (using also S_z (t) and N=6) across the entire parameter space, which is, once again, very precise. In panel d), the system is trained with the exact evolution operator (U). Symmetric and broken phases are well determined and the critical line (yellow line) is nicely reproduced, although a small part of the BS phase is predicted as symmetric (with low certainty). In panels e) and f), the Trotter evolution is considered, with n_T=4 and n_T=1, respectively. It is remarkable that these results are of the same quality as the exact ones, specially for n_T=1.  This fact is particularly relevant as it means that it is not necessary to correct the Trotter infidelities, saving computational resources and, thus, reducing the total infidelity incurred by lowering the number of required quantum gates.



Second, we use the correlation operator of Eq. (<ref>) to train the CNN to predict the phase diagram. This option is motivated by the results obtained in Refs. <cit.>.
Choosing different correlation functions can yield slightly different results and, due to the nature of deep learning algorithms, it is usually not possible to intuitively guess which one may be the best option. We stress that, while the following results were achieved using the  C_z(1,2,t) correlation function, analogous results were achieved with all possible combinations of i,j. The correlations C_y(i,j,t) and C_x(i,j,t) might also be viable probes. We follow the same procedure just described for S_z, replacing its evolution for the time evolution of C_z(1,2,t) as input of the CNN. In addition, the initial state |Φ_0 ⟩ = | ↑↓↑↓↑↓ ⟩ was arbitrarily chosen. Note that, since σ_z^i acts on a particular qubit, one cannot compute the correlation function in the collective spin basis. 

The results for this second choice of correlator are shown in Fig. <ref>.
We reach the same conclusions as for S_z(t), but with an improvement in the ability of the CNN to correctly classify the quantum phase of the system and precisely identify the critical point. For the three selected paths, the CNN only shows uncertainty for the points closest to the critical points in paths 1 and 2 (λ = 0.20 and λ = 0.18), which is the best possible guess. 

This ML approach to phase classification has another advantage, which is its robustness against errors introduced before the training. This is generally true for properly constructed CNNs <cit.>, but they become particularly robust for data such as the one at hand. The evolution of C_z(1,2,t) at each (λ,α) value has a particular shape, and the CNN learns the features of these shapes to classify them. The variation of the parameter λ has a very prominent effect on the shape of this evolution, while the parameter α has a relatively small effect in comparison. With the addition of  noise, data points at different values of α for the same value of λ become difficult to distinguish, but the CNN is still able to, overall, differentiate points at different values of λ. Moreover, the QPT depends mostly on λ, so even if the details on α are lost, the CNN can still classify most of the phases correctly, albeit with a shift in the exact location of the critical line. 

Having said that, one can train the Neural Network with the same kind of noisy data that one is expecting to predict. This allows the CNN to, in a way, learn the noise pattern, recovering a precise prediction of the QPT. Both of these issues are shown in the predictions of Fig. <ref>, where a random Gaussian noise with a standard deviation σ = 0.2 was added to the input data. Panel a) shows the results for a CNN trained with noiseless data, when predicting the phases of noisy data. We find that the predicted phase transition shifts, and the predictions of the symmetric phase enter the area of the BS phase, but the overall phase classification remains quite accurate. In contrast, panel b) shows a CNN trained on noisy data. Here, the misclassification issue is fixed and the phase transition is once again accurately predicted. We note, however, that infidelities in quantum gates  and other sources of errors in quantum simulators are biased, unlike random noise. It remains unclear whether the infidelities of a real quantum simulator would affect the predictions of the CNN, although it is safe to assume that significantly smaller errors will have a significantly smaller impact.





 §.§ Supervised classification with partial information


The main and most obvious disadvantage of supervised ML algorithms is that they require computing or labeling the data with prior knowledge of the correct classification. For that reason, one of the main goals of many supervised ML algorithms is to achieve a good enough generalization using a reduced dataset. 

We explore the generalization ability of our setup as follows.
First, we train our CNN only with data corresponding to the Lipkin model of Eq. (<ref>), that is, with α=0 in the ELM Hamiltonian (<ref>). Then, we try to predict the quantum phase diagram for the ELM, with α≠ 0. Note that the regular Lipkin Hamiltonian exhibits only a second-order QPT, while the ELM presents a first-order QPT for values of α≠ 0 <cit.>.

Figure <ref> shows the predicted quantum phase diagram using this approach. Since the CNN had no information about points with α≠ 0, we may expect  the predictions to become inaccurate as we move towards large values of α. However, the CNN predictions remain surprisingly accurate for large values of α. We take this as a sign of the generalization capabilities of CNN's. We stress that these results are not trivial. We are able to predict a first-order QPT in a fundamentally different Hamiltonian, with only data from the second-order QPT of the known Lipkin model. The prediction with only partial information has a slightly smaller categorical accuracy (98.86%) when predicting the quantum phase of the system, compared to the results of Fig. <ref> (99.85%). In other words, using the complete quantum phase space for training provides marginally more accurate predictions. The difference is, however, not easily observed in Fig. <ref>. 





 §.§ Unsupervised classification



Finally, we are interested in correctly predicting the phase diagram of the system without any prior theoretical knowledge about it. To achieve this, the use of unsupervised methods, which do not employ labeled training data, becomes necessary. Of course, with this approach, new challenges arise. The lack of labeled data for training means that these algorithms rely on the definition of some sort of similarity between the data.  


As will be shown below, unsupervised ML will allow to distinguish points near the critical area from those that lay well apart. This will provide a way to deduce the quantum phase diagram, but requires some previous knowledge of the model at hand, i.e., one needs to know where different phases exist. Of course, for more complex quantum phase diagrams with more than a single QPT, a different approach is likely to be necessary. 


There are many possible unsupervised ML algorithms that could be used, namely, clustering, factor analysis, principal component analysis, and embedding techniques.  
Analogously, different observables have been studied as input for the unsupervised algorithms, concluding that the evolution of the functions C_z(1,2,t) and S_z(t), which were used in the supervised ML, are unsuitable for the unsupervised cases, being impossible to obtain conclusive results. This was expected as part of the limitations of unsupervised ML, which cannot easily learn specific patterns or features to classify. 



We consider a different alternative and use as input the different correlation functions that arise from all possible combinations of i,j, and ν in Eq. (<ref>). This means that, in a system of size N=6, for a single data point, i.e. a single value of α and λ for the Hamiltonian, 45 correlation functions are computed (note that C_ν(i,i) = 0 and C_ν(i,j) = C_ν(j,i)). We use these correlations for a clustering algorithm known as Fuzzy C-means <cit.>, which is a variation of the popular K-means clustering algorithm <cit.>. 

K-means <cit.> is a general-purpose clustering algorithm in which, given a set of n samples X⃗ = {x⃗_i, …, x⃗_n } and a number of clusters K, cluster centroids C⃗ = {c⃗_1, …, c⃗_K } are computed to minimize the  objective function,

    f(C⃗) = ∑_i^n ∑_j^k w_i,j (|| x⃗_i - c⃗_j ||)^2,

where w_i,j∈{ 0, 1 } is a partition matrix that tells us whether sample i belongs to cluster j (w_i,k=1) or not (w_i,k=0). We implicitly assume that x⃗_i and c⃗_j correspond to multidimensional vectors.
Once the number of clusters has been fixed, cluster centroids are initialized, commonly either randomly or through some initialization function that may speed up the convergence, and each sample is assigned to the closest centroid. Then, the mean position of the samples assigned to each cluster defines the position of the cluster centroid. The samples are once again assigned to the new closest centroids and the process is repeated until convergence is reached. 

Fuzzy C-means is very similar to the K-means algorithm but 
the partition function w_i,j∈ [0,1] is given by

    w_i,j =1/∑_j^'^k ( || x⃗_i - c⃗_j ||/|| x⃗_i - c⃗_j^' ||)^2/m-1,

where m ∈ (0,∞) is a parameter that controls the “fuzziness” of the clusters. 
These weights are a measure of the degree of membership of the sample i to the cluster j and now appear in the partition function, Eq. (<ref>), to the power of m, as (w_i,j)^m. Note that, contrary to K-means clustering, each sample can now be partially a member of several clusters simultaneously, to a certain degree given by the weight w_i,j. In the limit m→ 1, the membership converges to either 0 or 1 and we recover the K-means algorithm. The bigger the value of the fuzziness parameter m, the larger the overlap between clusters and, as a matter of fact, in the limit m→∞, all samples belong to all clusters equally. A benefit of this fuzzy approach is that it allows to distinguish the points that are more representative of each cluster and those that are in-between.



In Fig. <ref>, the value of the cluster membership for each point of the phase diagram is presented, using the Fuzzy C-means algorithm with k=2 clusters and a fuzziness parameter m=4. Since  k=2, w_i,1 = 1 - w_i,2. Note that this algorithm is not capturing the quantum phases of the system, but whether a point is close to the critical area or the phase is well defined, without distinguishing between the symmetric or the broken phase. 
Thanks to this kind of soft clustering, we can identify the points that have a higher (or lower) degree of membership as the points closest to the QPT of the system. In fact, the dark green line in Fig. <ref> corresponds to the maximum value of the degree of membership for each value of α, which is close to the theoretical critical points for each value of α (light green line). Notably, for a system with a single QPT and two phases, this prediction of the QPT line allows for a relatively precise reconstruction of the quantum phase diagram, even if the unsupervised ML algorithm is unable to explicitly distinguish the corresponding quantum phases.

In general, predicting the location of the QPT with unsupervised ML methods can become a difficult task, since the system is expected to behave very similarly when measured on two points if both of them are near the critical point, even though they may correspond to different phases. This is specially true for small system sizes, like the ones used on the supervised ML methods before, where due to finite size effects the QPT becomes smooth instead of sharp.



§ SUMMARY AND CONCLUSIONS


In this work, we have explored the possibility of implementing the extended Lipkin model (ELM) in a analog-digital quantum simulator. First, we have presented the model and shown the connection with the nuclear interacting boson approximation (IBA). We have paid special attention to the connection between both models in the large-N limit where the phase diagram is identical in the two cases. In both models, spherical and deformed phases exist and, depending on the trajectory followed in the phase space, first- and second-order quantum phase transitions appear. Furthermore, we have used the ADAPT-VQE algorithm to calculate the ground state energy. To this end, the Jordan-Wigner mapping is not needed, because the ELM in the spin representation is directly mapped into the quantum simulator. The results are very close to the exact ones, but the convergence strongly depends on the position of the Hamiltonian in the parameter space. Next, a proposal for the implementation of the time evolution operator of the ELM on a quantum platform based on the Digital-Analog Quantum Computation (DAQC) paradigm is presented. Calculations of the expected fidelities are also presented.
Finally, we present a Machine Learning study to identify the phases of the ELM. This study proposes, as proxies to identify the phase of the system, the time evolution of ⟨ S_z (t) ⟩ and a two-point correlator C_z(1,2,t), Eq. (<ref>).  Our analysis has three different flavours. We first perform a supervised classification task with full data. Then, we carry out the same exercise with partial, noisy information. Lastly, we also use an unsupervised classification exercise based on the C-means algorithm. In all cases, the results obtained capture the position of the different phases and can predict rather precisely the critical area. 

These studies lay the groundwork for the implementation of more intricate nuclear structure models on quantum platforms. In particular, the use of the ADAPT-VQE method to calculate the ground state energy or even the excited states opens the possibility of dealing with nuclear systems with large Hilbert spaces <cit.>.

The proposed supervised method for determining the phase diagram of the system (see Sec. <ref>) requires the calculation of the evolution of a well chosen observable for an arbitrary wave function that does not correspond to the ground state of the system. This is a very relevant fact, because the calculation of the ground state in medium and large size many-body systems can be very challenging. The method nonetheless provides accurate results using very simple states and taking advantage of Machine Learning techniques. Moreover, the use of unsupervised methods (see Sect. <ref>) has provided very promising results allowing to locate the position of the critical area in the phase diagram with minimal biases.



§ ACKNOWLEDGMENTS

This work was partially supported 

through the projects PID2022-136228NB-C21, PID2022-136228NB-C22, PID2020-114687GB-I00, PID2020-118758GB-I00 funded by MCIN/AEI/10.13039/50110001103 and by “ERDF A way of making Europe”; 
“Ramón y Cajal" grant RYC2018-026072 funded by MCIN/AEI/10.13039/5011­00011033 and by  “ESF Investing in your future”; the “Consolidación Investigadora” Grant CNS2022-135529 
funded by MCIN/AEI/10.13039/­501100011033 and by the “European Union NextGenerationEU/PRTR” and the “Unit of Excellence María de Maeztu 2020-2023" award to the Institute of Cosmos Sciences, Grant CEX2019-000918-M funded by MCIN/AEI/10.13039/501­100011033.

This work has also been financially supported by the Ministry for Digital Transformation and of Civil Service of the Spanish Government through the QUANTUM ENIA project call - Quantum Spain project, and by the “European Union NextGenerationEU/PRTR” within the framework of the “Digital Spain 2026 Agenda”. Resources supporting this work were provided by the CEAFMC and the Universidad de Huelva High Performance Computer (HPC@UHU) funded by ERDF/MINECO project UNHU-15CE-2848. 


