Estimating truncation effects of quantum bosonic systems using sampling algorithms

To simulate bosons on a qubit- or qudit-based quantum computer, one has to regularize the theory by truncating infinite-dimensional local Hilbert spaces to finite dimensions. In the search for practical quantum applications, it is important to know how big the truncation errors can be. In general, it is not easy to estimate errors unless we have a good quantum computer. In this paper, we show that traditional sampling methods on classical devices, specifically Markov Chain Monte Carlo, can address this issue for a rather generic class of bosonic systems with a reasonable amount of computational resources available today. As a demonstration, we apply this idea to the scalar field theory on a two-dimensional lattice, with a size that goes beyond what is achievable using exact diagonalization methods. This method can be used to estimate the resources needed for realistic quantum simulations of bosonic theories, and also, to check the validity of the results of the corresponding quantum simulations.


Introduction
The Hilbert space for bosons is infinite-dimensional.To simulate bosons on a qubit-or qudit-based quantum computer, one has to introduce a finite-dimensional approximation of the theory by truncating the Hilbert space [1,2,3,4,5,6,7,8].Sometimes this truncation is referred to as digitization and it is one of the necessary steps in the construction of efficient quantum algorithms to simulate quantum gauge theories [9,10,11,12,13], and to estimate quantum resources [14,15,16].We will use truncation and digitization to mean the same thing in the rest of the manuscript, and this should not be confused with the discretization of the spacetime continuum on the lattice.Some errors associated with digitization are inevitable, and one needs to know how big they can be.While previous studies on truncated bosonic Hilbert spaces [1,2,3,4,5,6,7,8] have looked at the scaling of the errors with the digitization spacing (or equivalently, the number of qubits), this does not seem to be a straightforward task in general.In fact, when the bosonic systems investigated are too large, one would need to run the digitized theory on a good (noiseless) quantum computer and also efficiently compute the same observable with classical algorithms in order to benchmark the result.So far, quantitative assessment of the digitization errors using classical devices is limited to rather small systems because the only known approaches applicable to generic theories are based on the explicit construction of Hilbert space whose dimensionality increases exponentially with respect to the system size.A thorough numerical study of digitization errors for a ϕ 4 lattice scalar theory was done in Ref. [6], where the authors studied a single-site (one bosonic degree of freedom) and a twosite (two bosonic degrees of freedom) system.However, it was difficult to go to larger lattice sizes, due to the exponential growth of memory resources with the number of lattice sites (bosonic degrees of freedom).In this paper, we will resolve this issue for a rather generic class of bosonic theories and reduce the required resources from exponential to polynomial when estimating expectation values.We point out that this is the only known method able to do so for a rather generic class of bosonic systems in arbitrary dimensions.As a modest demonstration, we will study the digitization effect on a two-dimensional 16-site lattice model (16 bosonic degrees of freedom).
It is widely believed that the truncation effect vanishes exponentially, e.g., the low-lying energy eigenvalues have correction suppressed exponentially with respect to the truncation level [17,4,3,5,18,19,20,21].Although this scaling has been verified for concrete examples, there is no proof applicable to a wide class of theories.Furthermore, there is an implicit assumption, i.e., the wave function decays exponentially fast as |x| becomes large.Note also that, even if the exponential accuracy of the digitization is valid, the precise value of the error for each theory is not immediately clear.However, it is practically impossible to study the digitization effects via the explicit construction of the Hilbert space except for small systems, e.g.lattice field theories with a few local degrees of freedom and less than a dozen lattice sites [5,6].It is a frustrating situation that prevents us from quantitative resource estimates for quantum simulations.A related issue is that it is not easy to check the validity of the results of the quantum simulation with specific digitization on NISQ devices: the presence of noise can quickly invalidate the results without a working errorcorrection strategy.It would be practically useful if we can do some calculations on classical devices to establish robust truncation error quantification which, as a byproduct, offers us some ways to cross-check quantum simulations on NISQ devices. 2n this paper, we show that a Markov Chain Monte Carlo (MCMC) method can be used to estimate the digitization effects in a rather generic class of bosonic systems.MCMC methods [22] (see Ref. [23] for an elementary introduction) are straightforward in the original bosonic theory before digitization.Indeed, one can use the Euclidean path-integral method combined with MCMC to study some features of the quantum systems as long as there is no sign problem.If similar simulations are doable for the digitized theory, one can estimate the digitization effects.
A caveat is that the absence of the sign problem in the Euclidean path integral without digitization does not necessarily guarantee the absence of the sign problem in the digitized theory.In this work, we will show that we can apply MCMC to a wide class of theories without having the sign problem, at least for a certain digitization scheme.As a result, we are able to properly assess the amount of digitization effects even for large systems, well beyond the limits of exact diagonalization techniques.
We consider the generic system of N bos bosons consisting of coordinate variables xi (i = 1, 2, • • • , N bos ) and conjugate momentum variables pi that satisfy the canonical com-mutation relation3 [x j , pk ] = iδ jk . (1) We assume that the Hamiltonian is given by Ĥ = 1 2 where ) is a real function bounded from below and represents the potential energy as a function of the coordinate variables only.This is a rather generic class and it even includes SU(N ) Yang-Mills theories described by using the orbifold lattice formulation [24,25].We use the truncation in the coordinate basis defined in Sec. 2. As we will see, this scheme admits the MCMC simulations without sign problem.This paper is organized as follows.In Sec. 2, we introduce the digitization scheme associated with the coordinate basis.Firstly, the case of a single bosonic variable is explained, and then it is generalized to the case of a generic number of variables.The digitization of (2 + 1)-dimensional scalar Quantum Field Theory (QFT) on a lattice is explained as a concrete example.In Sec. 3, the Monte Carlo technique is introduced and some numerical experiments are conducted.Sec. 4 is devoted to concluding remarks and discussion of future directions.

Truncation scheme
In this section, we specify the truncation scheme.

Coordinate-basis truncation for the single-boson system
Here, we introduce digitization in the coordinate basis [2,1,5], also called the field amplitude basis [6,8].Let us start with a review of the single boson field example.Let {|x⟩ |x ∈ R} be the coordinate basis for this particle that satisfies A simple way to digitize it is to introduce a cutoff to the value of the eigenvalue x as and introduce Λ points, The digitization parameters Λ, a dig , and R should be sent to infinity, zero, and infinity, respectively, to recover the action of the original operator x.By using |n⟩ to denote |x(n)⟩, we can write The momentum operator p appears in the Hamiltonian only in the form of p2 .A convenient way of regularizing it is The dimension of the Hilbert space is Λ and the truncated Hamiltonian is expressed as a Λ × Λ matrix.For each concrete example, we can diagonalize the Hamiltonian up to a rather large value of Λ and confirm that the digitization effects below a fixed energy scale disappear exponentially as ∼ e −cΛ with some c > 0. In terms of the number of qubits q, the truncation level is Λ = 2 q , and hence the suppression of the digitization effects is expected to be doubly exponential, ∼ e −c•2 q .

Coordinate-basis truncation for multi-boson system
Next, let us consider the more interesting case where we have multiple bosonic degrees of freedom.Suppose there are N bos variables ⃗ Note that we can take different R, Λ, and a dig for each bosonic coordinate x i ; here we use the same values for simplicity.For each i = 1, 2, • • • , N bos , the momentum operator pi appears in the Hamiltonian only in the form of p2 i .This is a natural extension of the single-particle case presented in the previous section.Again we choose a regularization The dimension of the Hilbert space is Λ N bos and the truncated Hamiltonian is expressed as a Λ N bos × Λ N bos matrix.Unless N bos is relatively small, it is difficult to directly determine the energy eigenvalues.Still, it is believed that the digitization effects below a fixed energy scale disappear exponentially as ∼ e −cΛ with some c > 0 [5,6].As far as we know, there is no proof for this scaling.(Here, an implicit assumption is that the wave function decays exponentially fast as |x| becomes large.)Note also that, even if the exponential accuracy of the digitization is valid, the precise dependence of the errors on a dig for each theory is not immediately clear.Therefore, it is important to have a method to make a quantitative estimation of the truncation effect that is applicable to a wide class of theories.Our solution to this problem in presented later in Sec. 3.

Scalar QFT on spatial lattice
As a particularly important example and prototypical QFT, we consider a scalar QFT on a d-dimensional spatial lattice.We consider the square lattice with equal lattice spacing a lat in all d directions.The lattice Hamiltonian is defined by Fields φ and π are dimensionless, and they correspond to the fields in the continuum theory by φ = a φcont.and π = a (d+1)/2 lat πcont. .Parameters such as mass are also made dimensionless, e.g., m lat = a lat × m.A vector ⃗ n lat ∈ Z d labels the lattice sites.This is different from ⃗ n used for the digitization, which we now denote by ⃗ n dig .The canonical commutation relation is imposed, i.e., The operators φ and π are the same as x and p in the previous sections: here we have used a different notation to more easily connect with the traditional symbols used in the physics community.μ is the unit vector along the µ-th dimension of the spatial lattice (µ = 1, • • • , d).This is different from the unit vector î used for the digitization (i = 1, • • • , N bos ).As infrared regularization, we introduce the periodic boundary condition with period L to all directions.Typically, the continuum limit (a lat → 0) is taken by fixing the physical volume a lat L. The number of bosonic degrees of freedom is N bos = L d .Hence the dimension of the truncated Hilbert space is Λ N bos = Λ L d .
We digitize φ just as before, by introducing R, Λ and a dig .We use the same digitization parameters for all lattice points for sake of simplicity.Note that a lat → 0 and a dig → 0 do not necessarily commute: the correct order is a dig → 0 first, then a lat → 0.

Monte Carlo estimate for truncation effect
In this section, we show that the MCMC methods [22] (see Ref. [23] for an elementary introduction) can be used to estimate the digitization effects.Numerical demonstrations will be provided as well.A short review of the MCMC methods is provided in Appendix A. The crucial point is that, as long as the problem under consideration reduces to an average over non-negative weights, the MCMC methods allow efficient computations.

Formulation and algorithm
Let us estimate the truncation effect in the coordinate-basis scheme by using Markov Chain Monte Carlo simulation.For N bos variables ⃗ x = (x 1 , • • • , x N bos ), we define integers (8).For a Hamiltonian defined by (2), we consider the thermal partition function defined by Here, β is related to the temperature T by β = 1 T .The trace is over the truncated Hilbert space.We rewrite it as where ⃗ n (j+1) , we can truncate the Fourier expansion of e −∆• i p2 i 2 .If we keep order-∆ terms, we obtain The crucial point that enables us to use MCMC methods is that the final form in ( 15) is non-negative for any ⃗ n (j) and ⃗ n (j+1) , if 1 − N bos ∆ a 2 dig > 0. In other words, we can write the partition function as a sum of non-negative weights.
Of course, when the digitization is removed and x i are treated as continuous variables, we can easily rewrite the partition function in the form of the standard Feynman path integral, which is well known to be amenable to MCMC simulations.In some sense, we are solving a simple problem in a complicated manner, to estimate the digitization artifact.On the other hand, we are using a mature classical numerical simulation technique to understand an issue arising in the new growing field of quantum simulations for quantum field theories.
Note that ∆ a 2 dig has to be small for (15) to be a good approximation.Therefore, as a dig becomes smaller, we need smaller ∆ and hence more Trotterization steps.For this reason, the classical simulation cost is sensitive to a dig , but not to Λ.In the demonstrations shown in the later sections, we take Λ and R very large so that the finite-R effect is negligible and we can focus on the finite-a dig effect.
Via the Monte Carlo simulation, expectation values of functions of ⃗ x can be estimated from stochastic samples.We write the expectation value as Here, ⃗ x and ⃗ n are related by (8).We use the approximation (15) for the simulation.In the limit of ∆ → 0, the approximation becomes exact, and expectation values at finite temperature can be obtained including the digitization effects.By dialing the temperature, the digitization effects at various energy scales can be studied.

Metropolis algorithm
In the Metropolis algorithm, the chain of configurations {⃗ n (1) , • • • , ⃗ n (K) } is generated in such a way that the probability distribution of configurations is proportional to the right-hand side of (13), where the approximation ( 15) is understood.More explicitly, the probability distribution is proportional to e −S(⃗ n (1) ,••• ,⃗ n (K) ;∆) , where exp −S(⃗ n (1) Note that lim In the following, we describe a naive update move for the Metropolis algorithm.The algorithm will start from an initial point, or configuration of all variables.For the initial configuration, we take all ⃗ n (j) to be the same and between 0 and Λ − 1. (Practically, we should take n ∼ Λ 2 , x ∼ 0.) We perform the following procedure to 1. Proposal: as a candidate for the new value of n 2. The candidate n ′ is automatically rejected (i.e., n (j) i remains unchanged) unless the following three conditions are satisfied: 3. Metropolis test: the candidate n ′ is accepted (i.e., n i becomes n ′ ) with the probability min 1, e −δS , where δS is the increment of the action.Otherwise the candidate n ′ is rejected (i.e., n (j) i remains unchanged).This is 'one sweep'.We repeat many sweeps and collect successive configurations along the Markov Chain.The above naive approach is not efficient because ⃗ x is changed only locally (i.e., only one time slice at each step), while ⃗ x should slowly change along the imaginary-time circle following dominant configurations.
In classical MCMC simulations of spin systems, local update rules for the Metropolis algorithm are known to perform poorly, in particular in regions of metastability or close to phase transitions.Inspired by cluster algorithms [26], we improve the algorithm by allowing the simultaneous update of a cluster as follows: 1. Choose 1 ≤ B ≤ B max randomly, where B labels the size of the "clusters" (or blocks) 2. Proposal: as a candidate for the next configuration, we vary n are proposed with probability 1 2 for each.Note that we use the same sign for all l's.3. The candidate n ′ is automatically rejected unless the following three conditions are satisfied for all l = j, • • • , j + B: |⃗ n ′(j) − ⃗ n (j−1) | = 0 or 1.
4. Metropolis test: the candidate n ′ is accepted with the probability min 1, e −δS , where δS is the increment of the action.Otherwise, the candidate n ′ is rejected.
B max can be any number between 1 and K.The optimal value depends on the detail of the Hamiltonian and digitization parameters and in principle it can be found via a standard analysis of the autocorrelation time (see e.g., Ref. [23]).As an example, one would first try the algorithm with B max = 1 and measure the autocorrelation time of the observable of interest.Then, one would increase B max by 1 unit and measure again the autocorrelation time.For some larger value of B max the autocorrelation time will start decreasing because the update sweeps become more efficient at finding new configurations.For the simulations reported in this paper, we set B max = K 2 , except for Fig. 2 in which B max = 1 is used for comparison.We chose B = 1, 2, • • • , B max with equal probability.

Simulation cost
Suppose that R is sufficiently large so that the support of the wave function is mostly contained in the truncated Hilbert space.As long as this loose condition is satisfied, the simulation cost is not sensitive to R. On the other hand, the cost is sensitive to a dig , because it appears in the expansion parameter ∆ a 2 dig in (15).To keep the expansion parameter small, we need to scale the number of Trotter steps K proportionally to a −2 dig .Therefore, the number of variables in the Monte Carlo simulation scales as This can be considered the scaling of the memory size, because we need this amount to store all the variables.Note that this memory scaling is much better than 2 N bos nq , which is required for the direct computation of the eigenvalues with exact diagonalization methods.For practical applications, n q does not increase with N bos (each bosonic degree of freedom is represented with the same number of qubits), and hence the necessary memory size increase only linearly with the system size.In addition to the memory scaling, the computational cost for one sweep scales polynomially with N bos and a dig .
Strictly speaking, it is difficult to establish the exponential suppression of the truncation effect unless a clear exponential behavior sets in at a small truncation level because exponentially large statistics and hence exponentially large computational costs are needed to have exponentially small statistical error. 4However, to establish that the error decays faster than a certain power, only a polynomially large cost is needed.

Single boson (N bos = 1)
As a sanity check of the method, let us consider the case of a single boson.In this case, the exact diagonalization method can be used to determine the low-energy spectrum rather accurately even for large truncation level Λ.Therefore, we can confirm that the values obtained from the Monte Carlo simulations described in the previous sections are correct.
We considered the Hamiltonian with the following quartic interaction: Note that we can consider negative m 2 as well.We studied λ = 1.0, m 2 = ±1.0 and a dig = 0.3, 0.5, 0.7.The step size ∆ is set to 0.001 so that ∆ a 2 dig is small and the error associated with the truncation of the Taylor expansion of e −∆• p2 2 is suppressed.To focus on the errors coming from a dig and not from R, we took R and Λ to be very large (specifically, Λ = 2001 and R = 1000a dig ).
The results are summarized in Table 1.'Exact' values are obtained by exactly diagonalizing the corresponding Hamiltonian in the truncated coordinate basis, with a sufficiently large value of Λ such that the first four nonzero digits are obtained precisely.Both methods give consistent results by taking into account the stochastic (statistical) error from the MCMC simulation.The simulation history for a dig = 0.5, m 2 = +1.0 is shown in Fig. 1.For  this plot, we set the maximum cluster size that is updated simultaneously to be B max = K 2 .The advantage of updating a large cluster simultaneously can be understood by comparing this plot with Fig. 2 corresponding to a simulation with B max = 1.To demonstrate the validity of the method, we consider the free theory with the potential: where m lat = a lat m in the Hamiltonian defined by (10).We focus on quantities that can be analytically computed for the case of the infinite-dimensional local Hilbert space (the full QFT without digitization effects) and we reproduce those results by taking the limit of a dig → 0 in our Monte Carlo simulations.(We use this particular setup because the analysis is simpler and the essence of the MCMC approach can be conveyed efficiently.We will comment on the cases without analytic results later in this section.) The Fourier transform on a lattice is defined by where We can write the free Hamiltonian in terms of φ⃗ q and π⃗ q as Ĥlat, free = ⃗ n where Each mode contributes ω ⃗ q 2 to the ground-state energy.The zero-point fluctuation of each mode is This exact value should be obtained in the limits ∆ → 0, R → ∞, a dig → 0 and T → 0. Note that we need to take both the digitization spacing to zero, and the energy scale to zero.On the other hand, at finite temperatures we have: To reproduce this result, we only take the limits ∆ → 0, R → ∞ and a dig → 0. We use this relation for the numerical demonstrations that follow.Note that the digitization and Fourier transform do not commute in general.Therefore, even for the free theory, the estimation of the digitization effect is not trivial.
The ground-state wave function for each mode φ⃗ q is Gaussian with the width 1 √ ω lat,⃗ q .For the higher-momentum modes, the widths are narrower, and smaller a dig will be needed for better approximation.
For numerical demonstration, we study a two-dimensional 4 × 4 lattice.The number of bosons is N bos = 4 × 4 = 16 and the dimension of the Hilbert space Λ N bos increases so quickly with Λ that numerical analysis on classical devices is practically impossible if the Hilbert space is constructed explicitly.
The exact value without digitization (∆ → 0 followed by a dig → 0) is 1 2 tanh(0.5)≃ 1.081977 and 1 6 tanh(1.5)≃ 0.184131.Simulations are conducted by using multiple (n stream ) streams with different random seeds.Each stream has n step ≥ 3.93×10 5 steps.The auto-correlation length for each ⃗ q is estimated by computing the integrated auto-correlation time for each stream, and for the largest length d ⃗ q , the initial 10d ⃗ q steps are discarded as a burn-in period, and the following steps are split into 10d ⃗ q -step samples.d ⃗ q < 10 3 for a ≤ 0.8, and for a = 0.9 and 1.0, stream sizes are greater than 2200d ⃗ q .The number shown in the parenthesis indicates the unbiased estimation of the standard deviation from the averages of at least ten streams.
We took ∆ 2a 2 dig to be 0.01 at maximum.Therefore, we expect the error associated with Trotterization is at most a few percent.The values we obtained are shown in Table 2.In Fig. 3, we plot Exact−MC Exact , where 'Exact' is the exact value (31) which should be obtained in the limit of ∆ → 0, R → ∞ and a dig → 0, and 'MC' is the numerical result in Table 2.The horizontal axis is 1 a dig .We can see that the error decreases to Exact−MC Exact ∼ 0.01, which is more or less the expected value of the Trotterization effect.
Very roughly, the width of the distribution of φ⃗ q is estimated by ⟨ φ⃗ q φ−⃗ q ⟩, which is √ 1.0819 ≃ 1.04 for ⃗ q = (q x , q y ) = (0, 0) and √ 0.1841 ≃ 0.43 for ⃗ q = (q x , q y ) = (π, π).A natural expectation is that the error associated with the digitization disappears exponentially fast once a dig becomes smaller than these values.Qualitatively, we can see such an exponential decrease in Fig. 3.
For the example discussed above, we knew the analytic results in the limit of a dig → 0.Even without knowing the analytic result, the digitization error analysis is straightforward, because our method works for any R and a dig .For simplicity, let us assume that R = ∞ and focus on the finite-a dig effect, as above.For various values of a dig , we can calculate  2, Exact−MC Exact is plotted, where 'Exact' is the exact value (31) which should be obtained in the limit of ∆ → 0, R → ∞ and a dig → 0 and 'MC' is the numerical results in Table 2.The horizontal axis is 1 a dig .We performed a fit of the data in Table 2 using the function φ⃗ q φ−⃗ q = Ae −B/a dig + C. In this plot we fix C to the exact value 1.081977 for ⃗ q = (0, 0) and 0.184131 for ⃗ q = (π, π).This is different from Fig. 4, in which C is treated as a fitting parameter.We obtained A = −0.079(11), B = 0.983(66) for ⃗ q = (0, 0) from a dig = 0.25, 0.30, 0.40, 0.50 and A = −0.0430(69),B = 0.794(58) for ⃗ q = (π, π) from a dig = 0.20, 0.25, 0.30, 0.40.2, we performed a fit using the function φ⃗ q φ−⃗ q = Ae −B/a dig + C. Unlike Fig. 3, here we treated C as a fitting parameter.
⟨O⟩, where O is observable under consideration, say ϕ 2 .Then, we can determine the a dig -dependence by fitting the numerical results at finite a dig values.If the analytic result at a dig = 0 is known, such a fit is easier (has less free parameters), but the fit can be conducted even without knowing the value at a dig = 0.In Fig. 4 we show how to perform a fit that includes the a dig = 0 value of the observable as a free parameter.We compare the fitted result with the known a dig = 0 result and find that they are statistically compatible.Alternatively, even if an analytical result is not available, one can determine the value at a dig = 0 (without any digitization of the degrees of freedom) by performing the standard Euclidean path integral computation via MCMC.

Conclusion and discussion
In this paper, we introduced an MCMC-based technique to determine the digitization effects in a class of bosonic systems, which have the Hamiltonian of the form (2), in the coordinate-basis truncation scheme.This technique enables us to study the expectation values of various operators at finite temperature including their digitization effects.By dialing the temperature, various energy scales can be studied.As a prototypical QFT application, we studied the (2 + 1)-dimensional scalar QFT regularized on a lattice.To check its convergence to the right answer, we studied the weak-coupling limit.The inclusion of the interactions is straightforward, as we have demonstrated in the case of a single-boson system.Our method can be used to estimate the resources needed for realistic quantum simulations, and also, to check the validity of the results of the quantum simulations.As a specific example relevant for the NISQ era, we can consider a variational algorithms to determine low-energy quantum states of a multi-dimensional lattice system.By using the MCMC method, one determines the ground-state energy of the digitized theory, including the scaling of the digitization effects.Whether this value can be reproduced is a good test of variational algorithms on NISQ devices when no other classical methods can be used due to the exponentially increasing requirements with the lattice size.Once we could find a quantum algorithm and quantum device that can pass this benchmark test, then we can use it to study more complicated observables which cannot be accessed by the MCMC method, such as excited-state energies or time-dependent correlation functions.
In this paper, we took the 'infrared'-cutoff parameter R very large and focused on the effects due to finite a dig .In near-term quantum simulations the finite-R effects can also be relevant.It is straightforward to study small values of R and estimate the finite-R effects, by using the same simulation technique we have proposed.
To apply the method introduced in this paper, each bosonic variable must take values in R.An interesting class of theories of this kind is Hermitian multi-matrix models, which are important in the context of the holographic duality.A nontrivial issue for those models is how to reconcile the digitization of the gauge field with the existence of the gauge symmetry.However, the study of the corresponding ungauged models [27,28,29], which are also relevant for the duality, is straightforward.It is not necessarily a problematic issue if one considers a quantum simulation on the extended Hilbert space, which is a rather common approach.For such an approach to work, one has to make sure that the low-energy states of the ungauged model in the extended Hilbert space are correctly described.Lattice gauge theories form another interesting class.Although unitary variables which appear in the Kogut-Susskind formulation [30] of lattice gauge theory cannot be studied by using the formulation given in this paper, it is straightforward to study the orbifold lattice [24,25].
Whether it is possible to study other digitization schemes by Monte Carlo methods is an interesting question.For SU(N ) gauge theories, digitization schemes based on discrete subgroups have been proposed [11,12,13,31].The effects of such schemes can be studied by using the standard lattice simulations [14,16] and in principle, it is possible to take the continuum limit along the time directions both for gauged and ungauged theories.While such digitization schemes lack systematic un-digitization limits, they might be sufficiently good tools in the NISQ era.Another interesting approach that can be studied by Monte Carlo is the use of non-commutative geometry.For example, Refs.[32,7] studied the truncation of the target space of the O(3) nonlinear sigma model by fuzzy sphere.
Another potentially useful direction is the use of variational Monte Carlo methods, specifically with the neural quantum states.Such an approach, which was confirmed to be valid for some theories before digitization [33,18,34,35], may work even with digitization, and the introduction of fermions is straightforward at least conceptually.The introduction of the singlet constraint is also straightforward, for example by adding the square of the generators of gauge symmetry to the action as a penalty term.The MCMC technique introduced in this paper could be used for cross checking purposes.

Figure 1 :
Figure 1: Simulation history of T = 0.1, a dig = 0.5, ∆ = B max = K 2 = 5000.The exact value obtained by exact diagonalization is 0.2539 in this case.

3. 3
Scalar QFT (N bos = L d ) Let us consider the scalar QFT on a square d-dimensional lattice of L sites in each direction.The number of bosons is N bos = L d , which quickly grows and makes it impossible to construct the Hilbert space explicitly on classical devices.

Figure 3 :
Figure 3: By using the values in Table 2, Exact−MC

Figure 4 :
Figure 4:By using the values in Table2, we performed a fit using the function φ⃗ q φ−⃗ q = Ae −B/a dig + C. Unlike Fig.3, here we treated C as a fitting parameter.