Griffiths-McCoy singularity on the diluted Chimera graph: Monte Carlo simulations and experiments on the quantum hardware

The Griffiths-McCoy singularity is a phenomenon characteristic of low-dimensional disordered quantum spin systems, in which the magnetic susceptibility shows singular behavior as a function of the external field even within the paramagnetic phase. We study whether this phenomenon is observed in the transverse-field Ising model with disordered ferromagnetic interactions on the quasi-two-dimensional diluted Chimera graph both by quantum Monte Carlo simulations and by extensive experiments on the D-Wave quantum annealer used as a quantum simulator. From quantum Monte Carlo simulations, evidence is found for the existence of the Griffiths-McCoy singularity in the paramagnetic phase. The experimental approach on the quantum hardware produces results that are less clear-cut due to the intrinsic noise and errors in the analog quantum device but can nonetheless be interpreted to be consistent with the existence of the Griffiths-McCoy singularity as in the Monte Carlo case. This is the first experimental approach based on an analog quantum simulator to study the subtle phenomenon of Griffiths-McCoy singularities in a disordered quantum spin system, through which we have clarified the capabilities and limitations of the D-Wave quantum annealer as a quantum simulator.


I. INTRODUCTION
Spin systems with disorder exhibit a number of unusual properties, a most notable example of which is the spin-glass state with randomly frozen spin configurations caused by disorder and frustration [1]. Even in the absence of frustration, disorder alone can lead to unexpected behaviors, and the Griffiths singularity is one of the most prominent examples [2]. It was shown by Griffiths [3] that the magnetic susceptibility of a randomly diluted ferromagnetic system shows singularities as a function of the external magnetic field within the paramagnetic phase. This singularity originates in the existence of rare but extremely large ferromagnetic clusters which behave almost like a pure ferromagnetic system responding strongly to the external field if the temperature is below the transition temperature of the non-diluted ferromagnetic system. The Griffiths singularity in the susceptibility is, however, a very weak essential singularity and is hard to observe experimentally.
Quantum effects significantly enhance the strength of this singularity, leading to divergences of the linear and nonlinear susceptibilities [2], a phenomenon known as the Griffiths-McCoy singularity [4,5]. Numerical, theoretical, and experimental investigations have been conducted on this subject, particularly in low-dimensional systems where the singular behavior of physical quantities is expected to appear most prominently [2,[6][7][8][9][10][11][12][13].
In the present paper we study this problem in the ferromagnetic transverse-field Ising model on the diluted quasi-twodimensional Chimera graph with disordered interactions by quantum Monte Carlo simulation on a classical computer and by quantum hardware simulations on the D-Wave quantum annealer. There are several reasons to motivate this direction of investigation. First, there exist few numerical or theoretical studies on the Griffiths-McCoy singularity for (quasi-)twodimensional disordered ferromagnets without frustration [11] although spin-glass systems have been investigated relatively extensively [6][7][8][9][10]12]. One expects that the disordered ferromagnet and spin glasses would behave qualitatively in the same way as far as the Griffiths-McCoy singularity is concerned because only disorder is relevant to this phenomenon, not the existence of frustration, but it nevertheless makes sense to confirm this conjecture explicitly numerically and experimentally by an additional concrete example. Second, it is an interesting exercise to compare the results of numerical simulations with data from the analog quantum simulator, the latter of which can be regarded as an experimental apparatus because real physical phenomena corresponding to the theory are expected to take place on the chip of the device if the latter operates as designed.
We find that the Griffiths-McCoy singularity exists in the present problem from numerical simulations. In contrast, data from the D-Wave device include a good amount of uncertainties due to noise, systematic bias, and other imperfections, but the results can be understood to be compatible with the existence of Griffiths-McCoy singularity, providing the first experimental test of the existence of this subtle singularity on the analog quantum simulator. This paper is organized as follows. We summarize basic facts about the Griffiths-McCoy singularity and list two quantities to be measured, linear and nonlinear susceptibilities, by numerical and experimental approaches in Sec. II. Methods and results of numerical simulations by quantum Monte Carlo are described in Sec. III. Experiments on the D-Wave device are discussed, and results are compared with those from nu-merical simulations in Sec. IV. Section V discusses the results. Additional details are described in Appendixes.

II. GRIFFITHS-MCCOY SINGULARITY
In this section, we first explain the basic idea of the original Griffiths singularity for the classical ferromagnetic Ising model on a randomly diluted lattice, and then discuss the quantum version, the Griffiths-McCoy singularity for the transverse-field Ising model.
Let us consider the classical ferromagnetic Ising model on a regular lattice, e.g., the two-dimensional square lattice, with a finite transition temperature between the paramagnetic and ferromagnetic phases. Suppose that we remove each bond (interaction) randomly with probability 1 − p and keep the bond with probability p. The ferromagnetic phase becomes gradually unstable as p decreases from 1 and the transition temperature T c (p) decreases as p decreases. At the percolation threshold p c for the given lattice, e.g., p c = 1/2 for the square lattice, the transition temperature reaches zero, T c (p c ) = 0.
Griffiths [3] proved that, in the temperature range within the paramagnetic phase but below the transition temperature at p = 1, T c (p) < T < T c (1), the magnetic susceptibility χ(h) as a function of the external field h is singular at h = 0. This behavior exists for any 0 < p < 1 both above and below the percolation threshold p c . The intuitive reason is that there exist very large clusters of almost ferromagnetically-ordered spins even for p < p c and they respond strongly to the external field if T < T c (1) because, in the infinitely large system, the spontaneous magnetization is positive for h > 0 and is negative for h < 0, i.e., a discontinuity at h = 0 (an infinite susceptibility). Since the probability of the existence of very large clusters is exponentially small as a function of their size, the resulting Griffiths singularity in the susceptibility is very weak, an essential singularity, and is therefore hard to detect experimentally.
Introduction of quantum effects significantly enhances the strength of the singularity, known as the Griffiths-McCoy singularity [2,4,5]. If we formulate the statistical mechanics of the transverse-field Ising model in terms of the Suzuki-Trotter decomposition [14], Ising spins along the Trotter (imaginary time) axis are coupled strongly by ferromagnetic interactions for moderate and weak values of the transverse field Γ. This strong ferromagnetic coupling along the Trotter direction causes the ferromagnetically coupled clusters to extended to the Trotter direction, greatly enhancing the response to the external field. Indeed, at zero temperature, the linear and nonlinear susceptibilities χ and χ nl are known to behave as [2] χ h d/z −1 , where d is the spatial dimension of the lattice and z is effective (running) dynamical exponent dependent on the value of Γ. Equations (2) and (3) indicate that χ and χ nl diverge at h = 0 if d/z < 1 and d/z < 3, respectively. The former inequality is indeed satisfied on the two-dimensional square lattice with the interaction and the transverse field chosen uniformly randomly, and the exponent z diverges at the transition point [11]. In one dimension, the singularity is stronger with z being divergent for any Γ larger than the transition point [15]. In three dimensions, the singularity is weak if any [12]. These properties hold for other types of disorder in the interactions, not just simple dilution of ferromagnetic interactions, including the Ising spin glass model under a transverse field, because only the existence of randomness is relevant in the sense of renormalization group and frustration does not play an essential role in the Griffiths-McCoy singularity.
In order to measure d/z to determine whether or not the linear and nonlinear susceptibilities diverge, it is useful to plot the histogram of measured values of local linear and nonlinear susceptibilities, which are known to behave as [7,8,10] log P (χ loc ) − d z which holds for large values of the variables. The slope of the plot gives the exponent. The local susceptibility χ loc is defined as where m i is the local magnetization computed from the quantum-mechanical expectation value of the spins and h i is the local longitudinal field applied only to site i. The local susceptibility χ ii depends on the site index i, and we collect the statistics of this quantity over i and random samples to generate the histogram P (χ loc ). The local nonlinear susceptibility χ nlloc is defined similarly as the third derivative of m i with respect to h i . In the paramagnetic phase, the correlation length is short and the ordinary (global) linear and nonlinear susceptibilities are expected to follow the same formulas as Eqs. (4) and (5). The local quantities are used often in numerical simulations because a larger number of samples can be generated than for the global quantities, leading to better statistics.

III. QUANTUM MONTE CARLO SIMULATION
In this section we describe the methods and results of quantum Monte Carlo simulations of the Griffiths-McCoy singularity on the diluted chimera graph as depicted in Fig. 1, in which some of the interactions on the original chimera graph are removed. This dilution of the chimera graph is expected to help us enhance the parameter range of the Griffiths-McCoy singularity since lower-dimensional systems with less degree of connectivity tend to show stronger effects of this singularity [12]. The chimera graph is chosen because the topology is directly realized on the D-Wave chip without the necessity of embedding. The above dilution of the chimera graph is not random, and it has no effect on the Griffiths-McCoy singularity except for the enhancement of parameter range as written above.
Randomness causing the singularity is in the choice of the interaction strength at the bonds remaining on the graph of Fig. 1, meaning that the J ij are uniformly chosen from the set {0, −0.2, −0.4, −0.6, −0.8, −1.0}. This problem on the quasi-two-dimensional diluted chimera graph has not been studied so far, and it is interesting to investigate it both by classical simulation, i.e., quantum Monte Carlo, and by a direct quantum simulation on the D-Wave quantum annealer, which may be regarded as an experiment on the analog quantum simulator. The present section concerns the former classical simulation.

A. Method
We perform quantum Monte Carlo simulations based on the Suzuki-Trotter decomposition with parameters listed in Table  I. The total number of spins we deal with for each Monte Carlo Table I for the definition of each symbol. For each random instance, we store the following physical variables: the absolute total magnetization, the squared total magnetization, the fourth moment of the total magnetization, the squared local magnetization, and the fourth moment of the local magnetization, where the brackets ... denote the statistical-mechanical average (i.e, the Monte Carlo average), which is expected to reduce to the quantum-mechanical average in the lowtemperature limit, and m stands for the sum of spin values over all spatial sites and along the Trotter direction, The quantity m i is the spin of local site averaged over the Trotter direction, We employ a GPU-based algorithm to accelerate the Monte Carlo simulations. We apply finite-size scaling to the analysis of the critical point and critical exponents through the Binder ratio g [16]: the global susceptibility χ, and the magnetization [ |m| ], where [· · · ] denotes the average over instances of randomness in interactions. As for the exponent d/z , we generate the histograms P (χ loc ) ad P (χ nlloc ) of the local susceptibility χ loc and the nonlinear local susceptibility χ nlloc , These quantities are expected to show the behavior described in Eqs. (4) and (5). Each histogram is generated from N × N rand samples, where N = 8L 2 is the number of sites. We also plot similar histograms for the global linear and nonlinear susceptibilities, to compare their behavior with the corresponding data for local susceptibilities. Each histogram is generated from N rand samples. We use the Python module "scipy.optimize.curve fit" in SciPy [17] for the data fittings. I: Parameters for the quantum Monte Carlo simulation. N rand is the number of random instances, M is the number of Trotter slices, L is the square root of the number of chimera units (N = 8L 2 is the total number of sites), Nstep is the number of Monte Carlo steps, βmax and βmin are, respectively, the maximum and minimum values of the inverse temperature, N β is the number of inverse temperatures, Γmin and Γmax are, respectively, the minimum and maximum values of the transverse field, and NΓ is the number of transverse field values.

B. Results
We first determine the transition point between the paramagnetic and ferromagnetic phases and then move on to the Griffiths-McCoy singularity within the paramagnetic phase. Figure 2(a) shows the Binder ratio as a function of Γ with the inverse temperature β = 20. We employ − log(1 − g) instead of the naive Binder ratio g in order to obtain a good resolution near the critical point [18]. The figure indicates that the phase transition point is located around Γ 1.7. Figure 2(b) shows the result of a finite-size scaling analysis of the same data. It is observed that the data with different sizes collapse onto the same curve for Γ c = 1.75(4) and ν = 1.4(2). We apply the same finite-size scaling to each temperature, and the results for the transition point and the critical exponent are summarized in Fig. 3. From Fig. 3(a), we observe that Γ c grows linearly as the temperature T approaches zero. Linear fitting shows that Γ c can be determined to be about 1.78 in the limit of zero temperature. We also see that the critical exponent ν does not clearly depend on the temperature and can be estimated as ν 1.5. One may notice that the error bars are relatively large near zero temperature. This large uncertainty comes from the fact that the data collapse in Fig. 2(b) does not depend very much on the values of Γ c and ν: the data collapse is robust against the change of the values of Γ c and ν, especially the latter. We think it reasonable to assume Γ c = 1.78(2) based on the data in the temperature range T > 0.1 where the data are relatively stable. Estimation of ν and other critical exponents such as β and γ involves large uncertainties as was the case in three-dimensional spin glasses [19], as detailed in Appendix A. At least, the estimated transition point Γ c = 1.78 gives consistent results in finite-size scaling of other physical quantities such as the susceptibility and magnetization. See Appendix A.

Histograms of susceptibilities
In this section we estimate the exponent d/z , which is critical to determine the existence of the Griffiths-McCoy singularity, from the data of local and global susceptibilities. Figure 4 shows the histogram of the local susceptibility P (χ loc ) in a paramagnetic region far from the critical point (Γ = 1.895). In this region, the size dependence of the data is weak except for the tail of the distribution where the probability is small due to insufficient statistics. Since Eq. (4) is valid for large χ loc , we have to carefully choose the range to fit the data to Eq. (4). Using the data in the range between slightly below the peak and P (χ loc ) 10 −4 , d/z is estimated to be around 11.2. According to the discussion in Sec. II, we conclude that both linear and nonlinear susceptibilities exhibit no divergence at this value Γ = 1.895 since d/z > 3. Figure 5 is for Γ = 1.79, closer to the critical point Γ c = 1.78. We observe that the slope clearly depends on the system size due to the finite-size effect and the slope tends to be shallower as the system size increases. To extract the information on d/z , we use the data for the largest size L = 12 since this is expected to give a lower bound of the slope in the large-size limit. We thus conclude that the exponent d/z is smaller than about 4.0.

Local susceptibility
The histogram P (χ loc ) for data in the ferromagnetic phase Γ c = 1.4 is in Fig. 6. We find the behavior is quite different from previous cases for the paramagnetic phase. In the ferromagnetic phase, P (χ loc ) a slightly increasing function especially for the large system size.
We analyze the nonlinear local susceptibility χ nlloc similarly. The results for Γ = 1.895 (far from the critical point)  Fig. 4. The black solid line shows a fit to the data with the largest system size L = 12 whose slope is −5.40 ± 0.33, which indicates that d/z is at least smaller than about 4.40 ± 0.33.   Fig. 9. Although there exist uncertainties in these values, it is useful to take into account the fact that the plotted values of the exponent d/z give upper bounds. We then observe the plausibility that d/z becomes smaller than the threshold d/z = 3 before the critical point is reached, suggesting that the nonlinear susceptibility diverges within the paramagnetic phase. This behavior is consistent with the previous study for the case of continuous distributions of random ferromagnetic interactions and random transverse field for a system on the square lattice [11]. More subtle is the divergence of the linear susceptibility since it is difficult to determine from the data whether or not d/z becomes smaller than 1 in the paramagnetic phase Γ > Γ c . . The black solid line shows the fitting to the data whose slope is −4.02 ± 0.06, which indicates that d/3z is around 3.02 ± 0.06. . The black solid line shows the fitting to the data with the largest system size L = 12 whose slope is −1.29±0.08, which indicates that d/3z is smaller than 0.29±0.08. The exponent d/z as a function of Γ measured at the inverse temperature β = 50. The data in circles and crosses are taken from χ loc and χ nlloc , respectively. The vertical dashed line shows the critical point in the zero temperature limit Γc = 1.78. The critical point corresponding to this finite-temperature data would be smaller than Γc = 1.78 for the zero-temperature value. The horizontal dotted line is for d/z = 3, where the nonlinear susceptibility starts to diverge. Notice that the data for d/z are upper bounds.

Global susceptibility
We next verify if the data for the global susceptibility are consistent with those for the local susceptibility. Figures 10,  11, and 12 show the histogram of the global linear and nonlinear susceptibilities. It is to be noticed that we have less data points than in the case of the local susceptibilities. The resulting value of d/z is plotted in Fig. 13.
From these results we confirm that the data for the global susceptibilities exhibit similar behavior to those for the local susceptibilities both in the histograms and the exponent z . This is important because the experimental data from the D-Wave device are available only for the global susceptibilities.
Lastly, we point out that the data for the global nonlinear susceptibility at Γ = 1.865 as shown in Fig. 14 gives the value of the exponent d/z 1.3 much smaller than the value indicated in Fig. 13, around 4. It may be due partly to insufficient statistics but further scrutiny is needed.
where s is the time parameter running from 0 to 1, and the time dependence of A(s) and B(s) is depicted in Fig. 15. To generate a state as close as possible to a quantum thermal equilibrium state on the D-Wave machine, we employ the annealpause-quench protocol where we first perform an anneal up to value s = s * , pause at this point for a while, then quench to s = 1 as rapidly as possible. By applying this protocol, the FIG. 14: P (χ nl ) which is similar to Fig. 7 but for the global nonlinear susceptibility at Γ = 1.865. The black solid line shows a fit to the data whose slope is −1.616 ± 0.25, which indicates that d/3z is around 0.616 ± 0.25 or d/z 1.85 ± 0.75, which seems to be quite small compared to the data from global linear susceptibility in Fig. 13. D-Wave device may return the spin configuration σ sampled from a distribution not far from the canonical ensemble with the Hamiltonian H(s * ) of Eq. (16) although arguments exist on the accuracy of this procedure [20,21].
We use the same diluted Chimera graph as in Sec. III for the experiment on the D-Wave machine. We set the amplitude of interactions to a half of that used in the quantum Monte Carlo method, which is expected to reduce the effect of the noise in the D-Wave machine since a large amplitude of interactions tends to amplify analog errors on the machine [22]. The parameters used throughout the experiments are listed in Table II. First we describe how to measure the magnetization. For each annealing process, we obtain a set of classical values of spins {σ 1 , σ 2 , · · · } following the anneal-pause-quench protocol with longitudinal field switched off and store the magneti- II: Parameters for experiments on the D-Wave machine. N rand is the number of random instances, N rand is the number of random instances for the histogram of susceptibility, Nrep is the number of annealing repetitions, L is the square root of the number of Chimera units (8L 2 is the total number of sites to be denoted as N ), s * min is the minimum value of annealing schedule, s * max is the maximum value of the annealing schedule, NΓ is the number of transverse field values, t1 is the anneal time in the anneal-pause-quench protocol, t2 − t1 is the pause time, t f − t2 is the quench time.
This annealing process is repeated N rep times to obtain a set of magnetization values m a (1 ≤ a ≤ N rep ). Given enough interval time (200µs) between consecutive annealing processes, we can assume that samples of the magnetization m a are uncorrelated with each other. The nth moment of the magnetization is calculated from these as Various physical quantities such as the Binder ratio are derived from the above nth moment of magnetization.
To estimate the linear and nonlinear susceptibilities χ and χ nl , we measure the magnetization as a function of the longitudinal field h. Linear and nonlinear susceptibilities are obtained by applying a polynomial fit to the magnetization curve, for each given instance of random interactions. Due the analog nature of the D-Wave device, naive experiments without noise mitigation produce data with very limited reliability, in particular in the present case of the detection of a delicate phenomenon. We therefore apply two kinds of techniques, calibration of individual flux bias and the standard gauge averaging, to reduce noise for more reliable results. A technical description of the former method is given in Appendix B.

B. Results
We first analyze the data to determine the transition point and then estimate the exponent d/z . shows that the magnetization is distributed around zero below a certain value of s * close to 0.39 and tend to have a peak near ±1 above this threshold point. This change of the shape of the histogram suggests the existence of a phase transition at around s * c 0.39. Also, especially for the small system size, we observe that some samples have large values close to ±1 even in the region supposed to be the paramagnetic phase. One of the possible reasons for this unexpected behavior is a systematic error coming from the quench process during the anneal-pause-quench protocol [20]: At the end of the protocol, the annealing parameter s is changed from s * to 1 as quickly as possible, which corresponds to a "switchingoff" of the transverse field. Although this quenching process is supposed to be performed quickly, the elapsed time in this process ( 1µs) may still be long enough to affect the final state of the system, leading to a broad distribution of the magnetization in the presumed paramagnetic phase, which should not be the case in theory. A similar effect was reported in Ref. [23] where the histogram of the magnetization exhibits a broad distribution in the paramagnetic phase, similarly to our case. The data for the small size L = 6 in the paramagnetic phase are strongly affected by this imperfection and we should take sufficient care in the analysis.

Binder ratio and averaged magnetization
The results of the Binder ratio and the averaged magnetization m av = [ |m| ] are shown in Fig. 17. We usually expect the Binder ratio with different system sizes to cross at the transition point as in Fig. 2. However, in the present case, we find no crossing point in Fig. 17, and both the Binder ratio and the averaged magnetization have finite values even in the paramagnetic region most prominently for the small system size. This behavior can be understood in the same way as in the case of the histogram of the magnetization: The relatively slow quenching process in the anneal-pause-quench protocol may allow the system to follow the decrease of the transverse field, driving the system toward ferromagnetic ordering. The large values of the Binder ratio and magnetization at small s * , especially for L = 6, should reflect the broad distribution of magnetization seen in Fig. 16. Nevertheless, we at least find that the magnetization for the largest system L = 16 is likely to have an inflection point around s * 0.39, remotely suggesting the existence of a phase transition point around this value. The horizontal axis s * denotes the pause point during the anneal-pause-quench protocol, corresponding to a finite Γ of the transverse field. The vertical axis is for the magnetization. Up to a certain point s * 0.39, the magnetization tends to be distributed around zero, which implies that the system is in the paramagnetic phase. In the region with s * larger than 0.39, the magnetization tends to have values of saturation ±1, and the system is in the ferromagnetic phase. Note that the color code is in logarithmic scale.  Figure 18(a) shows the global linear susceptibility χ for three system sizes. We find that there are peaks at around s * 0.39 and the position becomes smaller as the system size increases. Extrapolation to the infinite-size limit as shown in Fig. 18 gives s c = 0.386 ± 0.002, which is consistent with the data of the histogram of the magnetization in Fig. 16. Since the data, especially for the smallest size L = 8, are not necessarily very reliable, we choose not to go beyond the estimate of the approximate value of the transition point and avoid finite-size scaling analysis of critical exponents.

Global linear and nonlinear susceptibilities
We also measured the global nonlinear susceptibility as shown in Fig. 19(a), where the aspect ratio is the same as in Fig. 18, for convenience of comparison. Although the peak is a little bit narrower compared to the global linear susceptibility, the peak position is almost the same as in the global linear susceptibility, leading to the same critical point s c = 0.386 ± 0.002 in the infinite system size limit as shown in the bottom panel of Fig. 19.

Estimation of the exponent
We now analyze the histogram of susceptibilities to estimate the exponent d/z . Figure 20 shows the histogram of the global linear susceptibility with the pause point s * = 0.365 which is far from the critical point and is expected to be in the paramagnetic phase. Notice that when applying a linear fit to estimate d/z , we exclude the data with P (χ) ≤ 10 −2 partly because these data may not be reliable due to insuffi-   cient statistics. Another reason is given below. We find that a linear fit to the data does not seem too bad in the region of large susceptibility, especially for the largest system size.
We also observe a long tail of small susceptibilities at the left part of the graph, notably for the small size. This tail may be understood by considering the feature of the histogram of the magnetization in Fig. 16, where samples with values near saturation 1 exist even in the paramagnetic region. These samples may have small susceptibilities because the system does not respond to the field when the magnetization is close to saturation, resulting in the long tail of small susceptibility at The black solid line shows a fit to data whose slope is −12.3 ± 2.0, which indicates that the exponent d/z is around 11.3±2.0. The black solid line shows a fit to data whose slope is −9.7 ± 2.1, which indicates that the exponent d/z is around 8.7 ± 2.1. the left part of Fig. 20. Samples with large negative magnetization close to −1 would respond very strongly to the positive field h > 0, flipping the state from m −1 to m 1, and the tail of the distribution for very large χ would correspond to such cases. We therefore drop the data in the right-most tail of distribution from the analysis. In contrast, samples with small magnetization in the paramagnetic phase are likely to have reasonable properties, which would yield the moderately large susceptibility compared to the case with saturated magnetization. It therefore seems reasonable to consider that analyses of data with moderately large susceptibilities would give relatively reliable results. We then apply a linear fit to the data with large, but not too large, values of the susceptibility. We find that d/z is around 13.7 for the present pause point s * = 0.365.
The distributions P (χ) at the pause points s * = 0.369 and 0.375, which are closer to the critical point but still in the paramagnetic phase, are shown in Figs. 21 and 22, respectively. The resulting exponents are d/z 11.3 ± 2.0 and 8.7 ± 2.1 for s * = 0.369 and 0.375, respectively. In contrast, the data for larger s * are difficult to analyze to extract the exponent reliably. The data for the region above the critical  point, s * = 0.394, is shown on Fig. 23. No simple power-law decay is observed here. We find that P (χ) monotonically increases as χ increases. The same behavior is observed for the local linear susceptibility obtained by quantum Monte Carlo as shown in Fig. 6.
Similar behavior is observed in the global nonlinear susceptibility. Figure 24 shows the histogram P (χ nl ) for s * = 0.365. The histogram of the global nonlinear susceptibility with the pause point s * = 0.369 is shown in Fig. 25. We observe that d/3z decreases monotonically as in the global linear susceptibility, about 7.0 ± 1.1 (s * = 0.365) and 4.8 ± 0.9 (s * = 0.369).
The relation between the exponent d/z and the pause point s * is shown in Fig. 26. Although data points for s * beyond 0.375 are excluded because of low reliability of the estimation of the slope of the histogram, the tendency seems consistent with the assumption that the inequality d/z < 3 for the divergence of the nonlinear susceptibility is satisfied within the paramagnetic phase. . The black solid line shows a fit to data whose slope is −5.8 ± 0.9, which indicates that the exponent d/3z is around 4.8 ± 0.9.
where β phys 12mK is the physical temperature of the D-Wave chip and H(s) and H denote the Hamiltonian of the D-Wave device and the quantum Monte Carlo in Eqs. (16) and (1), respectively. Assuming that the two exponentials in Eq. as follows, or with physical units written explicitly, where h and k B denote the Planck constant and the Boltzmann constant, respectively, the former not to be confused with the external field. Figure 27 shows the relations Eqs. (24) and (25). From this figure we read that the critical point s * = 0.386 corresponds to the parameter β = 2.49 and Γ = 1.37 in the quantum Monte Carlo method, which is not far from Γ 1.6 with the temperature β = 2.49 (T = 1/β 0.4) according to Fig. 3(a). Although perfect quantitative agreement has not been expected, we have reached a reasonable degree of agreement.

V. DISCUSSION
We have carried out quantum Monte Carlo simulations and experiments on the D-Wave quantum annealer in order to in-vestigate if the Griffiths-McCoy singularity is observed in the transverse-field Ising model with random ferromagnetic interactions on the diluted Chimera graph. The results of quantum Monte Carlo indicate it to be very likely that there exists a parameter range within the paramagnetic phase where the local and global nonlinear susceptibilities diverge, implying the existence of the Griffiths-McCoy singularity. It is difficult to determine with confidence from our data whether or not the local and global linear susceptibilities diverge in the paramagnetic phase although it can well be the case if the present model belongs to the same universality class as the transversefield Ising model on the square lattice with uniformly random ferromagnetic interactions and uniformly random transverse field [11].
The data from the D-Wave device include a larger amount of uncertainties than those of quantum Monte Carlo due to the systematic bias in flux qubits and other sources of errors that are intrinsic to the analog quantum device. In particular, the distribution of magnetization has a significant amount of data points near saturation ±1 within the paramagnetic phase even after careful calibrations to cancel the ferromagnetic bias. Nevertheless, the moderately-large-value part of the distribution of linear and nonlinear susceptibilities can be considered fairly robust against the bias and errors because a state with saturated magnetization 1 has no room of further change toward larger values of magnetization and therefore will respond only weakly to the external field, contributing very little to the large-value part of the distribution of susceptibility. The very-large-value tail of the distribution may also be discarded for a similar reason and for insufficient statistics. With this observation in mind, we suppose that the exponent d/z estimated from the part of moderately large values of susceptibility is relatively reliable. In this way, we may conclude that the data can be interpreted to be consistent with the statement that d/z < 3 is satisfied in the paramagnetic phase, meaning the existence of the Griffiths-McCoy singularity. If this is indeed the case, the present study is the first case in which this very subtle phenomenon involving rare regions of ferromagnetic clusters, enhanced by quantum effects, has been observed in an analog quantum simulator. Further developments in device technologies and tools of data analyses will lead to improved reliability, and possibly promoting research activities toward quantum simulation of complex many-body systems by analog quantum simulators. thors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, DARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. The research of K.N. was supported by JSPS KAKENHI Grant Number JP18J13685.

Appendix A: Critical exponents
This appendix discusses the estimation of critical exponents β, γ, and z at the ferro-para transition point from the data of quantum Monte Carlo.
We have carried out finite-size scaling analyses of the global susceptibility χ and the magnetization |m| as in Figs. 28 and 29. Estimated values of the exponent ν and the transition point Γ c described in Sec. III B give a reasonable collapse of the data. The finite-temperature values of the exponents γ and β are extrapolated to zero temperature as in Fig. 30. Although the extrapolation suggests γ to be close to 0.94 and β to 1.0, these figures show large uncertainties in these estimates.
We next try to estimate the dynamical exponent z. The finite-size scaling of the Binder ratio is given by  Figure 31 shows this plot. From these figures it is apparently difficult to determine z without large ambiguities, and therefore we instead employ the following "mean square error" criteria. Figure 31 suggests that the functionf (β/L z ) can be approximated by a quadratic form, where a < 0 as seen in Fig. 31. By fitting the parameters a, b and c to minimize the mean square error between the quadratic function and the actual data, we estimate the optimal parameters a, b and c and the corresponding mean square error. The mean square error is depicted in Fig. 32, which shows z 1 is the best estimate, consistently with a previous study [24].

Appendix B: Calibration of individual flux bias
We describe the method of calibration of qubits in the D-Wave device and its consequences. Without any calibrations, as shown in Fig. 33, some of the flux qubits tend to have the spin-up direction while others tend to have the spin-down direction even without interactions and longitudinal field. This intrinsic bias can be modeled as an effective longitudinal local field, and we try to eliminate its effect by calibration. On the D-Wave device, this effective local field can be canceled to some extent through the "flux-biases" option, and the task is to search the optimal flux bias for each flux qubit such that the average of a single isolated spin becomes zero. The binary search algorithm can be used to find the optimal flux bias for each qubit. An example is listed in Algorithm 1.  set flux bias of ith qubit to h flux i for each qubit 8: set interactions Jij and local fields hi to zero 9: perform anneal-pause-quench protocol on the D-Wave machine 10: calculate thermal average of spin mi for each qubit

Zero interactions
We next show how these error mitigation techniques affect the data from the D-Wave device. In the case of no interactions between qubits, We first prepare the Hamiltonian with all interactions J ij set to zero. We next measure the spin configuration for each qubit averaged over 100 runs of annealing, which is expected to take the value near zero. Figure 33 shows m i for each site before and after applying the error mitigation techniques. Purple dots and green dots show the result before and after applying calibration, respectively. We observe that m i of the data without error calibration tend to have values far from zero even if the interactions are set to zero. After error calibration, m i has much smaller deviations and the values are closer to zero than the case without calibration.

Diluted Chimera graph
Next we apply the above technique to the diluted Chimera graph and see how the result changes by calibration. Figure  34 shows the histogram of the magnetization without calibration. We see that there are unphysical multiple peaks in the ferromagnetic region, where the magnetization is expected to have peaks close to ±1. Although we see a symptom of phase transition around s * = 0.35, it is quite hard to reliably extract information from this set of data. Figure 35 shows the data after calibration. We clearly observe peaks only at ±1 in the ferromagnetic region, as they should. From these data we confirm that the calibration process is effective and indispensable for reliable quantum simulations on the D-Wave device.