{'fraction_non_alphanumeric': 0.03508023877819986, 'fraction_numerical': 0.011551283045197301, 'mean_word_length': 4.257591196250255, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical study of the hard-core Bose-Hubbard Model on an Infinite Square Lattice

We present a study of the hard-core Bose-Hubbard model at zero temperature on an infinite square lattice using the infinite Projected Entangled Pair State algorithm [Jordan et al., Phys. Rev. Lett. 101, 250602 (2008)]. Throughout the whole phase diagram our values for the ground state energy, particle density and condensate fraction accurately reproduce those previously obtained by other methods. We also explore ground state entanglement, compute two-point correlators and conduct a fidelity-based analysis of the phase diagram. Furthermore, for illustrative purposes we simulate the response of the system when a perturbation is suddenly added to the Hamiltonian.


I. INTRODUCTION
The physics of interacting bosons at low temperature has since long attracted considerable interest due to the occurrence of Bose-Einstein condensation 1 . The Bose-Hubbard model, a simplified microscopic description of an interacting boson gas in a lattice potential, is commonly used to study related phenomena, such as the superfluid-to-insulator transitions in liquid helium 2 or the onset of superconductivity in granular superconductors 3 and arrays of Josephson junctions 4 . In more recent years, the Bose-Hubbard model is also employed to describe experiments with cold atoms trapped in optical lattices 5 .
As in most many-body systems, the theoretical study of interacting bosons cannot rely only on the few exact solutions available. Numerical results are also needed, but these are not always easy to obtain. Indeed, the exponential growth of the Hilbert space dimension in the lattice size (even after placing a bound on the number of bosons allowed on each of its sites) implies that exact diagonalization techniques are only capable of addressing very small lattices. Thus, in order to study the ground state properties of the Bose-Hubbard model on e.g. the square lattice, as is the goal of the present work, a number of more elaborate techniques, such as mean field theory, spin-wave calculations or quantum Monte Carlo are traditionally used (see e.g. 6 and references therein).
Recently, a new class of simulation algorithms for two-dimensional systems, based on tensor networks, has gained much momentum. The basic idea is to use a network of tensors to efficiently represent the state of the lattice. Specifically, the so-called tensor product states 8,9 (TPS) or projected entangled-pair states 10,11 (PEPS) are used to (approximately) express the d N coefficients of the wave function |Ψ of a lattice of N sites in terms of just N tensors, in such a way that only O(N ) coefficients are actually specified. After optimizing these tensors so that |Ψ represents e.g. the ground state of the system, one can then extract from them a number of properties, including the expected value of arbitrary local observables. Moreover, in systems that are invariant under translations, the tensor network is made of copies of a small number of tensors. This leads to an even more compact description that depends on just O(1) parameters. The later is the basis of the infinite PEPS (iPEPS) algorithm 12 , which addresses infinite lattices and can thus be used to compute thermodynamic properties directly, without need to resort to finite size scaling techniques.
In this work we initiate the exploration of interacting bosons in an infinite 2D lattice with tensor network algorithms. We use the iPEPS algorithm 12,13 to characterize the ground state of the hard-core Bose-Hubbard (HCBH) model, namely the Bose Hubbard model in the hard-core limit, where either zero or one bosons are allowed on each lattice site. Although no analytical solution is known for the 2D HCBH model, there is already a wealth of numerical results based on mean-field theory, spin-wave corrections and stochastic series expansion 6 . These techniques have been quite successful in determining some of the properties of the ground state of the 2D HCHB model, such as its energy, particle density or condensate fraction. Our goal in this paper is twofold. Firstly, by comparing our results against those of Ref. 6 , we aim to benchmark the performance of the iPEPS algorithm in the HCBH model. Secondly, once the validity of the iPEPS algorithm for this model has been established, we use it to obtain results that are harder to compute with (or simply well beyond the reach of) the other approaches. These include the analysis of entanglement, two-point correlators, fidelities between different ground states 14,15,16 , and the simulation of time evolution.
We note that the present results naturally complement those of Ref. 11 for finite systems, where the PEPS algorithm 10 was used to study the HCBH model in a lattice made of at most 11 × 11 sites.
The rest of the paper is organized as follows. Sect. II introduces the HCBH model and briefly reviews the iPEPS algorithm. Sect. III contains our numerical results for the ground state of the 2D HCBH model. These include the computation of local observables such as the energy per lattice site, the particle density and the condensate fraction. We also analyze entanglement, twopoint correlators and ground state fidelities. Finally, the simulation of time evolution is also considered. Sect. IV contains some conclusions.

II. MODEL AND METHOD
In this section we provide some basic background on the HCBH model, as well as on the iPEPS algorithm.

A. The Hard Core Bose-Hubbard Model
The Bose-Hubbard model 2 with on-site and nearest neighbour repulsion is described by the Hamiltonian where a † i , a i are the usual bosonic creation and annihilation operators,n i =ρ i ≡ a † i a i is the number (density) operator at site i, J is the hopping strength, µ is the chemical potential, and V 1 , V 2 ≥ 0. The four terms in the above equation describe, respectively, the hopping of bosonic particles between adjacent sites (J), a single-site chemical potential (µ), an on-site repulsive interaction (V 1 ) and an adjacent site repulsive interaction (V 2 ).
Here we shall restrict our attention to on-site repulsion only (V 2 = 0) and to the so-called hard-core limit in which this on-site repulsion dominates (V 1 → ∞). Under these conditions the local Hilbert space at every site describes the presence or absence of a single boson and has dimension 2. With the hard-core constraint in place, the Hamiltonian becomes where a † i , a i are now hard-core bosonic operators obeying the commutation relation, A few well-known facts of the HCBH model are: and it thus has a U (1) symmetry, corresponding to transforming each site l by e iφn l , φ ∈ [0, 2π).
(ii) Duality transformation.-In addition, the transformation a l → a † l applied on all sites l of the lattice maps H HC (µ) into H HC (−µ) (up to an irrelevant additive constant). Accordingly, the model is self-dual at µ = 0, and results for, say, µ > 0 can be easily obtained from those for µ < 0.
(iii) Equivalence with a spin model.-The HCBH model is equivalent to a quantum spin 1 2 model, namely the ferromagnetic quantum XX model, which is obtained from H HC with the replacements where σ x , σ y and σ z are the spin 1 2 Pauli matrices. In particular, all the results of this paper also apply, after a proper translation, to the ferromagnetic quantum XX model on an infinite square lattice.
(iv) Ground-state phase diagram.-The hopping term in H HC favors delocalization of individual bosons in the ground state, whereas the chemical potential term determines the ground state bosonic density ρ, For µ negative, a sufficiently large value of |µ| forces the lattice to be completely empty, ρ = 0. Similarly, a large value of (positive) µ forces the lattice to be completely full, ρ = 1, as expected from the duality of the model. In both cases there is a gap in the energy spectrum and the system represents a Mott insulator. When, instead, the kinetic term dominates, the density has some intermediate value 0 < ρ < 1, the cost of adding/removing bosons to the system vanishes, and the system is in a superfluid phase 2 . The latter is characterized by a finite fraction of bosons in the lowest momentum modẽ a k=0 ≡ (1/N ) i a i , that is by a non-vanishing condensate fraction ρ 0 , In the thermodynamic limit, N → ∞, a non-vanishing condensate fraction is only possible in the presence of off-diagonal long range order (ODLRO) 17 , or a † j a i = 0 in the limit of large distances |i − j|, given that (v) Quantum phase transition.-Between the Mott insulator and superfluid phases, there is a continuous quantum phase transition 2 , tuned by µ J .

B. The algorithm
The state |Ψ of the infinite square lattice is represented using a TPS 8 or PEPS 10 that consists of just two different tensors A and B that are repeated in a The coefficients of tensors A and B are determined with the iPEPS algorithm 12 . Specifically, the ground state |Ψ GS of the HCBH model is obtained by simulating an evolution in imaginary time according to H HC , exploiting that We have also used the iPEPS algorithm to simulate (real) time evolution starting from the ground state |Ψ GS and according to a modified Hamiltonian H (see Eq. 12), These simulations, as well as the computation of expected values of local observables from the resulting state, involve contracting an infinite 2D tensor network. This is achieved with techniques developed for infinite 1D lattice systems 7 , namely by evolving a matrix product state (MPS). An important parameter in these manipulations is the bond dimension χ of the MPS, which parameterizes how many correlations the latter can account for. We refer to 12 for a detailed explanation of the iPEPS algorithm. In what follows we briefly comment on the main sources of errors and on the simulation costs.
We distinguish three main sources of errors in the simulations, one due to structural limitations in the underlying TPS/PEPS ansatz and two that originate in the particular way the iPEPS algorithm operates: limits the amount of correlations the TPS/PEPS can carry. A typical state of interest |Ψ , e.g. the ground state of a local Hamiltonian, requires in general a very large bond dimension D if it is to be represented exactly. However, a smaller value of D, say D ≥ D Ψ for some value D Ψ that depends on |Ψ , often already leads to a good approximate representation, in that the expected values of local observables are reproduced accurately. However, if D < D 0 , then the numerical estimates may differ significantly from the exact values, indicating that the TPS/PEPS is not capable of accounting for all the correlations/entanglement in the target state |Ψ .
(ii) MPS bond dimension χ.-Similarly, using a finite MPS bond dimension χ implies that the contraction of the infinite 2D tensor network (required both in the simulation of real/imaginary time evolution and to compute expected values of local observables) is only approximate. This may introduce errors in the evolved state, or in the expected value of local observables even when the TPS/PEPS was an accurate representation of the intended state.
(iii) Time step.-A time evolution (both in real or imaginary time) is simulated by using a Suzuki-Trotter expansion of the evolution operator (e −itH or e −τ H ), which involves a time step (δt or δτ ). This time step introduces an error in the evolution that scales as some power of the time step. Therefore this error can be reduced by simply diminishing the time step.
The cost of the simulations scales as O(χ 3 D 6 + χ 2 D 8 d) (here we indicate only the leading orders in χ and D; the cost of the simulation is also roughly proportional to the inverse of the time step). This scaling implies that only small values of the bond dimensions D and χ can be used in practice. In our simulations, given a value of D (D = 2, 3 or 4), we choose a sufficiently large χ (in the range 10 − 40) and sufficiently small time step (δt or δτ ) such that the results no longer depend significantly on these two parameters. In this way the bond dimension D is the only parameter on which the accuracy of our results depends.
On a 2.4 GHz dual core desktop with 4 Gb of RAM, computing a superfluid ground state (e.g. µ = 0) with D = 2, χ = 20 and with δτ decreasing from 10 −1 to 10 −4 requires about 12 hours. Computing the same ground state with D = 3 and χ = 40 takes of the order of two weeks.

III. RESULTS
In this section we present the numerical results obtained with the iPEPS algorithm.
Without loss of generality, we fix the hopping strength J = 1 and compute an approximation to the ground state |Ψ GS of H HC for different values of the chemical potential µ. Then we use the resulting TPS/PEPS to extract the expected value of local observables, analyze ground state entanglement, compute two-point correlators and fidelities, or as the starting point for an evolution in real time.
In most cases we only report results for µ ≤ 0 (equivalently, density 0 ≤ ρ ≤ 0.5) since due to the duality of the model, results for positive µ (equivalently, 0.5 ≤ ρ ≤ 1) can be obtained from those for negative µ.

A. Local observables and phase diagram
Particle density ρ. -Fig. 2 shows the density ρ as a function of the chemical potential µ in the interval −4 ≤ µ ≤ 0. Notice that ρ = 0 for µ ≤ −4, since each single site is vacant. Our results are in remarkable agreement with those obtained in Ref. 6 with stochastic series expansions (SSE) for a finite lattice made of 32 × 32 and with a mean field calculations plus spin wave corrections (SW). We note that the curves ρ(µ) for D = 2 and D = 3 are very similar.
Energy per site ǫ.
- Fig. 2 also shows the energy per site ǫ as a function of the density ρ. This is obtained by computing ǫ(µ) and then replacing the dependence on µ with ρ by inverting the curve ρ(µ) discussed above. Again, our results for ǫ(ρ) are in remarkable agreement with those obtained in Ref. 6 with stochastic series expansions (SSE) for a finite lattice made of 32 × 32. They are also very similar to the results coming from mean field calculations with spin wave corrections (SW) of Ref 6 , and for small densities reproduce the scaling (valid only in the regime of a very dilute gas) predicted in Ref. 18 by using field theory methods based on a summation of ladder diagrams. Once more, the curves ǫ(ρ) obtained with bond dimension D = 2 and D = 3 are very similar, although D = 3 produces slightly lower energies.
Condensate fraction ρ 0 .-In order to compute the condensate fraction ρ 0 , we exploit that the iPEPS algorithm induces a spontaneous symmetry breaking of particle number conservation. Indeed, one of the effects of having a finite bond dimension D is that the TPS/PEPS that minimizes the energy does not have a well-defined particle number. As a result, instead of having a i = 0, we obtain a non-vanishing value a i = 0 such that In other words, the ODLRO associated with the presence of superfluidity, or a finite condensate fraction, can be computed by analysing the expected value of a l , where the phase ϕ is constant over the whole system but is otherwise arbitrary. The condensate fraction ρ 0 shows that the model is in an insulating phase for |µ| ≥ 4 (ρ = 0, 1) and in a superfluid phase for −4 < µ < 4 (0 < ρ < 1), with a continuous quantum phase transition occurring at |µ| = −4, as expected. However, this time the curves ρ 0 (ρ) obtained with D = 2 and D = 3 are noticeably different, with D = 3 results again in remarkable agreement with the SSE and SW results of Ref. 6 .

B. Entanglement
The iPEPS algorithm is based on assuming that a TPS/PEPS offers a good description of the state |Ψ of the system. Results for small D will only be reliable if |Ψ has at most a moderate amount of entanglement. Thus, in order to understand in which regime the iPEPS algorithm should be expected to provide reliable results, it is worth studying how entangled the ground state |Ψ GS is as a function of µ.
The entanglement between one site and the rest of the lattice can be measured by the degree of purity of the reduced density matrix ̺ 1 for that site, as given by the norm r of the Bloch vector r. If the lattice is in a product or unentangled state, then each site is in a pure state, corresponding to purity r = 1. On the other hand, if the lattice is in an entangled state, then the onesite reduced density matrix will be mixed, corresponding to purity r < 1. Accordingly, one can think of r as measuring the amount of entanglement between one site and the rest of the lattice, with less purity corresponding to more entanglement. Fig. 3 shows the purity r as a function of the chemical potential. In the insulating phase (µ ≤ −4), the ground state of the system consists of a vacancy on each site. In other words, it is a product state, r = 1. Instead, For µ > −4 the ground state is entangled. Several comments are in order: (i) The purity r(µ) for D = 3 is smaller than that for D = 2 by up to 3%. This is compatible with the fact that the a TPS/PEPS with larger bond dimension D can carry more entanglement.
(ii) Results for D = 2, 3 seem to indicate that the ground state is more entangled (r is smaller) deep into the superfluid phase (e.g. µ = 0) than at the continuous quantum phase transition µ = −4. This is in sharp contrast with the results obtained e.g. for the 2D quantum Ising model 12 , where the quantum phase transition displays the most entangled ground state. However, notice that in the Ising model the system is only critical at the phase transition whereas in the present case criticality extends throughout the superfluid phase. Each value of µ in the superfluid phase corresponds to a fixed point of the RG flow. That is, in moving away from the phase transition we are not following an RG flow. Therefore, the notion that entanglement should decrease along an RG flow 20 , as observed in the 2D Ising model, is not applicable for the HCBH model.
(iii) Accordingly, we expect that the iPEPS results for small D become less accurate as we go deeper into the superfluid phase (that is, as we approach ρ = 0.5). This is precisely what we observe: the curves ρ 0 (ρ) for D = 2 and D = 3 in Fig. 2 differ most at ρ = 0.5. Fig. 3 also shows the entanglement entropy for the reduced density matrix ̺ L (L = 1, 2, 4) corresponding to one site, two contiguous sites and a block of 2 × 2 sites respectively. The entanglement entropy vanishes for an unentangled state and is non-zero for an entangled state. The curves S(̺ L ) confirm that the ground state of the HCBH model is more entangled deep in the superfluid phase than at the quantum phase transition point.

C. Correlations
From a TPS/PEPS for the ground state |Ψ GS it is easy to extract equal-time two-point correlators. For illustrative purposes, Fig. (4) shows a connected two-point correlation function C(s), between two sites that separated s lattice sites along the horizontal directionx. The plot corresponds to a superfluid ground state, µ = 0, where C(s) displays an exponential decay, in spite of the fact that the Hamiltonian is gapless.
The results show that while for short distances s = 0, 1, 2 the correlator C(s) is already well converged with respect to D, for larger distances s the correlator still depends significantly on D. This seems to indicate that while the iPEPS algorithm provides remarkably good results for local observables already for affordably small values of D, a larger D might be required in order to also obtain accurate estimates for distant correlators.

D. Fidelity
Given two ground states |Ψ GS (µ 1 ) and |Ψ GS (µ 2 ) , corresponding to different chemical potential µ, the fidelity per site f 15 , defined through can be used as a means to distinguish between qualitatively different ground states 14,15 . In the above expression, N is the number of lattice sites and the thermodynamic limit N → ∞ is taken. Importantly, the fidelity per site f (µ 1 , µ 2 ) remains finite in this limit, even though the overall fidelity | Ψ GS (µ 1 )|Ψ GS (µ 2 ) | vanishes. In a sense, f (µ 1 , µ 2 ) captures how quickly the overall fidelity vanishes. Fortunately, the fidelity per site f (µ 1 , µ 2 ) can be easily computed within the framework of the iPEPS algorithm 16 . In the present case, before computing the overlap each ground state is rotated according to e iϕσz /2 , where ϕ is the random condensate phase of Eq. 8. In this way all the ground states have the same phase ϕ = 0. The fidelity per site f (µ 1 , µ 2 ) is presented in Fig. 5. The plateau-like behavior of f (µ 1 , µ 2 ) for points within the separable Mott-Insulator phase (µ 1 , µ 2 ≤ −4 or µ 1 , µ 2 ≥ 4) is markedly different from that between ground states in the superfluid region (−4 ≤ µ 1 , µ 2 ≤ 4), where the properties of the system vary continuously. Moreover, similarly to what has been observed for the 2D quantum Ising model 16 or in the 2D quantum XYX model 19 , the presence of a continuous quantum phase transition between insulating and superfluid phases in the 2D HCBH model is signaled by pinch points of f (µ 1 , µ 2 ) at µ 1 = µ 2 = ±4. That is, the qualitative change in ground state properties across the critical point is evidenced by a rapid, continuous change in the fidelity per lattice site as one considers two ground states on opposite sides of the critical point and moves away from it.

E. Time evolution
An attractive feature of the algorithms based on tensor networks is the possibility to simulate (real) time evolution. A first example of such simulations with the iPEPS algorithm was provided in Ref. 13 , where an adiabatic evolution across the quantum phase transition of the 2D quantum compass orbital model was simulated in order to show that the transition is of first order.
The main difficulty in simulating a (real) time evolution is that, even when the initial state |Ψ(0) is not very entangled and therefore can be properly represented with a TPS/PEPS with small bond dimension D, entanglement in the evolved state |Ψ(t) will typically grow with time t and a small D will quickly become insufficient. Incrementing D results in a huge increment in computational costs, which means that only those rare evolutions where no much entanglement is created can be simulated in practice.
For demonstrative purposes, here we have simulated the response of the ground state |Ψ GS of the HCBH model at half filling (ρ = 0.5 or µ = 0) when the Hamiltonian H HC is suddenly replaced with a new Hamiltonian H given by where γ = 0.2 and, importantly, the perturbation V respects translation invariance. As the starting point of the simulation, we consider a TPS/PEPS representation of the ground state with bond dimension D = 2, obtained as before through imaginary time evolution. Fig. 6 shows the evolution in time of the expected value per site of the energies H HC and H , as well as the density ρ and condensate fraction ρ 0 . Notice that the expected value of H should remain constant through the evolution. The fluctuations observed in H , of the order of 0.2% of its total value, are likely to be due to the small bond dimension D = 2 and indicate the scale of the error in the evolution. The simulation shows that, as a result of having introduced a perturbation V that does not preserve particle number, the particle density ρ oscillates in time. The condensate fraction, as measured by | a l | 2 , is seen to oscillate twice as fast.

IV. CONCLUSION
In this paper we have initiated the study of interacting bosons on an infinite 2D lattice using the iPEPS algorithm. We have computed the ground state of the HCBH model on the square lattice as a function of the chemical potential. Then we have studied a number of properties, including properties that can be easily accessed with other techniques 6 , as is the case of the expected value of local observables, as well as properties whose computation is harder, or even not possible, with previ-ous techniques.
Specifically, using a small bond dimension D = 2, 3 we have been able to accurately reproduce the result of previous computations using SSE and SW of Ref. 6 for the expected value of the particle density ρ, energy per particle ǫ and condensate fraction ρ 0 , throughout the whole phase diagram of the model, which includes both a Mott insulating phase and a superfluid phase, as well as a continuous phase transition between them. Interestingly, in the superfluid phase the TPS/PEPS representation spontaneously breaks particle number conservation, and the condensate fraction can be computed from the expected value of the annihilation operator, ρ 0 = | a l | 2 .
We have also conducted an analysis of entanglement, which revealed that the most entangled ground state corresponds to half filling, ρ = 0.5. This is deep into the superfluid phase and not near the phase transition, as in the case of the 2D quantum Ising model 12 . Furthermore, inspection of a two-point correlator at half filling showed much faster convergence in the bond dimension D for short distances than for large distances. Also, pinch points in plot of the fidelity f (µ 1 , µ 2 ) were consistent with continuous quantum phase transitions at µ = ±4.
Finally, we have also simulated the evolution of the system, initially in the ground state of the HCHB model at half filling, when a translation invariant perturbation is suddenly added to the Hamiltonian. Now that the validity of the iPEPS algorithm for the HCBH model (equivalently, the quantum XX spin model) has been established, there are many directions in which the present work can be extended. For instance, one can easily include nearest neighbor repulsion, V 2 = 0, (corresponding to the quantum XXZ spin model) and/or investigate a softer-core version of the Bose Hubbard model by allowing up to two or three particles per site.