Near-term distributed quantum computation using mean-field corrections and auxiliary qubits

Consider the class of spin models:

$$\begin{align} H = -\sum_{\langle i,j \rangle} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} - \sum_{i = 1}^{N} h_{i} \hat{S}^{\left(i\right)}_x . \end{align} \tag{ 1 }$$

Here, $\hat{S}^{(i)}_{\alpha}$ and $\hat{S}^{(j)}_{\beta}$ are spin-1/2 spin operators acting on sites i and j, where $\alpha, \beta \in \{x, y, z\}$. The coefficient $J_{ij}^{\,\alpha, \beta}$ gives the strength and sign of the interaction. For concreteness and without loss of generality, we have selected transverse fields h_i to point along the x-axis. The Hamiltonian acting strictly within some sub-system f will neglect any operators acting outside of f, yielding

$$\begin{align} H^{\left(\,f\,\right)} = -\sum_{\langle i,j \rangle \in f} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} - \sum_{i\in f} h_{i} \hat{S}^{\left(i\right)}_x , \end{align} \tag{ 2 }$$

the bare Hamiltonian that acts within a fragment f when no corrections are included.

Clearly, the simple exclusion of interactions that span the interface between f and its environment (i.e. the fragmented evolution of f under $H^{(f)}$) will, in general, poorly approximate the evolution of the sub-system under the full Hamiltonian. The fragment qubits will behave as a closed system without external interactions. Although generally these interactions cannot be exactly simulated without modeling all of the system's spins on a single fragment, we introduce a mean-field to partially capture the action of each missing interaction. Mean-field methods have frequently been used to simplify the simulation and study of quantum systems, and statistical physics [44, 58, 59]. Here, the strength and sign of the introduced mean-field correction is informed by the measurement of the corresponding environment spin, while the correction's axis is determined by that of the corresponding interaction's spin operator that would act within fragment f. The resulting mean-field corrected Hamiltonian is given by:

$$\begin{align} H^{\left(\,f\,\right)}_\textrm{MF} = -\sum_{\substack{\langle i,j \rangle \in I,\\i \in f}} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \langle \hat{S}^{\left(j\right)}_{\beta} \rangle - \sum_{\langle i,j \rangle \in f} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} - \sum_{i \in f}h_{i} \hat{S}^{\left(i\right)}_x . \end{align} \tag{ 3 }$$

The strength and direction of the mean-fields appearing in $H^{(f)}_\textrm{MF}$ should be updated regularly to reflect the current state of the environment spins. Physically, this requires regular mean-field measurements of the fragments. Evolution must therefore be reset to the initial state in order to proceed by one time step dt, with each new mean-field measurement being stored to progress the evolution. The process of incrementing the time evolution by one time step per simulation is commonly implemented in order to track the time dynamics of an observable [60], resulting in a complexity that scales as $\mathcal{O}(N_t^2)$ in the number of time steps N_t .

For nearest-neighbor spin models (e.g. the transverse field Ising model [49]), the selection of a target spin for each auxiliary is somewhat trivial, as at most two qubits interact with a fragmented section of the chain. The choice of auxiliary qubit encoding may be unclear for more general systems. Here, we present a method for auxiliary target qubit selection that yields, on average, the optimal auxiliary qubit encoding. Specifically, we consider how auxiliary target selection affects the simulation error to the first non-vanishing order in dt. This simulation error arises from the omission of interactions forming the interface of some particular fragment f and the remaining environment spins E. The full derivation of the leading error is provided in appendix A.2; here, we sketch the derivation and build on the result.

To derive the leading error, we define the fidelity between the evolved state of a fragmented system and that of the full system to be:

$$\begin{align} F\left(t\right) = |\langle \Psi| U^{\dagger}\left(t\right) U_{I}^{\left(\,f\,\right)}\left(t\right) |\Psi \rangle|^2 . \end{align} \tag{ 4 }$$

The unitary operator $U(t) = \exp(-iHt)$ evolves the system exactly under the full Hamiltonian H, while $U_{I}^{(f)}(t) = \exp(-iH_\textrm{I}^{(f)}t)$ evolves the system under a fragmented Hamiltonian. Notably, the Hamiltonian $H_\textrm{I}^{(f)}$ around which this discussion revolves is not the Hamiltonian that acts only within fragment f; rather, $H_\textrm{I}^{(f)}$ additionally includes all interactions between qubits that are external fragment f. Thus, $H_\textrm{I}^{(f)}$ lives in the Hilbert space of the full system, neglecting only the interactions crossing the interface of fragment f in order to isolate the fragmentation error associated with f. In appendix A.2, equation (4) is expanded for short times t to understand how the evolution error $\epsilon(t) = 1 - F(t)$ depends on the strength of the neglected interactions. Through the use of Taylor expansion and the Baker–Campbell–Hausdorff (BCH) formula [61], we arrive at the first non-vanishing correction to the fidelity:

$$\begin{align} F\left(t\right) \approx 1 - \text{var}\left( H - H_\textrm{I}^{\left(\,f\,\right)}\right) t^2 , \end{align} \tag{ 5 }$$

where $\text{var}(\mathcal{O})$ is the quantum variance of operator $\mathcal{O}$. The error $\epsilon(t) = 1 - F(t)$ is thus given by $\text{var}( H - H_\textrm{I}^{(f)}) t^2$ for short times t.

The form of the short-time error provides a simple rule for choosing the target auxiliary qubits for fragment f to minimize error; namely, select the environment qubit(s) whose interactions contribute most significantly to the variance $\text{var}(H - H_\textrm{I}^{(f)})$. This choice will minimize the short-time error of evolving the state by the fragmented Hamiltonian, which will lead to higher fidelity performance, on average (see section 3.3). Moreover, if the auxiliary selection is updated sufficiently often, the selection becomes exact as the short-time error dominates from the time of one auxiliary encoding to the next.

Although the final form of the short-time evolution error provides insight into optimal auxiliary selection, the procedure for estimating a particular qubit's contribution to the error within the distributed framework is less straightforward. For a general spin Hamiltonian, this variance is given by:

$$\begin{align} \begin{aligned} \text{var}\left(H - H_\textrm{I}^{\left(\,f\,\right)}\right) & = \text{var}\left( -\sum_{\langle i,j \rangle \in I} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta}\right)\\ & = \sum_{\langle i,j \rangle \in I} \sum_{\langle i^{^{\prime}},j^{^{\prime}} \rangle \in I} \sum_{\alpha, \beta} \sum_{\alpha^{^{\prime}}, \beta^{^{\prime}}} J_{ij}^{\,\alpha, \beta} J_{i^{^{\prime}}j^{^{\prime}}}^{\,\alpha^{^{\prime}}, \beta^{^{\prime}}} \left( \left\langle \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} \hat{S}^{\left(i^{^{\prime}}\right)}_{\alpha^{^{\prime}}} \hat{S}^{\left(j^{^{\prime}}\right)}_{\beta^{^{\prime}}}\right\rangle - \left\langle \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} \right\rangle \left\langle \hat{S}^{\left(i^{^{\prime}}\right)}_{\alpha^{^{\prime}}} \hat{S}^{\left(j^{^{\prime}}\right)}_{\beta^{^{\prime}}} \right\rangle \right) . \end{aligned} \end{align} \tag{ 6 }$$

Estimating the full variance of equation (6) requires 4-point correlation measurements. If the distributed simulators are linked solely via classical channels, correlation measurements are only accessible when all relevant qubits are local to a single fragment. This implies that two auxiliary qubits—one targeting j and one targeting j'—must already be placed within the fragment in order to access the required 4-point correlator measurements. For N_E environment qubits, there are $\mathcal{O}(N_{E}^{2})$ combinations, requiring $\mathcal{O}(N_{E}^{2})$ copies of the system in order to estimate all required 4-point correlators, undermining (although not necessarily precluding) the motivations for fragmented quantum simulation with such a technique.

The correlator calculation simplifies significantly when the variance is calculated with respect to a known product state, but a new issue arises: for many spin model Hamiltonians, the variance will vanish for certain initial product states. In fact, for the case of the transverse field Ising model [49], this quantity vanishes for all computational basis states, providing no insight into the proper auxiliary choice.

We propose a two-part solution that addresses these issues. First, we propose a proxy v(a) that estimates the contribution of one potential auxiliary a to the variance:

$$\begin{align} v\left(a\right) = \sum_{\langle i,j \rangle \in I} \sum_{\langle i^{^{\prime}},j^{^{\prime}} \rangle \in I} \sum_{\alpha, \beta} \sum_{\alpha^{^{\prime}}, \beta^{^{\prime}}} J_{ij}^{\,\alpha, \beta} J_{i^{^{\prime}}j^{^{\prime}}}^{\,\alpha^{^{\prime}}, \beta^{^{\prime}}} \delta_{j,a} \delta_{j^{^{\prime}},a} \left( \left\langle \hat{S}^{\left(i\right)}_{\alpha}\hat{S}^{\left(j\right)}_{\beta} \hat{S}^{\left(i^{^{\prime}}\right)}_{\alpha^{^{\prime}}} \hat{S}^{\left(j^{^{\prime}}\right)}_{\beta^{^{\prime}}}\right\rangle - \left\langle \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} \right\rangle \left\langle \hat{S}^{\left(i^{^{\prime}}\right)}_{\alpha^{^{\prime}}} \hat{S}^{\left(j^{^{\prime}}\right)}_{\beta^{^{\prime}}} \right\rangle \right) . \end{align} \tag{ 7 }$$

Inserting the two Dirac delta functions $\delta_{j,a} \delta_{j^{^{\prime}},a}$ eliminates the cross-terms in equation (6) that depend on multiple environment qubits. Thus, a single auxiliary is required to estimate v(a), and in total $\mathcal{O}(N_{E})$ partitions are required to acquire v(a) for all potential auxiliary targets. In addition to requiring fewer measurements, this proxy focuses on a's contribution to the variance while neglecting the cross-terms that involve contributions from other potential auxiliary qubits. Secondly, to avoid scenarios where the variance vanishes for initial product states, we suggest first evolving the system for one time step dt for a particular choice of a before estimating v(a). Although this procedure is more involved than calculating v(a) for the initial product state directly, the overhead remains linear in the number of potential auxiliary qubits.

We first focus on the scheme that involves no quantum information transfer. In this scheme, the auxiliary qubits are selected at the beginning of the simulation and fixed to target a single environment qubit throughout the evolution. We refer the reader to figure 1 and appendix A.1 for an in-depth look at how a system is fragmented for time evolution. As a representative example, consider the transverse field Ising model (TFIM) [48, 49] with a uniform transverse field:

$$\begin{align} H_\textrm{TFIM} = - J \sum_{\langle i,j \rangle} \hat{S}_{z}^{\left(i\right)} \hat{S}_{z}^{\left(j\right)} - h \sum_{i}^{N} \hat{S}_{x}^{\left(i\right)}. \end{align} \tag{ 8 }$$

This system has been studied in depth to better understand the physics of quantum phase transitions [62–64]. We consider the evolution of a quantum system initialized in the computational basis state $|\boldsymbol{0}\rangle$ under $H_\textrm{TFIM}$, implementing exact unitary evolution numerically using PennyLane [65] with mean-field measurements updated every $\textrm{d}t = 0.1 / J$. Figure 2 displays the scheme's performance for a 12-qubit model with nearest neighbor interactions of J = 1.0. The results presented in figure 2(a) are averaged over non-zero transverse fields h ranging from ±1, while figure 2(b) features the specific case of h = 1.0. In figure 2(a), the system is split into two fragments, each simulating six of the system qubits and some number of auxiliary qubits. The average of the quantity F_f is plotted, which we define as the fidelity between the reduced density operator of the system qubits within the fragment (tracing out any auxiliary qubits a) and the reduced density of the same system qubits for the exact evolution of the full system (tracing out all environment qubits forming E). We use the generalization of fidelity for density matrices [66] to enable the focused evaluation of the fragment sub-system:

$$\begin{align} F_f & = \left( \textrm{Tr} \textrm{} \sqrt{\sqrt{\rho_f}\rho^{ex}_f\sqrt{\rho_f}} \right)^2 , \end{align} \tag{ 9 }$$

$$\begin{align} \rho_f & = \textrm{Tr}_a \ \rho_{f+a} , \end{align} \tag{ 10 }$$

$$\begin{align} \rho^{ex}_f & = \textrm{Tr}_E \ \rho^{ex} . \end{align} \tag{ 11 }$$

For a short evolution time, the scheme captures the correct state of the system qubits within the fragment. This time can be extended by the inclusion of additional auxiliary qubits.

Zoom In Zoom Out Reset image size

In figure 2(b), we consider a specific instance of the TFIM with J = 1.0 and h = 1.0. To test the scheme, we split the system into smaller partitions with $N_f = 3$. When we make use of two auxiliary qubits and mean-field corrections, the local expectation values exhibit little error for several units of Jt, as expected from the fidelity results.

In figure 3, we examine the scaling performance for the nearest neighbor model, increasing the number of fragments simulated with increasing N (keeping N_f constant) with $N_a = 2$. In the left panel, the first fragment is considered. This fragment contains the boundary of the chain, and consequently, the fragment's interface consists of only one missing interaction. The right panel considers the second fragment, which is on the interior of the chain for N > 6 and consequently neglects two interactions, leading to reduced performance. This is manifested in the reduced fragment fidelity going from N = 6 to N = 9 for fragment 2. However, there is no such visible drop going from N = 9 to N = 12 due to the monogamy of entanglement [67]—that is, although the number of qubits in the system grows, the qubits that are most strongly entangled with each other remain local to one fragment, and thus the amount of lost information shrinks as N is further increased. We therefore expect our classical scheme to scale well with N for systems that are locally interacting, and to serve as a strong approximation for moderate evolution times.

Zoom In Zoom Out Reset image size

Next, we examine the case of selective quantum information transfer between quantum simulators, applicable when non-local operations are available, even if only in a limited capacity. In this case, the fragmentation scheme can be modified to include limited quantum information transfer (a quantum channel [9]). The role of the auxiliary qubits shifts from being bystanders confined to a single fragment to qubits that are physically shared between simulators through selective non-local interactions, accomplished through qubit shuttling [50] or teleportation [51, 52] (see figure 1(c)). If the simulations are being executed in parallel on a single quantum simulator [42, 43], this would only require a few additional SWAP gates to include a limited number of cross-simulation interactions. In addition to providing more complete information transfer, a quantum channel further enables the active correction of auxiliary encoding as the system evolves. The selected number of auxiliary qubits places a limit on the number of qubits that are physically teleported / shuttled to a fragment; however, which environment qubits play this role can be changed from one time step to the next depending on which potential auxiliary qubit(s) have the largest contribution to the most recent estimate of the short-time error. When quantum channels and synchronized measurements are available, all correlation measurements are accessible. The quantity v(a) can thus be estimated for any a at any time. As the potential auxiliary qubits' contribution to the variance shift, new auxiliary qubits can be selected—that is, we can make a new selection for which qubit(s) physically interact with a fragment native to a different simulator. If the time steps are sufficiently small such that the first non-vanishing order in the error dominates, then this becomes optimal even for long simulation times.

To evaluate this scheme numerically, we abstract away the details of information transport; this topic has been investigated by other research in the context of distributed time evolution [68]. In our actively updated simulation, the quantity v(a) is calculated for each potential auxiliary qubit at each time step to determine its contribution to the short-time error. This requires the estimation of correlators between each potential auxiliary and each fragment system qubit (requiring $\mathcal{O}(N_f N_E N_{\alpha\beta}^2)$ per fragment, where N_αβ is the number of αβ interaction types), but if there is only one kind of interaction (as is the case for the TFIM and other Ising-like models, with $\alpha = \beta = z$), all relevant correlators can be estimated from a set of full system snapshot measurements. The largest contributors are selected to be auxiliaries—numerically, this amounts to keeping the interactions between these qubits and the fragment qubits, while zeroing the $J^{\alpha\beta}_{ij}$ coefficients of all other environment-fragment interactions (see appendix A.3). Any zeroed interactions can be approximately included via mean-field corrections. At the next time step, the selection of zeroed interactions might change due to a change in the selected auxiliaries for each fragment, as dictated by the short-time error.

Figure 4 compares this scheme (labeled 'Q. Channel' to indicate the addition of quantum information transfer) to the previous scheme in section 3.1, which involves only classical information transfer ('C. Channel'). The graph plots the fragment fidelity F_f averaged over 100 transverse field Ising-like models with h = 1.0 and randomly generated graphs J_ij . Each edge ij exists with probability 0.5, and edge weights J_ij are sampled from a Gaussian distribution with mean µ = 0.0 and width σ = 1.0. Furthermore, we randomly select a computational basis state to initialize the fragmented system. Although both schemes outperform the case of no information transfer (in red), the complicated long-range nature of the Hamiltonians considered challenges the previous scheme, which only employs classical information transfer. In contrast, the quantum scheme preserves a large fragment fidelity, even at late simulation times.

Zoom In Zoom Out Reset image size

The benefit of using short-time error to inform auxiliary selection can be isolated by evaluating the performance of each auxiliary choice independently. Consider a system of N = 12 qubits, fragmented into two groups of $N_f = 6$. This leaves six environment qubits from the perspective of each fragment that could be targeted by an auxiliary qubit. Selecting two auxiliary qubits ($N_a = 2$), we rank the six potential choices for target auxiliary encoding according to the size of v(a). In figure 5, the six target encoding choices are divided into three groups of two based on v(a), and each option is explored for randomly generated transverse field Ising-like Hamiltonians with h = 1.0, as considered in the previous section. A total of 100 such Hamiltonians are generated and simulated; the averaged results are presented in figure 5, where v₀ corresponds to encoding the two environment qubits with the largest value for v(a). On the left, the results are plotted for the case of classical information transfer. Any separation between the fidelity curves corresponding to different auxiliary choices indicates that the v(a) metric meaningfully separates the potential auxiliary choices according to fidelity performance. The fact that the ordering corresponds to the ranked choice is evidence that using short-time error to select auxiliary encoding propagates to better performance at later times. In red, we consider random auxiliary encoding. The random performance roughly converges to the middle-ranked choice v₁ and can be thought of as the performance averaged over auxiliary encoding. In the center, the results are plotted for the case of additional quantum information transfer, without actively updating the auxiliary encoding. The results qualitatively match those of the classical case, with slightly better performance overall, consistent with figure 4. In the right panel, we consider the quantum channel with actively updated auxiliary encoding. In this case, v₀ (v₂) corresponds to selecting the two auxiliary targets with the largest (smallest) values for v(a) at each decision. The performance of v₀ marginally increases with the introduction of active updates, while the performance of v₂ marginally decreases. However, the random performance increases most markedly. Here, the rapid shuffling of auxiliary qubit encoding allows the fragments to quickly share information, leading to performance comparable to the optimal variance choice, v₀. The random, actively updated case has the added advantage of being measurement-efficient as it forgoes any variance estimation, but the highly frequent change of auxiliary encoding may lead to an overhead in qubit routing / swapping in order to be realized.

Zoom In Zoom Out Reset image size

Finally, we note that in the averaged results presented in figure 5, the mean-field corrected simulation (plotted with a dashed line) outperforms the corresponding simulation that fully neglects these interface interactions for every case considered. Appendix B investigates the use of mean-field corrections to reduce simulation error through a numerical study.

We focus on pre-training brickwork circuits with a limited number of layers to constrain entanglement growth between fragments. We note that this is analogous to restricting ourselves to short evolution times to ensure more accurate results from a fragmented scheme. By utilizing a shallow circuit for pre-training, the error in the output of the ensemble of fragmented circuits is kept manageable. Although a circuit ansatz with high complexity is often necessary for interesting VQE applications in order to provide enough expressivity to reach the ground state [93, 94], fragmentation-based pre-training is still beneficial through the use of a layer-wise approach [95]. If a shallow brickwork circuit is placed ahead of a more expressive PQC ansatz, the brickwork layers can first be optimized using the fragmented approach. These layers serve to bring the state of the system to have some ground state overlap. The full circuit VQE can then be performed, initializing the leading brickwork layers of the circuit with the pre-trained parameter values and initializing the remaining gates of the ansatz to approximately act as identity—specifically, we choose to randomly initialize these parameters to be small values bounded by ±ε (with $\varepsilon = 10^{-5}$ for our results), to balance maintaining the optimized action of the initial layers after pre-training while avoiding training issues associated with a true identity initialization [73, 96]. The overall circuit layout is outlined in figure 7. Although we employ a fully-entangling architecture in the latter portion of the circuit to increase expressibility, we note that one might limit the number of gates connecting fragments, as is explored in [97], to efficiently create a fully distributed architecture.

Zoom In Zoom Out Reset image size

To evaluate our pre-training method, we focus on random Ising-like models. In section 5.2.1, we present the numerical performance of the scheme for the classical case of zero transverse field (h = 0). Having established the advantage of the approach, in section 5.2.2 we derive its success as stemming from the mean-field corrective terms included in the loss function, which shift the global minimum of the collective fragmented circuit to coincide with that of the full optimization problem. Finally, in section 5.2.3, we use perturbation theory and numerical simulation to demonstrate that our approach remains beneficial for $|h| \gt 0$, in the regime of a weak transverse field.

5.2.1. VQE results for MaxCut

We first benchmark the scheme using randomly generated classical Ising Hamiltonians, where the all-to-all J_ij interactions are sampled from a Gaussian distribution (mean µ = 0.0, width σ = 1.0) and the transverse field h is fixed to be zero. Finding the ground state of such models can be mapped to weighted MaxCut problems with positive and negative weights that are randomly generated [87]. The MaxCut problem is a graph partitioning problem that is known to be NP-hard, with applications ranging from network optimization to circuit design [98]. For the circuit ansatz, a fixed number of brickwork layers is used ($L_s = 4$) to keep this portion of the circuit shallow, while the all-to-all entangling portion of the circuit is made up of $L_c = N$ layers. The parameterized single qubit rotations within each layer are selected to be one rotation about x followed by one rotation about y, and the entangled gates are selected to be controlled z (CZ) rotations. All simulations are performed numerically using PennyLane [65], using the Adam optimizer built into PyTorch for all optimization [99, 100]. Lastly, note that parameter convergence (evaluated every 100 iterations) is used as the stopping criterion for both the fragmented circuit and full circuit optimization, with a maximum of 5000 iterations permitted.

To assess the performance of pre-training using circuit fragmentation, the same circuit is optimized using random initial values (referred to as 'vanilla VQE'). In this scenario, the parameters that correspond to the shallow, brickwork portion of the circuit are drawn from a uniform distribution ranging between $[0, 1)$. The all-to-all entangling portion of the gate is initialized near identity, as is done for the pre-training scheme. Figure 8 provides a case-by-case comparison between fragment-initialized VQE and vanilla VQE for 500 such models, for circuits of up to 15 qubits. In the top panels, the final percent error $\epsilon = (E - E_0) / |E_0|$ (where E₀ is the true ground state energy) is plotted for both approaches, along with the geometric mean of the results. The geometric mean of the fragment-initialized final error lies roughly three orders of magnitude below that of the vanilla VQE, with this gap growing even larger with increasing system size. For the larger system sizes, the vanilla VQE struggles to find a solution having $\epsilon \lt 10^{-2}$, while the fragment-initialized approach reaches $\epsilon \sim 10^{-7}$ for the same problem Hamiltonian. Moreover, using the same stopping criterion, the fragment-initialized VQE reaches this solution in fewer iterations ($N_\textrm{iter}$), decreasing the average number by nearly an order of magnitude, as illustrated by the bottom panels of figure 8. After successful pre-training, the parameters of the stitched-together circuit produce a loss that is already in the neighborhood of the minimum, so fewer iterations are required to reach convergence. We note that this corresponds to a significant reduction in the required execution time on quantum hardware, which can be prohibitively long to fully train a VQA [101]. For this simulation, we employ a batched optimization of T different fragmented circuits performed in parallel. In this approach, T different circuit partitions are randomly generated, and the resulting ensembles of fragmented circuits are optimized in parallel. The trained ensemble with lowest error is then used to initialize the full circuit. See appendix D.1 for an in-depth description of this approach and discussion of performance with batch size T.

Zoom In Zoom Out Reset image size

5.2.2. Solving maxcut with mean-field terms

5.2.3. Mean-field terms as first-order perturbation corrections

We now use perturbation theory to elucidate our technique of replacing multi-qubit interactions with mean fields when h ≠ 0. In the previous section, it is established that the ground state and ground state energy of an Ising-like Hamiltonian with zero transverse fields remain unchanged when one or more of the interactions are replaced by the corresponding mean-field approximation term. Following a similar argument, one can further establish that the computational basis states are stationary states of the mean-field corrected Hamiltonian $H_\textrm{MF}(|\psi\rangle)$, and therefore $H_\textrm{MF}(|\psi\rangle)$ and the unaltered Hamiltonian H share the same spectrum and set of eigenstates (although this term is used loosely for $H_\textrm{MF}(|\psi\rangle)$, as the dependence on $|\psi\rangle$ causes the stationary Schrödinger equation to deviate from a linear eigenvalue problem).

In this section, we consider adding a small transverse field to the classical Ising-like model, propelling the problem into the quantum domain. The first-order corrections to the ground state $|x^* \rangle$ and ground state energy E_g are computed using perturbation theory. The case of the mean-field corrected Hamiltonian $H_\textrm{MF}(|\psi\rangle)$ is treated with a version of perturbation theory modified to accommodate mean-field terms, and notably, the same first-order corrections to $|x^* \rangle$ and E_g are recovered. For a full derivation, please refer to appendix C.

Adding a transverse field, the unaltered Hamiltonian containing all interactions is given by:

$$\begin{align} H = H_{0} + H_\textrm{I} + \lambda V, \end{align} \tag{ 12 }$$

where H₀ contains the intra-fragment interactions:

$$\begin{align} H_{0} = -\sum_{\langle i, j \rangle \notin I} J_{ij} \hat{S}_{z}^{\left(i\right)} \hat{S}_{z}^{\left(j\right)} , \end{align} \tag{ 13 }$$

H_I contains the inter-fragment interactions:

$$\begin{align} H_\textrm{I} = -\sum_{\langle i, j \rangle \in I} J_{ij} \hat{S}_{z}^{\left(i\right)} \hat{S}_{z}^{\left(j\right)} , \end{align} \tag{ 14 }$$

V contains the perturbing transverse field:

$$\begin{align} V = -h \sum_{i} \hat{S}_{x}^{\left(i\right)} , \end{align} \tag{ 15 }$$

and λ is a perturbation parameter. We remind the reader that the inter-fragment interactions $H_\textrm{I}$ are those that will be replaced by mean-field corrections.

In contrast, the mean-field corrected Hamiltonian denoted $H_\textrm{MF}$ is given by:

$$\begin{align} H_\textrm{MF}\left(|\psi\rangle\right) = H_{0} + H_\textrm{I, MF}\left(|\psi\rangle\right) + \lambda V , \end{align} \tag{ 16 }$$

where the form of the Hamiltonian now depends on the state of the system due to the mean-field corrections:

$$\begin{align} H_\textrm{I, MF}\left(|\psi\rangle\right) = -\sum_{\langle i, j \rangle \in I} J_{ij} \hat{S}_{z}^{\left(i\right)} \langle \psi | \hat{S}_{z}^{\left(j\right)} | \psi \rangle . \end{align} \tag{ 17 }$$

Before any corrections can be computed, it is imperative to establish the correct zeroth order energies and eigenstates for each Hamiltonian. Following perturbation theory, the zeroth order eigenstates of H and $H_\textrm{MF}$ generally equal those of the unperturbed counterparts (that is, taking h = 0); these coincide with the set of computational basis states $\{|x\rangle \}$—including the unperturbed ground state, $|x^* \rangle$. However, there are degeneracies in the unperturbed Hamiltonians, and thus, degenerate perturbation theory is required.

When the unperturbed spectrum contains degeneracies, the proper linear combinations of the unperturbed eigenstates forming the degenerate subspace must be determined; these are the states that the perturbed eigenstates approach as $h \rightarrow 0$. The unperturbed Ising-like model possesses $\mathbb{Z}_2$ symmetry. Practically, this means that for each eigenstate $|x\rangle$, the 'flipped' eigenstate $|\bar{x}\rangle : = \bigotimes_i X_i |x_i\rangle$ is degenerate. For the unaltered Ising-like model H, the proper zeroth order eigenstates for the degenerate subspace containing the ground state are given by $|\pm_{x^{*}}\rangle = \frac{1}{\sqrt{2}} (|x^{*}\rangle \pm |\bar{x}^*\rangle)$. The transverse field will break the ground state degeneracy of H, and the positive superposition $|+_{x^{*}}\rangle$ is preferred by the ground state.

Shifting attention to the mean-field corrected Hamiltonian $H_\textrm{MF}$, the stationary Schrödinger equation is no longer linear in $|\psi\rangle$, and the linearity that characterizes quantum mechanics no longer applies. The notion of finding proper linear combinations is not an appropriate procedure due to the problem's nonlinearity. In particular, superpositions of degenerate eigenstates can yield different energies for $H_\textrm{MF}$ and thus effectively exist outside the degenerate subspace.

To illustrate this, consider a single mean-field factor, $\langle\hat{S}_{z}^{(i)}\rangle$, such as those within $H_\textrm{MF}$. While the expectation value of this quantity with respect to a computational basis state $|x^*\rangle$ yields

$$\begin{align} \langle x^* | \hat{S}_{z}^{\left(i\right)} | x^* \rangle = \frac{1}{2} \left(-1\right)^{x^*_{i}} , \end{align} \tag{ 18 }$$

evaluating the same term with respect to $|+_{x^*}\rangle$ leads to the term vanishing as

$$\begin{align} \begin{aligned} \langle +_{x^*} | \hat{S}_{z}^{\left(i\right)} | +_{x^*} \rangle & = \frac{1}{2} \left( \langle x^* | \hat{S}_{z}^{\left(i\right)} | x^* \rangle + \langle x^* | \hat{S}_{z}^{\left(i\right)} | \bar{x}^* \rangle + \langle \bar{x}^* | \hat{S}_{z}^{\left(i\right)} | x^* \rangle + \langle \bar{x}^* | \hat{S}_{z}^{\left(i\right)} | \bar{x}^* \rangle\right) \\ & = \frac{1}{4} \left( \left(-1\right)^{x^*_{i}} + \left(-1\right)^{\bar{x}^*_{i}} \right) \\ & = 0 . \end{aligned} \end{align} \tag{ 19 }$$

Notably for the ground state of the unperturbed Hamiltonian $|x^*\rangle$, this means that the pure computational basis states $|x^*\rangle, |\bar{x}^*\rangle$ are energetically preferred over any linear combination of them. Thus, for $H_\textrm{MF}$, the computational basis states remain the proper zeroth order eigenstates with a perturbative transverse field.

After establishing the zeroth order eigenstates and eigenenergies ($|k^{(0)}\rangle$ and $E_{k}^{(0)}$, respectively) of conventional Hamiltonians such as equation (12), perturbation theory proceeds by expanding $|k\rangle$ and E in λ in the stationary Schrödinger equation and equating orders of λ:

$$\begin{align} \begin{aligned} & \left( H_{0} + H_\textrm{I} + \lambda V \right) \left( |k^{\left(0\right)}\rangle + \lambda|k^{\left(1\right)}\rangle + \lambda^2 |k^{\left(2\right)}\rangle + \cdots \right) \\ & \quad = \left( E_{k}^{\left(0\right)} + \lambda E_{k}^{\left(1\right)} + \lambda^2 E_{k}^{\left(2\right)} + \cdots\right) \left( |k^{\left(0\right)}\rangle + \lambda|k^{\left(1\right)}\rangle + \lambda^2 |k^{\left(2\right)}\rangle + \cdots \right) . \end{aligned} \end{align} \tag{ 20 }$$

Following this procedure for H and carefully treating the degeneracy, the first-order energy correction $E_{k}^{(1)}$ vanishes and the first-order eigenstate correction takes the form:

$$\begin{align} | k^{\left(1\right)} \rangle = \sum_{m \notin D_k} \frac{\langle m^{\left(0\right)} |V| k^{\left(0\right)} \rangle}{E_{k}^{\left(0\right)} - E_{m}^{\left(0\right)}} |m^{\left(0\right)} \rangle , \end{align} \tag{ 21 }$$

where D_k represents the degenerate subspace that $|k^{(0)}\rangle$ occupies.

To derive the analogous correction to $H_\textrm{MF}$, we employ a modified approach to perturbation theory that can accommodate the nonlinearity of the stationary Schrödinger equation. In particular, the expanded form of $|k_\textrm{MF}\rangle$ is explicitly inserted into the state-dependant terms of $H_\textrm{MF}$ prior to equating orders of λ to compute the corrections. Following this procedure, the first order energy correction $E_{k, \textrm{MF}}^{(1)}$ again vanishes, and the first order eigenstate correction takes on an identical form to that of H:

$$\begin{align} | k_\textrm{MF}^{\left(1\right)} \rangle = \sum_{m \notin D_k} \frac{\langle m_\textrm{MF}^{\left(0\right)} |V| k_\textrm{MF}^{\left(0\right)} \rangle}{E_{k, \textrm{MF}}^{\left(0\right)} - E_{m, \textrm{MF}}^{\left(0\right)}} |m_\textrm{MF}^{\left(0\right)} \rangle . \end{align} \tag{ 22 }$$

There is one crucial difference between equations (21) and (22): the zeroth order eigenstates $|m^{(0)} \rangle$ and $|m_\textrm{MF}^{(0)} \rangle$. For the full Hamiltonian, each $|m^{(0)} \rangle$ has the form $|\pm_{y} \rangle \propto (|y\rangle \pm |\bar{y}\rangle)$, while each $|m_\textrm{MF}^{(0)} \rangle$ is a single computational basis state, $|y\rangle$. This leads to the fidelity between the first order ground state $|\psi_{g, \textrm{MF}} \rangle \propto |x^{*}\rangle + | x^{*(1)}_\textrm{MF}\rangle$ and that of the full Hamiltonian $|\psi_{g} \rangle \propto |+_{x^{*}}\rangle + |+_{x^{*}}^{(1)}\rangle$ to be F = 0.5 rather than perfect unity (see appendix C.3). Nonetheless, half overlap provides significant information about the true ground state for pre-training.

In figure 9, we examine the mean performance of the pre-training scheme as a function of the transverse field strength h. The case of neglecting mean-field terms during pre-training is also considered to highlight the vital role these corrections play at small values of h. When h = 0, the model is classical, and the global fragmented Hamiltonian with mean-field corrections shares the same ground state and ground state energy of the full model, as discussed in section 5.2.2. This leads to remarkable pre-training performance, even as N is increased. The fragment initialization that neglects mean-field terms is not guaranteed to share the ground state of the full model—in fact, the two are likely to be orthogonal. The pre-training will still feature low entanglement, which likely explains why the mean-field free initialization scheme outperforms random initialization, but overall, the average error exceeds that of the mean-field corrected case by orders of magnitude, particularly as N is increased. When a small transverse field is added to the model, the half overlap between the fragmented and full Hamiltonian ground states provided by including mean-field corrections leads to error orders of magnitude smaller than the other approaches. Only as h is further increased—entering the regime where first-order perturbation theory is inadequate to describe the ground state—do the approaches begin to perform comparably, with increasing error and required iterations.

Zoom In Zoom Out Reset image size

For added clarification, here we consider a toy model of N = 6 qubits fragmented into two groups such that $N_f = 3$ (see figure A1, a more explicitly labeled version of figure 1(b)). One auxiliary qubit is included for each fragment ($N_a = 1$). The auxiliary qubit of Fragment 1, a₁, targets qubit 5 in Fragment 2. Similarly, the auxiliary qubit of Fragment 2, a₂, targets qubit 1 in Fragment 1. In the figure, the interactions mediated by the auxiliary qubits are explicitly labeled to illustrate how each auxiliary qubit plays the role of one qubit from the environment. Which environment qubit is selected for this role will depend on the quantum variance $\text{var}( H - H_\textrm{I}^{(f)})$ (see section 2.2). The mean-field corrections applied to a₁ are included in the figure to account for the missing interactions with qubits 4 and 6 that the target qubit (qubit 5) would participate in if all qubits were present. Although only one such arrow appears in the figure, these corrections exist for each qubit participating in fewer interactions than it would in the full system, regardless of whether the qubit is categorized as a fragment qubit or an auxiliary targeting an environment qubit.

Zoom In Zoom Out Reset image size

Consider the error caused by the omission of interactions that form the interface of some particular fragment f and the remaining environment spins E. This Hamiltonian which we denote $H_\textrm{I}^{(f)}$ is subtly different than the previously introduced $H^{(f)}$, as it includes operators acting in the space of E in order to isolate the error caused by the section of the interaction interface produced by a particular fragment f:

$$\begin{align} \begin{aligned} H_\textrm{I}^{\left(\,f\,\right)} & = H^{\left(\,f\,\right)} + H^{\left(E\right)} \\ & = -\sum_{\langle i,j \rangle /\in I} \sum_{\alpha, \beta} J_{ij}^{\,\alpha, \beta} \hat{S}^{\left(i\right)}_{\alpha} \hat{S}^{\left(j\right)}_{\beta} - \sum_{i}^{N} h_{i} \hat{S}^{\left(i\right)}_x . \end{aligned} \end{align} \tag{ A.1 }$$

The fidelity of the state evolved under $H_\textrm{I}^{(f)}$ with that evolved under H is given by:

$$\begin{align} F\left(t\right) = |\langle \Psi| U^{\dagger}\left(t\right) U_\textrm{I}^{\left(\,f\,\right)}\left(t\right) |\Psi \rangle|^2 , \end{align} \tag{ A.2 }$$

where $U(t) = \exp(-iHt)$, $U_\textrm{I}^{(f)}(t) = \exp(-iH_\textrm{I}^{(f)}t)$, and $|\Psi_0\rangle$ is the current state of the system [66].

To analyze the fidelity metric, we can combine the product $U^{\dagger}(t) U_\textrm{I}^{(f)}(t)$ into a single exponential $\exp({Z})$. This final exponential argument is obtained using the BCH formula [61], which provides an expression for Z in terms of the nested commutators of the individual exponential arguments of $U^{\dagger}(t)$ and $U_\textrm{I}^{(f)}(t)$:

$$\begin{align} \begin{aligned} Z & = it\left(H - H_\textrm{I}^{\left(\,f\,\right)}\right) + \frac{1}{2}t^2\left[H, H_\textrm{I}^{\left(\,f\,\right)}\right] + \frac{1}{12}t^3 \left( i\left[H, \left[H, H_\textrm{I}^{\left(\,f\,\right)}\right]\right] + \left(-i\right)\left[H_\textrm{I}^{\left(\,f\,\right)}, \left[H_\textrm{I}^{\left(\,f\,\right)}, H\right]\right] \right) \\[.5ex] & \quad - \frac{1}{24}t^4 \left[H_{I}^{\left(\,f\,\right)}, \left[H, \left[H, H_\textrm{I}^{\left(\,f\,\right)}\right]\right]\right] + \cdots \end{aligned} \end{align} \tag{ A.3 }$$

Notice that each subsequent term in the BCH formula has a higher order in t. To approximate the error, we keep the terms up to small orders in t—specifically, we can write $Z = \sum_{n = 1}^{\infty} z_{n} t^{n}$, where the coefficients z_n do not depend on t and generally do not commute with each other, and truncate the sum after some n. We note that the first and second order z_n are given by:

$$\begin{align} z_1 & = i\left(H-H_\textrm{I}^{\left(\,f\,\right)}\right) , \end{align} \tag{ A.4 }$$

$$\begin{align} z_2 & = \frac{1}{2}\left[H, H_\textrm{I}^{\left(\,f\,\right)}\right] . \end{align} \tag{ A.5 }$$

Taylor expanding the exponential $\exp({Z})$ to second order in t:

$$\begin{align} \begin{aligned} U^{\dagger}\left(t\right) U_{I}^{\left(\,f\,\right)}\left(t\right) & = \exp\left( \sum_{n = 1}^{\infty} z_{n} t^{n} \right) \\[.5ex] & = 1 + \sum_{n = 1}^{\infty} z_{n} t^{n} + \frac{1}{2}\left( \sum_{n = 1}^{\infty} z_{n} t^{n} \right) ^2 + \cdots \\[.5ex] & = 1 + z_{1}t + t^2 \left(z_{2} + \frac{1}{2} z_{1}^2\right) + \mathcal{O}\left(t^3\right) . \end{aligned} \end{align} \tag{ A.6 }$$

The fidelity is then the modulus square of the expectation of the above expression, taken with respect to the initial state of the full system:

$$\begin{align} \begin{aligned} F\left(t\right) &\approx \big|1 + \langle z_{1} \rangle t + t^2 \left( \langle z_{2} \rangle + \frac{1}{2}\langle z_1 ^2 \rangle \right) \big|^2 \\[.5ex] & = 1 + t \left( \langle z_{1} \rangle + \langle z_{1} \rangle ^*\right) + t^2 \left( |\langle z_{1} \rangle|^2 + \langle z_{2} \rangle + \langle z_{2} \rangle^* + \frac{1}{2}\langle z_{1}^2 \rangle + \frac{1}{2}\langle z_{1}^2 \rangle^*\right) + \mathcal{O}\left(t^{3}\right) . \end{aligned} \end{align} \tag{ A.7 }$$

Plugging in the expressions for z₁ and z₂:

$$\begin{align} F\left(t\right) &\approx 1 + t\left(i \left\langle H - H_\textrm{I}^{\left(\,f\,\right)} \right\rangle -i \left\langle H - H_\textrm{I}^{\left(\,f\,\right)} \right\rangle ^* \right) + t^2 \left( |\left\langle H - H_\textrm{I}^{\left(\,f\,\right)} \right\rangle|^2 + \frac{1}{2} \left\langle \left[H, H_\textrm{I}^{\left(\,f\,\right)}\right] \right\rangle + \frac{1}{2} \left\langle \left[H, H_\textrm{I}^{\left(\,f\,\right)}\right] \right\rangle ^* \right. \nonumber\\[.5ex] & \left. \quad +\ \left(i\right)^2 \frac{1}{2}\left\langle \left(H - H_\textrm{I}^{\left(\,f\,\right)}\right)^2 \right\rangle + \left(-i\right)^2 \frac{1}{2} \left\langle \left(H - H_\textrm{I}^{\left(\,f\,\right)}\right)^2 \right\rangle ^*\right) . \end{align} \tag{ A.8 }$$

Notice that, as both H and $H_\textrm{I}^{(f)}$ are Hermitian, $\langle H - H_\textrm{I}^{(f)} \rangle ^* = \langle H \rangle ^* - \langle H_\textrm{I}^{(f)} \rangle ^* = \langle H \rangle - \langle H_\textrm{I}^{(f)} \rangle$. Thus, the terms first order in t cancel with one another. The second-order terms can be simplified by expanding them:

$$\begin{align} \begin{aligned} F\left(t\right) &\approx 1 + t^2 \left( |\left\langle H - H_\textrm{I}^{\left(\,f\,\right)} \right\rangle|^2 + \frac{1}{2} \left(\left\langle H H_\textrm{I}^{\left(\,f\,\right)}\right\rangle - \left\langle H_\textrm{I}^{\left(\,f\,\right)} H\right\rangle\right) + \frac{1}{2} \left(\left\langle H H_\textrm{I}^{\left(\,f\,\right)}\right\rangle^* - \left\langle H_\textrm{I}^{\left(\,f\,\right)} H\right\rangle^*\right)\right. \\[.5ex] & \left. \quad - \frac{1}{2} \left(\langle H^2 \rangle - \left\langle H H_\textrm{I}^{\left(\,f\,\right)} \right\rangle - \langle H_{I_f} H\rangle + \left\langle \left(H_\textrm{I}^{\left(\,f\,\right)}\right)^{2} \right\rangle\right) \right. \\[.5ex] & \left. \quad -\ \frac{1}{2} \left(\langle H^2 \rangle^* - \left\langle H H_\textrm{I}^{\left(\,f\,\right)} \right\rangle^* - \left\langle H_\textrm{I}^{\left(\,f\,\right)} H \right\rangle^* + \left\langle \left(H_\textrm{I}^{\left(\,f\,\right)}\right)^{2} \right\rangle^*\right)\right) . \end{aligned} \end{align} \tag{ A.9 }$$

Using the fact that $\langle A B \rangle^* = \langle B^{\dagger} A^{\dagger} \rangle$ to simplify, we arrive at a compact expression for the fidelity to second order in t:

$$\begin{align} \begin{aligned} F\left(t\right) & = 1 + t^2 \left( |\left\langle H - H_\textrm{I}^{\left(\,f\,\right)} \right\rangle|^2 - \left\langle \left(H - H_\textrm{I}^{\left(\,f\,\right)}\right)^2 \right\rangle \right) + \mathcal{O}\left(t^3\right) \\[.5ex]&\approx 1 - \text{var}\left( H - H_\textrm{I}^{\left(\,f\,\right)}\right) t^2 . \end{aligned} \end{align} \tag{ A.10 }$$

The error $\epsilon(t) = 1 - F(t)$ is thus given by $\text{var}( H - H_\textrm{I}^{(f)}) t^2$ for short times t.

When simulators are linked via a quantum channel, the states of the fragments no longer live in separate Hilbert spaces; they comprise the state of the collective system, which lives in the larger Hilbert space of size 2^N . To simulate this numerically, we evolve a modified Hamiltonian that lives in the full Hilbert space. This Hamiltonian involves adjusted connectivity matrices $J_{ij}^{\,\alpha\beta}$ to match the connectivity provided by the current auxiliary encoding. That is, defining a set S_f for each fragment containing the list of fragment qubits and target qubits for auxiliaries of fragment f:

$$\begin{align} S_f = \left\{i \in f\,\right\} \cup \left\{a_f\,\right\} , \end{align} \tag{ A.11 }$$

we implement the union of connectivity matrices $\cup_f J_{ij}^{\,\alpha\beta, (f)}$ for each αβ interaction type, defined by:

$$\begin{align} J_{ij}^{\,\alpha\beta, \left(\,f\,\right)} = \begin{cases} J_{ij}^{\,\alpha\beta} & \text{if } i, j \in S_f\\ 0 & \text{otherwise} . \end{cases} \end{align} \tag{ A.12 }$$

Any interactions still zeroed in the union $\cup_f J_{ij}^{\,\alpha\beta, (f)}$ are approximately included via mean-field corrections. If the auxiliary encoding is updated, the union $\cup_f J_{ij}^{\,\alpha\beta, (f)}$ will change due to the change in selected auxiliaries $\{a_f\}$ of each fragment.

In this section, we report results for the time simulation of the Heisenberg model. We consider a 1D Heisenberg model with uniform coupling J, governed by the following Hamiltonian:

$$\begin{align} H_\textrm{Heis} = - J \sum_{\langle i,j\rangle \in n.n.}\left( \hat{S}_{x}^{\left(i\right)} \hat{S}_{x}^{\left(j\right)} + \hat{S}_{y}^{\left(i\right)} \hat{S}_{y}^{\left(j\right)} + \hat{S}_{z}^{\left(i\right)} \hat{S}_{z}^{\left(j\right)} \right). \end{align} \tag{ A.13 }$$

Note that this model features significantly larger interaction as compared to the Ising model previously considered. To adjust for this increased interaction, we modify the effective simulation time to be $\sqrt{3J^2}t$. Figure A2, presents the fragment fidelity for the fragmented evolution of the Heisenberg model. Repeating the setup of figure 2(a), we split a system of 12 qubits into 2 fragments and evolve the system with different numbers of auxiliary qubits, both with and without mean-field corrections. The results provided are averaged over 100 different initial states, sampled from computational basis states. The general trend again matches that of figure 2(a), with progressively better performance for larger N_a . However, we note that the mean-field corrections provide less benefit due to the increased interaction likely leading to a more entangled structure that cannot be mimicked by mean-field corrections. This concept is explored further in appendix B.

Zoom In Zoom Out Reset image size

In this section, we follow the degenerate perturbation theory procedure outlined in [56]. When all interactions are included, the unperturbed Hamiltonian is given by the sum of H₀ and $H_\textrm{I}$. The eigendecomposition of this operator defines the zeroth order eigenenergies and eigenstates:

$$\begin{align} H_{0} + H_\textrm{I} = \sum_{k} E_{k}^{\left(0\right)} |k^{\left(0\right)}\rangle \langle k^{\left(0\right)} | . \end{align} \tag{ C.7 }$$

The unperturbed Hamiltonian is classical, with computational basis states for eigenstates. It also possesses $\mathbb{Z}_2$ symmetry, such that each eigenstate $|x\rangle$ and its 'flipped' version $|\bar{x}\rangle$ are degenerate, including the ground state $|x^{*}\rangle$. It is necessary to find the correct linear combinations of the degenerate states to properly calculate the higher order corrections. Typically, this is accomplished by diagonalizing the perturbing Hamiltonian V in the subspace of degenerate eigenstates, which yields the proper zeroth order eigenstates as well as the first order energy corrections, $\{E_{k}^{(1)}\}$. The particular V of this problem—a global transverse field—vanishes in the subspace of $|x\rangle$, $|\bar{x}\rangle$, requiring the first order energy corrections $\{E_{k}^{(1)}\}$ to vanish but providing no insight into the correct zeroth order eigenstates.

To calculate the first order eigenstate corrections, the correction is split into one within the degenerate space D_k and one in the space outside D_k , which we define as $\bar{D}_k$. For the latter correction, we have:

$$\begin{align} P_{\bar{D}_{k}} | k^{\left(1\right)} \rangle = \sum_{m \notin D_{k}} \frac{\langle m^{\left(0\right)} |V| k^{\left(0\right)} \rangle}{E_{k}^{\left(0\right)} - E_{m}^{\left(0\right)}} |m^{\left(0\right)} \rangle , \end{align} \tag{ C.8 }$$

where $P_{\bar{D}_{k}}$ is the projects out the degenerate subspace. Following [56], the correction within the degenerate subspace is given by:

$$\begin{align} P_{D_{k}} | k^{\left(1\right)} \rangle = \sum_{\substack{k^{^{\prime}} \in D_{k}, \\ k^{^{\prime}} \neq k}} \frac{P_{D_{k}}|k^{^{\prime}}{\left(0\right)}\rangle}{E_{k}^{\left(1\right)} - E_{k^{^{\prime}}}^{\left(1\right)}} \langle k^{^{\prime}}{\left(0\right)} | V P_{\bar{D}_k} \frac{1}{E_{k}^{\left(0\right)} - H_0 - H_\textrm{I}} P_{\bar{D}_k} V |k^{\left(0\right)} \rangle . \end{align} \tag{ C.9 }$$

The first order energy difference in the denominator of the first fraction of this expression is singular, because the degeneracy within D_k has not been lifted by the first order energy correction. However, this apparent issue provides us with the correct zeroth order eigenstates required to cancel the singular denominator: $\{|k^{(0)}\rangle\}$ must be selected to diagonalize the object $W : = V P_{\bar{D}_k} (E_{k}^{(0)} - H_0 - H_\textrm{I})^{-1} P_{\bar{D}_k} V$ within the subspace of degenerate states, such that the quantity $\langle k^{^{\prime}}{(0)} | W |k^{(0)} \rangle$ vanishes for all $k^{^{\prime}} \neq k$ within D_k . With this choice for $\{|k^{(0)}\rangle\}$, the first order eigenstate correction within the degenerate space vanishes, and the full first order eigenstate correction takes the form:

$$\begin{align} | k^{\left(1\right)} \rangle = \sum_{m \notin D_{k}} \frac{\langle m^{\left(0\right)} |V| k^{\left(0\right)} \rangle}{E_{k}^{\left(0\right)} - E_{m}^{\left(0\right)}} |m^{\left(0\right)} \rangle . \end{align} \tag{ C.10 }$$

One can show that $\langle x | W | x \rangle = \langle x | W | \bar{x} \rangle = \langle \bar{x} | W | x \rangle = \langle \bar{x} | W | \bar{x} \rangle$ using the $\mathbb{Z}_2$ symmetry of the unperturbed Hamiltonian. Thus, the correct zeroth order eigenstates are the symmetric and antisymmetric superpositions of $|x\rangle$ and $|\bar{x}\rangle$, given by $|\pm_{x}\rangle = (1/\sqrt{2})(|x\rangle \pm |\bar{x}\rangle)$.

Throughout this derivation, we will take advantage of the particular form of $H_\textrm{I, MF}(|\psi\rangle)$ to make simplifications.

Consider a modified Schrödinger equation that takes into account the state-dependence of $H_\textrm{MF}(|\psi\rangle)$:

$$\begin{align} \left(H_{0} + H_\textrm{I, MF}\left(|k_\textrm{MF}\rangle\right) + \lambda V\right) |k_\textrm{MF}\rangle = E_{k,\textrm{MF}}|k_\textrm{MF}\rangle. \end{align} \tag{ C.11 }$$

In conventional perturbation theory, the eigenstates and eigenvalues ($|k_\textrm{MF}\rangle$ and $E_{k, \textrm{MF}}$, respectively) are expanded about their unperturbed counterparts, $|k_\textrm{MF}^{(0)}\rangle$ and $E_{k, \textrm{MF}}^{(0)}$. However, the modified Schrödinger equation written above is no longer an eigenvalue problem; the dependence on $|k_\textrm{MF}\rangle$ is non-linear. In the spirit of conventional perturbation theory, we will proceed in the usual manner, taking special care to include the non-linearity. Expanding $|k_\textrm{MF}\rangle$ and $E_{k, \textrm{MF}}$ in orders of λ:

$$\begin{align} |k_\textrm{MF}\rangle = |k_\textrm{MF}^{\left(0\right)}\rangle + \lambda|k_\textrm{MF}^{\left(1\right)}\rangle + \lambda^2 |k_\textrm{MF}^{\left(2\right)}\rangle + \cdots , \end{align} \tag{ C.12 }$$

$$\begin{align} E_{k, \textrm{MF}} = E_{k, \textrm{MF}}^{\left(0\right)} + \lambda E_{k,\textrm{MF}}^{\left(1\right)} + \lambda^2 E_{k, \textrm{MF}}^{\left(2\right)} + \cdots . \end{align} \tag{ C.13 }$$

Replacing $|k_\textrm{MF}\rangle$ and $E_{k, \textrm{MF}}$ in the modified Schrödinger equation:

$$\begin{align*} &\left( H_{0} + H_\textrm{I, MF}\left(|k_\textrm{MF}^{\left(0\right)}\rangle + \lambda|k_\textrm{MF}^{\left(1\right)}\rangle + \lambda^2 |k_\textrm{MF}^{\left(2\right)}\rangle + \cdots\right) + \lambda V \right) \left( |k_\textrm{MF}^{\left(0\right)}\rangle + \lambda|k_\textrm{MF}^{\left(1\right)}\rangle + \lambda^2 |k_\textrm{MF}^{\left(2\right)}\rangle + \cdots \right) \\ & \quad = \left( E_{k, \textrm{MF}}^{\left(0\right)} + \lambda E_{k, \textrm{MF}}^{\left(1\right)} + \lambda^2 E_{k, \textrm{MF}}^{\left(2\right)} + \cdots\right) \left( |k_\textrm{MF}^{\left(0\right)}\rangle + \lambda|k_\textrm{MF}^{\left(1\right)}\rangle + \lambda^2 |k_\textrm{MF}^{\left(2\right)}\rangle + \cdots \right). \end{align*}$$

Now, we can isolate and equate the various orders of λ to calculate the perturbations. If we first consider the zeroth order, we recover the unperturbed spectrum:

$$\begin{align} \left(H_{0} + H_\textrm{I, MF}\left(|k_\textrm{MF}^{\left(0\right)}\rangle\right)\right) |k_\textrm{MF}^{\left(0\right)}\rangle = E_{k, \textrm{MF}}^{\left(0\right)} |k_\textrm{MF}^{\left(0\right)}\rangle. \end{align} \tag{ C.14 }$$

Moving to higher orders, it is helpful to write the explicit form of $H_\textrm{I, MF}(|\psi\rangle)$ to properly count the orders of λ. The first order terms take the form:

$$\begin{align} \begin{aligned} & H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle + V |k_\textrm{MF}^{\left(0\right)}\rangle -\sum_{\langle i, j \rangle \in I} J_{ij} \left( \sigma_{z}^{\left(i\right)} \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle \right) |k_\textrm{MF}^{\left(1\right)}\rangle \\ & \quad -\sum_{\langle i, j \rangle \in I} J_{ij} \left( \sigma_{z}^{\left(i\right)} \langle k_\textrm{MF}^{\left(1\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle + \sigma_{z}^{\left(i\right)} \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(1\right)} \rangle \right) |k_\textrm{MF}^{\left(0\right)}\rangle \\ & \qquad = E_{k, \textrm{MF}}^{\left(0\right)} |k_\textrm{MF}^{\left(1\right)}\rangle + E_{k, \textrm{MF}}^{\left(1\right)} |k_\textrm{MF}^{\left(0\right)}\rangle . \end{aligned} \end{align} \tag{ C.15 }$$

The operators $\sigma_{z}^{(j)}$ are diagonal in the computational basis $\{|k_\textrm{MF}^{(0)}\rangle\}$. By definition, corrections to $|k_\textrm{MF}^{(0)}\rangle$ such as $|k_\textrm{MF}^{(1)}\rangle$ will be orthogonal to $|k_\textrm{MF}^{(0)}\rangle$; therefore, terms of the form $\langle k_\textrm{MF}^{(1)} | \sigma_{z}^{(j)} | k_\textrm{MF}^{(0)} \rangle$ will vanish. The first order equation thus simplifies to:

$$\begin{align} H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle + V |k_\textrm{MF}^{\left(0\right)}\rangle -\sum_{\langle i, j \rangle \in I} J_{ij} \left( \sigma_{z}^{\left(i\right)} \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle \right) |k_\textrm{MF}^{\left(1\right)}\rangle = E_{k, \textrm{MF}}^{\left(0\right)} |k_\textrm{MF}^{\left(1\right)}\rangle + E_{k, \textrm{MF}}^{\left(1\right)} |k_\textrm{MF}^{\left(0\right)}\rangle . \end{align} \tag{ C.16 }$$

To determine the first order energy shift, we project the above equation with the unperturbed eigenstate $\langle k_\textrm{MF}^{(0)}|$:

$$\begin{align} \begin{aligned} & \langle k_\textrm{MF}^{\left(0\right)}| H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle + \langle k_\textrm{MF}^{\left(0\right)}|V |k_\textrm{MF}^{\left(0\right)}\rangle -\sum_{\langle i, j \rangle \in I} J_{ij} \langle k_\textrm{MF}^{\left(0\right)}|\sigma_{z}^{\left(i\right)} |k_\textrm{MF}^{\left(1\right)}\rangle \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle \\ & \quad = E_{k, \textrm{MF}}^{\left(0\right)} \langle k_\textrm{MF}^{\left(0\right)}|k_\textrm{MF}^{\left(1\right)}\rangle + E_{k, \textrm{MF}}^{\left(1\right)} \langle k_\textrm{MF}^{\left(0\right)}|k_\textrm{MF}^{\left(0\right)}\rangle . \end{aligned} \end{align} \tag{ C.17 }$$

The terms originating from $H_\textrm{I, MF}(|\psi\rangle)$ as well as the first term on the right-hand side of the equation will vanish due to the previously discussed fact that $|k_\textrm{MF}^{(1)}\rangle$ will be orthogonal to $|k_\textrm{MF}^{(0)}\rangle$. Let us examine the first term on the left-hand side, making use of the definition of the unperturbed spectrum:

$$\begin{align*} \langle k_\textrm{MF}^{\left(0\right)}| H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle & = \langle k_\textrm{MF}^{\left(0\right)}| E_{k, \textrm{MF}}^{\left(0\right)} - H_\textrm{I, MF}\left(|k_\textrm{MF}^{0}\rangle\right) |k_\textrm{MF}^{\left(1\right)}\rangle \\ & = E_{k, \textrm{MF}}^{\left(0\right)} \langle k_\textrm{MF}^{\left(0\right)} | k_\textrm{MF}^{\left(1\right)}\rangle + \sum_{\langle i, j \rangle \in I} J_{ij} \langle k_\textrm{MF}^{\left(0\right)}|\sigma_{z}^{\left(i\right)} |k_\textrm{MF}^{\left(1\right)}\rangle \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle \\ & = 0 . \end{align*}$$

Thus, this term vanishes for the same reason, and the first order energy shift from conventional perturbation theory is recovered:

$$\begin{align} E_{k, \textrm{MF}}^{\left(1\right)} = \langle k_\textrm{MF}^{\left(0\right)}|V |k_\textrm{MF}^{\left(0\right)}\rangle . \end{align} \tag{ C.18 }$$

Of course, this energy correction vanishes for all k, in the same way that the first order energy shift of the full Hamiltonian vanishes.

To determine the first order eigenstate shift $|k_\textrm{MF}^{(1)}\rangle$, we project the first order equation given in equation (C.15) with the unperturbed eigenstate $\langle m_\textrm{MF}^{(0)}|$, where m ≠ k:

$$\begin{align} \begin{aligned} & \langle m_\textrm{MF}^{\left(0\right)}| H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle + \langle m_\textrm{MF}^{\left(0\right)}|V |k_\textrm{MF}^{\left(0\right)}\rangle -\sum_{\langle i, j \rangle \in I} J_{ij} \langle m_\textrm{MF}^{\left(0\right)}|\sigma_{z}^{\left(i\right)} |k_\textrm{MF}^{\left(1\right)}\rangle \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle \\ & \quad = E_{k, \textrm{MF}}^{\left(0\right)} \langle m_\textrm{MF}^{\left(0\right)}|k_\textrm{MF}^{\left(1\right)}\rangle + E_{k, \textrm{MF}}^{\left(1\right)} \langle m_\textrm{MF}^{\left(0\right)}|k_\textrm{MF}^{\left(0\right)}\rangle . \end{aligned} \end{align} \tag{ C.19 }$$

The only term guaranteed to vanish is the second term on the right-hand side, as $\langle m_\textrm{MF}^{(0)}|k_\textrm{MF}^{(0)}\rangle = 0$ for m ≠ k. The term $\langle m_\textrm{MF}^{(0)}| H_{0} |k_\textrm{MF}^{(1)}\rangle$ can also be simplified, with careful attention to the definition of the zeroth order energy $E_{m, \textrm{MF}}^{(0)}$:

$$\begin{align} \begin{aligned} \langle m_\textrm{MF}^{\left(0\right)}| H_{0} |k_\textrm{MF}^{\left(1\right)}\rangle & = \langle m_\textrm{MF}^{\left(0\right)}| E_{m, \textrm{MF}}^{\left(0\right)} - H_\textrm{I, MF}\left(|k_\textrm{MF}^{0}\rangle\right) |k_\textrm{MF}^{\left(1\right)}\rangle \\ & = E_{m, \textrm{MF}}^{\left(0\right)} \langle m_\textrm{MF}^{\left(0\right)} | k_\textrm{MF}^{\left(1\right)}\rangle + \sum_{\langle i, j \rangle \in I} J_{ij} \langle m_\textrm{MF}^{\left(0\right)}|\sigma_{z}^{\left(i\right)} |k_\textrm{MF}^{\left(1\right)}\rangle \langle k_\textrm{MF}^{\left(0\right)} | \sigma_{z}^{\left(j\right)} | k_\textrm{MF}^{\left(0\right)} \rangle . \end{aligned} \end{align} \tag{ C.20 }$$

Notice that the mean-field terms in the expression above precisely cancel with the mean-field terms remaining in equation (C.19). Simplifying:

$$\begin{align} \langle m_\textrm{MF}^{\left(0\right)}|V |k_\textrm{MF}^{\left(0\right)}\rangle = \left(E_{k, \textrm{MF}}^{\left(0\right)} - E_{m, \textrm{MF}}^{\left(0\right)}\right) \langle m^{\left(0\right)}|k^{\left(1\right)}\rangle . \end{align} \tag{ C.21 }$$

Therefore, the first order correction to the eigenstate is identical to that of the full Hamiltonian with no mean-field corrections:

$$\begin{align} \begin{aligned} |k_\textrm{MF}^{\left(1\right)}\rangle & = \sum_{m} |m_\textrm{MF}^{\left(0\right)}\rangle \langle m_\textrm{MF}^{\left(0\right)} |k_\textrm{MF}^{\left(1\right)}\rangle \\ & = \sum_{m \notin D_{k}} |m_\textrm{MF}^{\left(0\right)}\rangle \frac{\langle m_\textrm{MF}^{\left(0\right)}|V |k_\textrm{MF}^{\left(0\right)}\rangle}{E_{k, \textrm{MF}}^{\left(0\right)} - E_{m, \textrm{MF}}^{\left(0\right)}} . \end{aligned} \end{align} \tag{ C.22 }$$

In this section, we derive the overlap between an eigenstate of an Ising-like Hamiltonian H with a small transverse field and the corresponding eigenstate of a mean-field corrected version of the Hamiltonian, $H_\textrm{MF}$, to leading order in perturbation theory. This overlap is approximately correct as long as the strength of the transverse field $|h|$ is small enough that first order perturbation theory is a good description of the eigenstates.

To first order, an eigenstate of each Hamiltonian is given by:

$$\begin{align} |\psi_{k}\rangle & = \mathcal{N}\left( |k^{\left(0\right)} \rangle + |k^{\left(1\right)} \rangle \right) , \end{align} \tag{ C.23 }$$

$$\begin{align} |\psi_{k,\textrm{MF}}\rangle & = \mathcal{N}_\textrm{MF}\left( |k_\textrm{MF}^{\left(0\right)} \rangle + |k_\textrm{MF}^{\left(1\right)} \rangle \right) , \end{align} \tag{ C.24 }$$

where $\mathcal{N}$, $\mathcal{N}_\textrm{MF}$ are normalization factors. The zeroth order eigenstates are normalized by definition, so re-normalization is required when any higher order corrections are included.

The overlap between the two states is given by their inner product:

$$\begin{align} \langle \psi_{k} | \psi_{k, \textrm{MF}} \rangle = \mathcal{N}^{*}\mathcal{N}_\textrm{MF} \left( \langle k^{\left(0\right)} | k_\textrm{MF}^{\left(0\right)} \rangle + \langle k^{\left(1\right)} | k_\textrm{MF}^{\left(1\right)} \rangle \right) , \end{align} \tag{ C.25 }$$

where the cross-terms $\langle k^{(1)} | k_\textrm{MF}^{(0)} \rangle$, $\langle k^{(0)} | k_\textrm{MF}^{(1)} \rangle$ have been dropped due to the fact that the first order corrections $| k^{(1)} \rangle$, $| k_\textrm{MF}^{(1)} \rangle$ lie outside the degenerate subspace of their corresponding zeroth order states $| k^{(0)} \rangle$ and $| k_\textrm{MF}^{(0)} \rangle$, which both live in the same degenerate subspace D_k .

As argued in previous sections, the zeroth order eigenstates $|k_\textrm{MF}^{(0)}\rangle$ are computational basis states $|x\rangle$, while the zeroth order eigenstates $|k^{(0)}\rangle$ are superpositions of two computational basis states related by $\mathbb{Z}_2$ symmetry, $|\pm_{x}\rangle = (1/\sqrt{2})(|x\rangle \pm |\bar{x}\rangle)$. Without loss of generality, we will assume $| k^{(0)} \rangle$ is the positive superposition $|+_{x}\rangle$, as is the case for the ground state given our conventions. The overlap between zeroth order eigenstates follows directly from these definitions: $\langle k^{(0)} | k_\textrm{MF}^{(0)} \rangle = \frac{1}{\sqrt{2}}$.

To find the overlap between first order corrections, we rewrite these corrections as a sum over degenerate subspaces D_m rather than over zeroth order eigenstates. This is possible because of the $\mathbb{Z}_2$ symmetry of each Hamiltonian, allowing the full set of eigenstates to be split into degenerate pairs.

First, we consider $| k^{(1)} \rangle$. The correction can be expressed abstractly as a weighted sum over states within each degenerate subspace D_m :

$$\begin{align} \begin{aligned} |k^{\left(1\right)}\rangle & = \sum_{D_{m} \neq D_{k}} c\left(m, k\right) |D_{m} \rangle \\ & = \sum_{D_{y} \neq D_{x}} c\left(y, x\right) |D_{y} \rangle . \end{aligned} \end{align} \tag{ C.26 }$$

In the second line, the labels are changed to explicitly highlight the relationship to computational basis states $x, y$. Using the previous definition in equation (21):

$$\begin{align} \begin{aligned} c\left(y, x\right) |D_{y} \rangle & = \frac{1}{\Delta E_{x, y}} \left( \langle +_{y} | V |+_{x}\rangle |+_{y}\rangle + \langle -_{y} | V |+_{x}\rangle |-_{y}\rangle \right) \\ & = \frac{1}{\Delta E_{x, y}} \left[ \frac{1}{2} \left( V_{yx} + V_{y\bar{x}} + V_{\bar{y}x} + V_{\bar{y}\bar{x}}\right) |+_{y}\rangle + \frac{1}{2} \left( V_{yx} + V_{y\bar{x}} - V_{\bar{y}x} - V_{\bar{y}\bar{x}}\right) |-_{y}\rangle \right] \\ & = \frac{1}{\Delta E_{x, y}} \left( V_{yx} + V_{y\bar{x}} \right) |+_{y}\rangle , \end{aligned} \end{align} \tag{ C.27 }$$

where $\Delta E_{x, y}$ is the energy difference between D_x and D_y , and $V_{xy} = \langle x | V | y \rangle$ is an element of the matrix V in the computational basis. We have made use of properties of V (namely, $V_{yx} = V_{\bar{y}\bar{x}}$ and $V_{y\bar{x}} = V_{\bar{y}x}$) to simplify.

Rewriting $| k_\textrm{MF}^{(1)} \rangle$:

$$\begin{align} \begin{aligned} |k_\textrm{MF}^{\left(1\right)}\rangle & = \sum_{D_{m} \neq D_{k}} c_\textrm{MF}\left(m, k\right) |D_{m, \textrm{MF}} \rangle \\ & = \sum_{D_{y} \neq D_{x}} c_\textrm{MF}\left(y, x\right) |D_{y, \textrm{MF}} \rangle . \end{aligned} \end{align} \tag{ C.28 }$$

The differing zeroth order states leads to a different expression for the weighted sum:

$$\begin{align} \begin{aligned} c_\textrm{MF}\left(y, x\right) |D_{y, \textrm{MF}} \rangle & = \frac{1}{\Delta E_{x, y}} \left( \langle y | V | x \rangle |y \rangle + \langle \bar{y} | V | x \rangle |\bar{y}\rangle \right) \\ & = \frac{1}{\Delta E_{x, y}} \left( V_{yx} |y \rangle + V_{\bar{y}x} |\bar{y}\rangle \right) . \end{aligned} \end{align} \tag{ C.29 }$$

The inner product between first order corrections is now straightforward:

$$\begin{align} \begin{aligned} \langle k^{\left(1\right)} | k_\textrm{MF}^{\left(1\right)} \rangle & = \sum_{D_{y} \neq D_{x}} c^{*}\left(y, x\right) c_\textrm{MF}\left(y, x\right) \langle D_y | D_{y,\textrm{MF}} \rangle \\ & = \sum_{D_{y} \neq D_{x}} \frac{1}{\left(\Delta E_{x,y}\right)^2} \langle +_{y}| \left( V_{yx} + V_{y\bar{x}} \right) \left( V_{yx} |y \rangle + V_{\bar{y}x} |\bar{y}\rangle \right) \\ & = \frac{1}{\sqrt{2}} \sum_{D_{y} \neq D_{x}} \frac{1}{\left(\Delta E_{x,y}\right)^2} \left( V_{yx} + V_{y\bar{x}} \right)^2 \\ & = \frac{1}{\sqrt{2}} \langle k^{\left(1\right)} | k^{\left(1\right)} \rangle . \end{aligned} \end{align} \tag{ C.30 }$$

Finally, the overlap between eigenstates to first order is given by:

$$\begin{align} \begin{aligned} \langle \psi_{k} | \psi_{k, \textrm{MF}} \rangle & = \frac{1}{\sqrt{2}} \mathcal{N}^{*}\mathcal{N}_\textrm{MF} \left( \langle k^{\left(0\right)} | k^{\left(0\right)} \rangle + \langle k^{\left(1\right)} | k^{\left(1\right)} \rangle \right)\\ & = \frac{1}{\sqrt{2}} \frac{\mathcal{N}_\textrm{MF}}{\mathcal{N}} . \end{aligned} \end{align} \tag{ C.31 }$$

Consider the form of the normalization coefficient $\mathcal{N}$ and $\mathcal{N}_\textrm{MF}$:

$$\begin{align} \mathcal{N} & = \sqrt{\frac{1}{1 + \langle k^{\left(1\right)} | k^{\left(1\right)} \rangle}} \nonumber\\ & = \sqrt{\frac{1}{1 + \sum_{D_{y} \neq D_{x}} \left( V_{yx} + V_{y\bar{x}} \right)^2 / \left(\Delta E_{x,y}\right)^2}} \end{align} \tag{ C.32 }$$

$$\begin{align} \mathcal{N}_\textrm{MF} & = \sqrt{\frac{1}{1 + \langle k_\textrm{MF}^{\left(1\right)} | k_\textrm{MF}^{\left(1\right)} \rangle}} \nonumber \\ & = \sqrt{\frac{1}{1 + \sum_{D_{y} \neq D_{x}} \left( |V_{yx}|^2 + |V_{y\bar{x}}|^2 \right) / \left(\Delta E_{x,y}\right)^2}}. \end{align} \tag{ C.33 }$$

In the computational basis, the matrix elements of the perturbing transverse field V are real and all have the same sign (with the sign depending on the sign of h). Thus, any product of two matrix elements will be a positive number or zero. Likewise, the squared energy difference $(\Delta E_{x,y})^2$ is positive. The denominator of $\mathcal{N}$ contains two additional factors from the cross-terms of $\big( V_{yx} + V_{y\bar{x}} \big)^2$ that do not appear in $\mathcal{N}_\textrm{MF}$; however, these cross-terms vanish as only V_yx or $V_{y\bar{x}}$ can be nonzero for the same y. It follows that $\mathcal{N}_\textrm{MF} = \mathcal{N}$, and:

$$\begin{align} \langle \psi_{k} | \psi_{k, \textrm{MF}} \rangle = \frac{1}{\sqrt{2}} . \end{align} \tag{ C.34 }$$

These results are evaluated numerically for a small system in figure C1. Here, we consider N = 6 qubits split into two fragments with $N_f = 3$. The qubits interact all-to-all, with J_ij drawn from a Gaussian distribution centered at zero with width 1.0. Three Hamiltonians are considered: the full Hamiltonian H, the fragmented Hamiltonian $H^{(f)}$ which neglects all interactions crossing the fragment interface, and finally the mean-field corrected Hamiltonian $H_\textrm{MF}^{(f)}$, which replaces interface interactions by mean-field terms. Sweeping across the transverse field strength h, the minimum energy state is calculated and compared to the ground state of H, $|\psi_g\rangle$. To calculate the minimum energy state of $H_\textrm{MF}^{(f)}$ (which is state dependent and thus cannot be computed by solving a linear eigenvalue problem), we perform gradient descent directly on the state vector, minimizing the energy cost. Additionally, the ground states of H and $H_\textrm{MF}^{(f)}$ are approximately computed using first order perturbation theory.

Zoom In Zoom Out Reset image size

As predicted, the minimum energy state of $H_\textrm{MF}^{(f)}$ has fidelity equal to 0.5 with $|\psi_g\rangle$ for small h. As h is increased beyond $~0.25$, the system enters a regime where first order perturbation theory is no longer sufficient to describe the state. Finally, we note that for small h, the minimum energy state of $H^{(f)}$ has negligible overlap with $|\psi_g\rangle$: when J_ij interactions are totally neglected, it is unlikely that the same $|x^*\rangle$ will minimize the energy.

The role of batch size T is explored in figure D1. Here, we consider the geometric average of the final error after fragmented pre-training, contrasting the results to those of vanilla VQE (which utilizes a combination of random and near identity initialization for brickwork and all-to-all layers, respectively) plotted as dashed lines. On average, even pre-training a single set of fragmented circuits reduces the final error by an order of magnitude or more. The performance gain due to fragmented pre-training grows even larger as T is increased, and notably, the amount of error reduction grows with system size. These results indicate that increasing the batch size T will continue to reduce the error on average. In general, larger batch sizes are necessary to produce the same average error as the system size is increased. This can be understood by the fact that the number of possible ways to partition the system can increase with N, although exact scaling depends heavily on how the fragmentation is constrained. Empirically, the scaling of averaged final error is piece-wise. For small T, the scaling decays according to a power law; the data for T < 10 is fit to a power law and the resulting fit is provided by the dotted lines in figure D1. At roughly T = 10, the averaged error begins to plateau, and increasing T beyond 10 provides diminishing returns.

Zoom In Zoom Out Reset image size