Variational quantum simulation of long-range interacting systems

Current quantum simulators suffer from multiple limitations such as short coherence time, noisy operations, faulty readout and restricted qubit connectivity in some platforms. Variational quantum algorithms are the most promising approach in near-term quantum simulation to achieve practical quantum advantage over classical computers. Here, we explore variational quantum algorithms, with different levels of qubit connectivity, for digital simulation of the ground state of long-range interacting systems as well as generation of spin squeezed states. We find that as the interaction becomes more long-ranged, the variational algorithms become less efficient, achieving lower fidelity and demanding more optimization iterations. In particular, when the system is near its criticality the efficiency is even lower. Increasing the connectivity between distant qubits improves the results, even with less quantum and classical resources. Our results show that by mixing circuit layers with different levels of connectivity one can sensibly improve the performance. Interestingly, the order of layers becomes very important and grouping the layers with long-distance connectivity at the beginning of the circuit outperforms other permutations. The same design of circuits can also be used to variationally produce spin squeezed states, as a resource for quantum metrology.

Variational Quantum Eigensolver (VQE) is one of the most established VQAs for generating the low-energy eigenstates, specially the ground state, of a many-body Hamiltonian [23,[60][61][62]. In this algorithm, a parameterized quantum circuit is used to generate a complex quantum state from a simple input wave function. Then the average energy of the desired Hamiltonian is measured at the output of the parameterized quantum circuit. This measured average energy is fed into a classical optimizer to be iteratively minimized. At each iteration the classical optimizer provides a new update for the parameters of the quantum circuit and a new set of measurements on the quantum simulator provides a new estimation for the average energy. As the average energy reaches its global minimum the output of the quantum circuit simulates the ground state of the system. VQE has also been generalized for simulating excited states through addition of penalizing terms to the cost function [60,[63][64][65], subspacesearch VQE [66] and exploiting symmetries [48,[67][68][69][70][71][72][73][74]. The complexity of any VQA, including the VQE, can be quantified through both quantum and classical resources which are required to successfully accomplish the target task. Quantum resources can be evaluated through the minimum depth of the quantum circuit. Since the most challenging ingredient of any quantum circuit is the twoqubit entangling gate, e.g. controlled-not gate, one can naturally use the number of such gates as a measure of quantum resources. The classical resources are determined by the complexity of the optimization. This can be quantified by the number of tunable parameters in the circuit and the number of iterations which is required for convergence of the algorithm.
Long-range interactions are very common and most fundamental forces in nature are naturally long-ranged. In these systems the strength of interaction between two particles decays as 1/d α , where d is the distance and α is an exponent which controls the strength of interaction. Smaller values of α make the interaction more long-ranged. Coulomb (with α = 1), dipole-dipole (with α = 3) and van der Waals (with α = 6) interactions are just a few examples of long-range forces in nature. The presence of long-range interaction in a many-body system can lead to the emergence of several interesting phenomena [75][76][77][78][79][80][81][82] with rich phase diagrams [83]. In particular, the generation of spin squeezing [84,85] through longrange interactions has both fundamental and practical implications. Fundamentally, they can be used for testing the foundations of quantum mechanics [86,87]. From a practical perspective, they can be used as resource for metrology purposes to achieve quantum enhanced precision [88].
An interesting subject to explore is the possibility of VQE simulation for generating the ground state of a longrange interacting system [14,16,89] and spin squeezed states [54,90]. Thanks to the presence of long-range couplings, the ground state of such systems shows significant correlation between distant particles and more multipartite entanglement. So far, analog VQE simulations have been used for generating the ground state of long-range interacting systems [14,16]. However, in such simulations the quantum simulator should already support the same type of long-range interaction between its qubits which may not be feasible. It is thus highly desirable to implement such simulation on digital quantum simulators, as universal computing machines rather than non-universal analog ones. Digital quantum simulators, in NISQ era, come with different levels of qubit connectivity and thus limited interactions between their qubits. In some physical setups, such as superconducting systems [6][7][8][9][10][11][12][13] or cold atoms in optical lattices [21][22][23], the connectivity is determined by the geometry of the simulator and is normally restricted to nearest neighbors. In other platforms, such as ion-traps [14][15][16] and Rydberg atoms [17][18][19][20], the connectivity is more versatile and in principle interactions can be induced between any pair of qubits. A natural open problem is whether the flexibility in qubit connectivity can help digital VQE simulation of the long-range interacting systems. Although, one may expect that longer-distance qubit connectivity should enhance the performance of the VQE simulation of long-range interacting systems, no demonstration still exist to quantitatively show this effect.
In this paper, we address this issue through exploring the VQE simulation of long-range XY model in transverse field. In this system, one can control the anisotropy γ to change the Hamiltonian from Ising in transverse field to anisotropic XX model. Moreover, by tuning the parameter α one can control the strength of long-range interactions. In addition, by controlling the ratio between the exchange interaction and the transverse field one can vary the phase of the system. Therefore, these three control parameters produce a large class of models with a rich phase diagram to be explored. We perform digital VQE simulation for the ground state of the system across its phase diagram assuming different levels of qubit connectivity. We find that the VQE simulation gets harder as the interaction is more long-ranged, i.e. smaller α. In addition, around the quantum phase transition achieving higher fidelity is more challenging and requires more iterations for the convergence of the classical optimizer. Our results show that having connectivity between more distant qubits not only improves the fidelity but also speeds up the convergence of the classical optimizer. Interestingly, grouping the circuit layers with long-distance connectivity at the beginning of the circuit can improve the VQE simulation in terms of both achievable fidelity and required optimization iterations.
The structure of the paper is as follows. We start by giving a brief review on VQE in Section II, showing the important ingredients of such algorithm. Then, in Section III, we introduce the long-range XY Hamiltonian which contains three individual parameters with which one can tune the long-range interaction, the ratio between the exchange interaction and the strength of the transverse field, as well as the anisotropy to change the model from Ising in transverse field to anisotropic XX model. In Section IV, we provide VQE simulation for the ground state of the long-range Ising model in transverse field using different levels of qubit connectivity. We further explore the impact of qubit connectivity in Section V. In Section VI, we provide VQE simulation results for the ground state of various models to demonstrate the generality of our strategy for circuit design. We then show that how our quantum circuits can be used for generating spin squeezing in Section VII. Finally, we conclude our results in Section VIII.

II. VARIATIONAL QUANTUM EIGENSOLVER
Here, we briefly review the VQE algorithm [23] for generating the ground state of a many-body system on a quantum simulator. Those who are familiar with the VQE algorithm can skip this section. A parameterized quantum circuit implements an operation described by a unitary operator U(ϕ), where ϕ = (ϕ 1 , ϕ 2 , · · · , ϕ L ) represents the tunable parameters of local rotations in the quantum circuit. For a given initial state |Ψ 0 ⟩, the output of the quantum circuit is given by |Ψ(ϕ)⟩ = U(ϕ)|Ψ 0 ⟩. In this algorithm, the output of the parameterized quantum circuit will be |Ψ(ϕ)⟩. A cost function C(ϕ) is then measured directly as the average of an observable H at the output of the quantum circuit. The circuit parameters ϕ are updated iteratively via a classical optimizer to minimize the measured cost function C(ϕ). The procedure stops once an optimal set of circuit parameters ϕ=ϕ * is reached. (b) The schematic picture of the impact of circuit depth on the VQE algorithm. For a very shallow circuit, the span set |Ψ(ϕ)⟩ will not be able to include the target state. Therefore, the final output of the VQE algorithm |Ψ(ϕ * )⟩ cannot have a good approximation over the target state. As the depth of the quantum circuit increases, the span set |Ψ(ϕ)⟩ will eventually contain the target state, so that the target state can be well approximated if the global minimum is reached during the optimization.
By varying ϕ, the quantum state |Ψ(ϕ)⟩ can span some part of the Hilbert space. As the quantum circuit gets deeper, and thus the number of parameters L increases, the spanned part enlarges until eventually it covers the whole Hilbert space. In the NISQ era, the goal is to keep the circuit as shallow as possible and yet being able to approximate the target state for a set of parameters ϕ * . Designing such an optimal quantum circuit is a nontrivial task and highly depends on the underlying problem. For instance, for simulating the ground state of a many-body system, a circuit based on adiabatic evolution becomes extremely deep, well beyond the capability of currently available NISQ simulators [48]. Variational quantum algorithms have been developed to achieve the required unitary operation on a shallow circuit. This is accomplished through dividing the complexity between a parameterized quantum circuit and a classical optimizer. As schematically shown in Fig. 1(a), the output of the quantum circuit will be |Ψ(ϕ)⟩. The target state has to be formulated variationally as a minimum of a cost function C(ϕ) which is directly measured as the average of an observable H at the output of the quantum circuit such that C(ϕ) = ⟨Ψ(ϕ)|H|Ψ(ϕ)⟩. In the VQE algorithm, the observable H should be the Hamiltonian of the system and the cost function becomes the average energy [23]. By classically optimizing the cost function C(ϕ) with respect to parameters ϕ, one can iteratively reach to an optimal set of parameters ϕ = ϕ * . Provided that the span set |Ψ(ϕ)⟩ includes the target state, i.e. the ground state |GS⟩, then the output of the quantum circuit |Ψ(ϕ * )⟩ well approximates the target state, namely |Ψ(ϕ * )⟩ ≈ |GS⟩. It is worth emphasizing that the optimization of the cost function C(ϕ) can only generate the ground state if the set of states spanned by |Ψ(ϕ)⟩ includes the target state |GS⟩. For a very shallow circuit, this may not be the case. By increasing the depth of the circuit, the span set |Ψ(ϕ)⟩ will eventually be large enough to contain the target state |GS⟩. The impact of the circuit depth on the span set |Ψ(ϕ)⟩ is schematically described in Fig. 1 The required quantum resources can be quantified by the depth of the circuit. On the other hand, the most complex ingredient of a quantum circuit is the two-qubit gates, such as controlled-not. Therefore, one can simply use the number of two-qubit gates R Q as a quantification for quantum resources. For classical resources, we have to focus on the complexity of the optimization which depends on both the number of parameters L and the number of iterations needed to converge to the minimum of the cost function. Therefore, a natural definition for classical resources can be where #Iteration represents the number of iterations that a classical optimizer needs to converge. The classical resource CR also quantifies the total number of measurements that one has to perform to fulfill the whole process. More precisely, the total number of measurements that one has to perform on the quantum simulator for the whole VQE simulation is twice CR. The factor 2 is because for each parameter the gradient descent demands two successive measurements. Unlike R Q which only depends on the number of qubits and the number of circuit layers, CR highly depends on the choice of the classical optimizer and initialization of the parameters ϕ. This is due to the fact that the number of iterations for the convergence of the optimizer depends on these two factors. In this paper, we use the gradient-based L-BFGS algorithm as our classical optimizer [91] and in order to have a sensible estimation for the number of iterations we repeat the procedure for 50 different random initialization and take the average required iterations for computing CR. Note that the results do not depend on the choice of optimizer, however the convergence speed of the L-BFGS algorithm is normally faster than other gradient-based optimizers, such as Adam [92]. This can be attributed to the ability of L-BFGS to approximate the inverse Hessian matrix, incorporating second-order information that enables more informed and efficient up- dates during each iteration. In contrast, first-order methods, including Adam, depend solely on the gradient or first-order information for their updates. The faster convergence means that fewer iterations are needed for the optimization part, which results in the reduction of classical resources. Finally, it is worth emphasizing that the convergence of the VQE algorithm with average energy as a cost function can only be successful if energy gap ∆E between the ground and the first excited state is larger than the standard deviation of energy, namely ∆E > ⟨H 2 ⟩ − ⟨H⟩ 2 .
If this condition is not satisfied, the classical optimizer cannot discriminate the two eigenstates and thus the VQE outcome will be an arbitrary superposition of them. One way to address this issue is to add extra terms to the cost function, based on the Hamiltonian symmetries whose values are different for the two eigenstates [73]. This will be further clarified in the following sections.

III. MODEL
We study the VQE simulation for generating the ground state of long-range interacting system which is described by the Hamiltonian where H 0 denotes the single particle interaction part of the Hamiltonian, e.g. interaction with external fields, and H k (with k > 0) represents the interaction between all qubits with distance k. In this paper, we consider long-range XY model in transverse field such that where N represents the number of qubits, γ is the anisotropy, H k (with k≥1) indicates the interaction between qubits and 0≤θ≤π/2 controls the phase of the system as cot(θ) is the ratio between the exchange interaction and the transverse field. The strength of long-range interaction is controlled by exponent α. The Hamiltonian in Eq. (3) covers a wide range of interactions as the parameters θ, γ and α vary. For instance, α = 0 represents a fully connected graph in which all the qubits interact equally and as α increases the interaction tends to become short ranged such that in the limit of α→∞ the Hamiltonian (3) becomes the conventional nearest neighbor XY model which can be solved analytically through Jordan-Wigner transformation [93]. Moreover, by controlling the anisotropy γ, one can change the Hamiltonian from Ising in transverse field, i.e. γ = 1, to anisotropic XX model, i.e. γ = 0. For any given α the long-range XY Hamiltonian (3) has a second order quantum phase transition in its ground state |GS⟩ at a special value of θ = θ c which separates a paramagnetic phase from an antiferromagnetic phase. The paramagnetic phase (θ < θ c ) is characterized by all spins aligned in the transverse field direction. For shortrange interacting systems, i.e. large α, the antiferromagnetic phase can be characterized by staggered magnetization at a proper direction depending on γ [93]. In the case of long-range, characterization of the antiferromagnetic phase becomes more complex [83]. As for the case of nearest neighbor interaction, i.e. α→∞, the quantum phase transition is well-known to be at θ c = π/4. For θ ≤ θ c the system is in a paramagnetic phase with a unique ground state aligned by the external magnetic phase in z direction. For θ > θ c the system is in an antiferromagnetic phase which becomes gapless in the thermodynamic limit. The Hamiltonian (3) supports a special symmetry with Z = ⊗ N k=1 σ z k such that [H, Z] = 0. Consequently, each eigenstate of the Hamiltonian has an eigenvalue of ±1 with respect to Z. In the antiferromagnetic phase, where the energy gap is small (in the thermodynamic limit it becomes zero), this symmetry can be used to enhance the quality of convergence for the VQE algorithm to distinguish between the two lowest energy eigenstates. In particular, in this paper, we focus on the ground state with +1 eigenvalue, namely Z|GS⟩ = |GS⟩.
In Ref. [83], the phase diagram of the system for γ = 1 (i.e. transverse Ising) as a function of θ and α has been explored through entanglement analysis using matrix product states for large system sizes N . We perform similar analysis to short systems that we use for our VQE analysis to see how the behavior of the system approaches its thermodynamic limit. Without loss of generality, we also fix γ = 1 and compute the ground state of the system for given values of α and θ. One can trace out half of the qubits from the right side of the system to obtain the reduced density matrix of the left side as ρ L = T r R |GS⟩⟨GS|, where T r R denotes the partial trace. The entanglement between the two halves of the system is then quantified through the von Neumann entropy S V = −T r [ρ L log(ρ L )]. In Fig. 2 we plot S V as a function of θ and α for the long-range XY Hamiltonian with system size N =12. The von Neumann entropy peaks at the criticality and thus can be used as an indicator for the quantum phase transition. Similar results can be found for other γ's.

IV. DIGITAL VQE SIMULATION OF THE GROUND STATE
In this section, we focus on digital VQE simulation of the ground state of the Hamiltonian (3) across its phase diagram as θ and α vary. In order to prevent the complexity of degenerate ground states in the antiferromagnetic phase, we include the Z symmetry in the cost function of the system [73] such that The second term in the cost function penalizes the eigenstates whose Z eigenvalue is −1 and thus targets the desired ground state across the whole phase diagram.
Regarding the quantum simulator, we consider an array of qubits in a one-dimensional geometry with different levels of connectivity. The circuits in which one can perform two-qubit entangling gates, e.g. controllednot gate, between a pair of qubits with distance k (with k ≥ 1) are called ansatz with k-distance connectivity. For instance, 1-distance (2-distance) connectivity means that controlled-not gates are only allowed between nearest neighbor (and next nearest neighbor) qubits. As a two-qubit unitary operation between qubits m and n, we use the following operation: which can be implemented by the circuit in Fig. 3(a) [94].
As the figure shows, to implement this operation one needs three controlled-not gates. However, one can simplify the circuit to only use two controlled-not gates for the specific value of γ=1, as shown in Fig. 3(b). The two qubit unitary operations N m,n (ϕ) is used to generate a VQE ansatz inspired by quantum approximate optimization algorithm [39]. Depending on the connectivity between the qubits, one can apply N m,n (ϕ) between qubits of different distances. The VQE circuit with k-distance can be described by the following unitary operation.
where ϕ = (ϕ N −k ) are the tunable angles and R z (ϕ (0) 1 , · · · , ϕ (0) N ) is a set of rotation gates acting on each qubit as Therefore, each U k (ϕ) contains L k = 2N − k parameters to be optimized. This includes N − k two-qubit parameters ϕ In order to converge to the ground state, we need to concatenate several of these layers. The number of layers M is chosen such that the ground state is achieved with high accuracy. Quantum circuits for U k (with k > 2) can be similarly produced. In general, any single layer with a quantum circuit which implements U k contains 3(N − k) controlled-not gates. For the specific value of γ=1, using the simplified circuit implementation, each layer will instead contain 2(N − k) controlled-not gates. Depending on the connectivity of qubits different combinations of U k 's can be used as our VQE circuit. Note that, all the VQE simulations are performed using the simulator MindQuantum [95] with the choice of the classical optimizer being the L-BFGS algorithm [91]. Moreover, we apply the layer-recursive method [48] for simulating circuits with the number of layers M >1 which will start the training by only one layer then a new layer will be added after the existing circuit is optimized. The initial parameters for the first layer are randomly sampled from a normal distribution. When a new layer is added, the optimized parameters from the last layer of the preceding circuit serve as the initial parameters for the subsequent layer. This method has been proven to provide a great improvement for the convergence of the VQE optimization [48]. We first focus on 1-distance connectivity. A circuit of M layers is described by U M = M l=1 U 1 (ϕ l ), where ϕ l represents the tunable parameters at layer l. One can quantify the performance of the VQE through fidelity of the real ground state and the output of the trained quantum simulator F=|⟨GS|Ψ(ϕ)⟩| 2 . In Fig. 4(a)-(c), we plot the infidelity 1 − F for simulating the ground state of the Hamiltonian (3) with γ = 1 (i.e. transverse Ising model) in a system of length N = 10 across its phase diagram for M = 2, M = 4 and M = 6 layers, respectively. As the number of circuit layer increases, it is clear from the figures that the performance of the VQE improves significantly. In particular, it converges to a low infidelity in most of the phase diagram with the circuit depth of M = 6. However, it fails to yield the ground state with low infidelity when α is less than 0.5 where the Hamiltonian tends to be fully connected. The optimizer computes the iterations until the improvement of infidelity is less than a threshold, here ∼ 10 −9 . We choose this small threshold to make sure that the optimizer truly saturates. In Figs. 4(d)-(f), we plot the corresponding optimization iterations needed for convergence to the obtained infidelities. An interesting feature can be seen in Fig. 4(f), whose corresponding infidelity in Fig. 4(c) shows good performance in the entire phase diagram. In fact, the required iteration is larger along the phase transition line, in particular when α is small (i.e. the system is more long-ranged). This can be attributed to the complexity of the ground state at criticality which naturally exhibits more entanglement and long-range correlations.
One would naturally wonder if higher distance connectivity between the qubits could be of any benefits for VQE simulation of the Hamiltonian (3). To explore this, we consider a circuit with both 2-and 1-distance connectivity. Unlike the case of 1-distance connectivity, in this case the order of the layers is another degree of freedom which can be optimized to get better performance. To investigate different configurations, we keep M = 6 layers for the transverse Ising model (i.e. γ=1) of size N = 10 qubits. Later, we will consider other γ's too. We first consider alternating layers, implemented by U 1 and U 2 . Our simulation shows that the performance is obtained for the configuration For the sake of brevity, we only show the results for this optimal configuration in Fig. 5(a) and Fig. 5(c) for the achieved infidelity 1−F and the corresponding iterations, respectively. The infidelity shows improvement in comparison with the circuit using only U 1 with the same number of layers, see Fig. 4(c). For the case of 3-distance connectivity, our simulations show that the best performance can be obtained using The infidelity and its corresponding iteration is shown in Fig. 5(b) and Fig. 5(d). The overall improvement of infidelity in compare to the circuit with only 1-distance gates, see Fig. 4(c), and the circuit with 1-and 2-distance gates, see Fig. 5(a), is evident. In order to have a quantitative comparison between these circuit designs with different qubit distance connectivity, in Table. I we com- pare the minimum achievable fidelity F min , the average achievable fidelity F avg , and the maximum classical resources CR max across the entire phase diagram. Additionally, the number of two-qubit gates R Q is included as a metric for assessing quantum resources. Given that F≈1 can be readily achieved with smaller θ values in the phase diagram, the maximum achievable fidelity F max is omitted in this discussion. Indeed, the information shown in Table. I provides a comparison over the worst and the average performance of different circuit designs in the entire phase diagram. Clearly, the circuit represented by U 1 U 1 U 2 U 2 U 3 U 3 achieves higher fidelity even at its worst point of the phase diagram. Interestingly, the number of controlled-not gates R Q is also minimum for this circuit configuration, showing the minimal demand of quantum resources. Similarly, the maximum classical resources CR max in the entire phase diagram is far less for this circuit compared to the others.
If we restrict ourselves to 1-and 2-distance gates, we also see that the overall performance of U 1 U 1 U 1 U 2 U 2 U 2 outperforms the other configurations such as U 1 U 2 U 1 U 2 U 1 U 2 . This can be seen in Table. I. While the minimum fidelities F min are pretty close to each other, the classical resources CR max for the former circuit is much less. This shows that by grouping the layers with long-distance connectivity at the beginning of the circuit, one can achieve better overall performance. Indeed, such circuit design is capable for creating long-distance correlations between distant qubits using fewer layers. This is particularly useful for long-range interacting systems  Table I. The minimum fidelity Fmin, and the average fidelity Favg achievable across the entire phase diagram for various circuit designs with 6 layers are presented, along with the maximum required classical resources, CRmax. Additionally, the number of two-qubit gates RQ is provided for each specific circuit design. Each circuits i1i2 · · · i k represents a quantum circuit with unitary operation of the type Ui k · · · Ui 2 Ui 1 . Figure 6. The required circuit layer M to obtain fidelity F≥0.95 for simulating the ground state of the transverse Ising Hamiltonian (i.e. γ=1) with α=0.5 and θ=0.425π, in various system sizes N as well as different maximum qubit connectivities k.
which inherently have larger values of long-range correlations. As a result, this specific design of the circuit achieves better performance for simulating the ground state of long-range interacting systems.

V. SCALING UP
In this section, we focus on the impact of qubit connectivity for simulation of larger system sizes. In particular, we consider the optimal circuit design in which the layers with the same qubit connectivity are grouped together and positioned in descending order. For instance, in a circuit with M = 5 layers and the maximum qubit connectivity k=3, we use U 1 U 2 U 2 U 3 U 3 . For a transverse Ising chain (i.e. γ=1) with α=0.5, which represents a strong long-range interaction, and θ=0.425π, which is close to the critical point, we simulate the ground state of the Hamiltonian to obtain fidelity F ≥ 0.95. This choice of α and θ represents the most difficult part of the phase diagram for the VQE convergence. By fixing qubit connectivity k, we determine how many circuit layers M is required to achieve the fidelity F ≥ 0.95. The results for various system sizes are shown in Fig. 6. As one can see, for a given qubit connectivity k, at least k layers are required. This guarantees that layers with all available connectivities, namely k, k − 1, · · · , 1, are required for achieving high fidelity VQE performance. While large k does not show any advantage for small system sizes, its benefit becomes very evident for large systems. In fact, when k is larger, the VQE procedure converges with fewer layers as the system size increases. This shows that having larger qubit connectivity is hugely beneficial for scaling up the VQE simulation of long-range interacting systems.

VI. GENERALITY WITH RESPECT TO γ
So far, we have focused on transverse Ising model with γ=1. In this section, we show that our strategy for circuit design is applicable for other values of γ. We have to note that for γ̸ =1 the circuit design for the two-qubit gate operation N m,n (ϕ) has to change as shown in Fig. 3(a). For a chain of size N =10 and a circuit with M =6 layers, we present the minimum achievable fidelity F min , average achievable fidelity F avg , and maximum classical resources CR max , obtained across the entire phase diagram, in Table. II for three different values of γ. As F≈1 can be readily attained with smaller θ values in the phase diagram, we have excluded the maximum achievable fidelity F max from this discussion. The superiority of the circuit U 1 U 1 U 2 U 2 U 3 U 3 is evident in achieving the best fidelity F min , as the worst outcome of the circuit in the entire phase diagram. Remarkably, this circuit also demands less classical resources in comparison with other circuit designs. This shows that grouping the layers with long-distance connectivity at the beginning of the circuit can significantly enhance the performance of the VQE in terms of both fidelity and required classical resources.

VII. SPIN SQUEEZING
Generating the eigenstates of a Hamiltonian is not the only application of variational quantum algorithms [31]. In principle, any problem that can be written variationally in terms of minimization of a cost function can be solved on a quantum simulator. This has led to a wide range of variational quantum algorithms for solving problems in linear algebra [96,97], combinatorial optimization [39,41] and machine learning [33][34][35][36][37][38]. In this section, we show that our circuits with k-distance qubit connectivity can indeed provide advantage for producing an important class of states with a non-classical feature called spin squeezing.
Spin squeezed quantum states [84,85] are an important class of entangled states which have fundamental and practical applications. Fundamentally, these states can be used for testing the foundations of quantum mechanics through violating Bell inequalities [86,87]. From a practical perspective, these states can be regarded as resource for quantum metrology protocols which can achieve quantum enhanced sensitivity [88]. Spin squeezing is defined for a set of spin particles whose collective spin has a net polarization in one direction while the fluctuations of the corresponding transverse spin component is suppressed. We define the collective spin components in x, y and z directions as S x,y,z tot = 1/2 k σ x,y,z k . Assuming that the net spin polarization is in the z direction one can quantify spin squeezing through r = max −10 log 10 ξ 2 , 0 in which the parameter ξ 2 is defined as [98] where, Var[•] stands for variance and angle β spans the collective spin in the transverse direction. To have nonzero squeezing r, one has to reach ξ 2 < 1, which requires strong spin polarization |⟨S z tot ⟩| 2 along the z direction and suppressed variance for the collective spin component in the transverse direction.
In practice, creating spin-squeezed states is very challenging. Here, we consider three different methods for generating such states. The first approach is to search for spin squeezing in the ground state of long-range Ising Hamiltonian in the transverse field. The second approach is to use the quench dynamics of such Hamiltonian for creating spin squeezed states. Finally, the third approach is to use a variational quantum circuit to generate such states using the parameter ξ 2 , in Eq. (9), as the cost function [54,90]. Let's consider N spin-1/2 particles interacting via transverse Ising Hamiltonian (3) with γ=1. For any given values of α and θ one can compute the squeezing parameter r for the ground state |GS⟩. In Fig. 7(a), we plot the spin squeezing r as a function of θ in the whole phase diagram, for various values of α in a system of N = 10 particles. Two main features can be observed: (i) decreasing α (i.e. making the interaction more long-ranged) monotonically enhances the spin squeezing; (ii) for a fixed α the spin squeezing takes its maximum around the phase transition point. It is worth emphasizing that the minimum of ξ 2 , in Eq. (9), takes place at β = 0 in the whole phase diagram for the ground state squeezing.
Although, the ground state of the system reveals squeezing, its value can be enhanced through quench dynamics [99]. In this case, we initialize the system in the state |Ψ(0)⟩ = |0, 0, · · · , 0⟩. In general, this state is not an eigenvector of the transverse Ising Hamiltonian (3) and thus evolves to |Ψ(t)⟩ = e −iHt |Ψ(0)⟩. By computing the squeezing parameter r with respect to quantum state |Ψ(t)⟩ one can explore the dynamics of spin squeezing as a function of time. In Figs. 7(b)-(c) we plot the squeezing parameter r(t) as a function of time for various choices of α when θ is fixed to θ = π/8 and θ = π/4, respectively. As the figures show r starts from 0 and increases by time to reach its maximum at an optimal time t = t * and then decreases again. The maximum spin squeezing r(t * ) quantifies the capacity of the Hamiltonian for dynamical production of spin squeezing. In Fig. 7(d), we plot the maximum obtainable squeezing r(t * ) as a function of θ and α for a system of length N = 10. Unlike the ground state, the optimal β which minimizes ξ 2 , in Eq. (9), varies in the dynamical case as θ and α change. It should be noted that using the dynamical approach results in more squeezing than the ground state of the same Hamiltonian, literally throughout the whole phase diagram.
As an alternative to the ground state or the dynamics of the transverse Ising Hamiltonian (3) with γ=1 for generating spin squeezed states, one can directly apply a variational quantum algorithm with ξ 2 , given in Eq. (9), as the cost function. Without loss of generality we fix β = 0. We use the same quantum circuits which we have used for simulating the ground state. The obtainable squeezing r from a circuit of M =6 layers with different designs and the initial state |0, 0, · · · , 0⟩ is shown in Ta-ble. III. For each circuit, the obtainable squeezing r is averaged over 100 repetition of the protocol and the uncertainty in the values of Table. III is determined through the standard deviation. Interestingly, the circuits which perform better for simulating the ground state also provide higher squeezing rate r. In addition, variational quantum approach for generating spin squeezed states results in higher squeezing than both the ground state and the dynamics of long-range interacting Hamiltonian across the whole phase diagram. For instance, the highest squeezing rate r which one can obtain by the ground state and the dynamics in the whole phase diagram is r = 4.4 and r = 5.1, respectively. Indeed, as Table. III shows, one can outperform both of these strategies using variational quantum algorithms.

VIII. CONCLUSION
The most promising approach to achieve quantum advantage with imperfect NISQ simulators is variational quantum algorithm. VQE, as one the most widely used VQAs, has been developed for generating the ground state of many-body systems on current quantum simulators. In this paper, we have focused on digital VQE simulation of a large class of long-range Hamiltonians with a quantum phase transition. We find out that by making the interaction more long-ranged, the VQE algorithm results in lower fidelities and demands more optimization iterations. The situation becomes worse as the Hamiltonian gets closer to its criticality. Interestingly, our simulations show that qubit connectivity in quantum simulators has a direct impact on the performance of VQE. With the possibility of applying two-qubit entangling gates between distant qubits the achievable fidelity is improved, in particular, for Hamiltonians which are more long-ranged. Interestingly, such enhancement in fidelity is combined with improvement in demanding both quantum and classical resources. The performance is further enhanced if the layers with larger distance connectivity grouped together and act before the layers with shorter distance connectivity.
As an alternative application, we have shown that the same design of circuits for variationally generating spin squeezed states, as resource for quantum metrology. The variational generation of spin squeezing outperforms the results from the ground state and quench dynamics in long-range interacting systems.
Current NISQ simulators provide various forms of qubit connectivity. In certain physical setups, such as superconducting devices, the qubit connectivity is determined by fabrication and is usually limited to nearest neighbors. In other systems, such as ion-traps and Rydberg atoms, the connectivity is more versatile and operations between distant qubits are also possible. This suggests that ion-traps and Rydberg atoms are more suitable for simulating long-range interacting systems.

IX. DATA AVAILABILITY
All the codes for training the circuits are available online [100], and the data presented on this paper will be provided upon reasonable requests from the authors.