Qubit vitrification and entanglement criticality on a quantum simulator

Abstract

Many elusive quantum phenomena emerge from a quantum system interacting with its classical environment. Quantum simulators enable us to program this interaction by using measurement operations. Measurements generally remove part of the entanglement built between the qubits in a simulator. While in simple cases entanglement may disappear at a constant rate as we measure qubits one by one, the evolution of entanglement under measurements for a given class of quantum states is generally unknown. We show that consecutive measurements of qubits in a simulator can lead to criticality, separating two phases of entanglement. Using up to 48 qubits, we prepare an entangled superposition of ground states to a classical spin model. Progressively measuring the qubits drives the simulator through an observable vitrification point and into a spin glass phase of entanglement. Our findings suggest coupling to a classical environment may drive critical phenomena in more general quantum states.

Introduction

The Born rule, which states that the outcome of a measurement performed on a quantum state is a random variable whose probability distribution is determined by quantum theory, governs information transfer from a quantum system to its classical environment. While the Born rule is at play constantly all around us since all matter is fundamentally quantum, its effects are only evidenced in bulk due to the astronomical number of measurement events that occur at macroscopic lengths and time scales. In contrast, a quantum simulator, a programmable array of qubits, is an otherwise isolated quantum system that we can couple to its environment at will with tailor-made measurements of all or some of the qubits. As such, it allows us to study in detail how the quantum characteristics of a system change as we progressively measure its components.

If we were to pick the state of an ideal quantum simulator uniformly at random from the set of all states accessible to it, we would get a volume-law state: a state whose entropy of entanglement of an extensive subsystem, measured in bits, is proportional to the number of qubits. Measuring a single qubit in such a state removes roughly one bit of entanglement. One therefore expects that the entanglement of a volume-law state should decrease as a simple linear function of the number of measured qubits, yielding a classical unentangled state once we measure all the qubits.

However, this is not always the case: the behaviour of entanglement in a quantum simulator can change dramatically and abruptly as we progressively measure its qubits, exhibiting the phenomenon of criticality¹. Critical behaviour indicates a transition between two distinct phases of entanglement. Previous work revealed entanglement phase transitions in different settings, namely, in random ensembles of quantum circuits in which qubits undergo measurement at a finite rate^{2,3,4,5,6,7,8} and in models with topological order^9,10,11.

Detecting and characterizing such critical behaviour in experiments is challenging though. First, current quantum processors are faulty, limiting the quantum gates we can reliably implement and the number of error-free measurements we can obtain. Second, measuring the entanglement entropy of an arbitrary quantum state is hard. The straightforward approach requires tomographic reconstruction of the quantum state, which is intractable for quantum systems with many components. Finally, although theoretical models for entanglement phase transitions have been introduced in the context of random circuit ensembles^2,3,5,8, these are only approximate in the experimentally relevant limit.

Here, we program two entanglement phases and the criticality between them on a quantum simulator of up to 48 superconducting qubits. We do this by implementing an ensemble of quantum circuits¹ that allow us to reliably generate volume-law states whose entanglement we can deliberately decrease with qubit measurements, experimentally determine the entanglement entropy, and capture the exact dependence of entanglement on measurements using a physical theory¹ (see “Measuring entanglement”). The pertinent physical theory is that of spin vitrification, i.e., the transition to a spin glass phase^12,13. We detect the vitrification point, which agrees with spin glass theory. Our work shows measurements alone can trigger entanglement criticality, suggesting a classical environment could induce critical behaviour in more general quantum states.

Results

Theory and model

Our experimental system is an array of superconducting qubits on a quantum simulator. To drive these qubits through an entanglement phase transition, we first execute quantum circuits that prepare highly entangled states. The circuits and theory that follow come from previous theoretical work¹ on entanglement phase transitions.

Each of our circuits implements a system of R linear equations on L Boolean variables using R + L = N qubits. (In practice, we use fewer qubits to implement our system on hardware. See “Circuit optimization” in the “Methods”.) We can write the system as the matrix equation

$$B{{{{{{{\bf{x}}}}}}}}={{{{{{{\bf{y}}}}}}}}\,{{{{{{\mathrm{mod}}}}}}}\,\,2\,,$$

(1)

where B is an R × L Boolean matrix whose rows represent equations and columns represent variables, as shown in Fig. 1a. If B_ij = 1, then equation i involves variable j. Otherwise, the entry is zero. By setting the elements of B according to some distribution, we get an ensemble of matrices. Each element of the parity vector ${{{{{{{\bf{y}}}}}}}}\in {\left\{0,\, 1\right\}}^{R}$ fixes the parity of an equation and each ${{{{{{{\bf{x}}}}}}}}\in {\left\{0,\, 1\right\}}^{L}$ is a solution to the system for a given y.

Fig. 1: An example of system (Eq. (1)) with L = 6 and R = 4.

a The matrix B of Eq. (1), where each column represents a variable and each row represents an equation. b The first linear equation in B compiled as a circuit. The gate labelled “1” is a CNOT followed by a SWAP, as shown in the inset of (c). We use this circuit design to construct the full quantum circuit for Bx = y. c The quantum circuit built from B. B_ij = 1 corresponds to the two-qubit gate shown in the inset, whereas B_ij = 0 corresponds to a SWAP gate. The output state ${|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}}}$ holds all y that yield solutions to Eq. (1) and the corresponding x for the matrix B. Throughout all panels, the variable labels are in light blue and the parity/equation labels are in dark blue.

To implement the system of Eq. (1) using a quantum circuit, we organize the qubits in the simulator in two registers, as sketched in Fig. 1b,c: a “variable” register consisting of L variable qubits and a “parity” register consisting of R parity qubits. We input variable qubits in the $|\!+\! \rangle$ state and parity qubits in the $\left|0\right\rangle$ state into a circuit from the ensemble defined above. The initial state is therefore ${|\psi \rangle }_{{{{{\rm{in}}}}}}={|\!+\rangle }^{\otimes L}{|0\rangle }^{\otimes R}$. The state ${|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}}}$ at the output of the circuit is an entangled equal superposition of solutions x for each possible y (see “Quantum state and entanglement entropy” in the “Methods” for the exact state). For each y, the set of solutions $\left\{{{{{{{{\bf{x}}}}}}}}\right\}$ is unique. Each parity qubit at the output holds the parity of the variables that appear in the corresponding row of B (see Fig. 1b). Because these variable qubits are in a superposition of classical states, so is the parity qubit: it is 0 or 1 depending on the state of the variable qubits. The variable qubits are thus entangled with the parity qubit and contribute one bit of entanglement.

What is the total entanglement between variable and parity qubits? We quantify this with the entanglement entropy S, which counts the number of bits of entanglement between two parts of a quantum system. To answer our question above, we need to know how many possible vectors y there are in the superposition. For each of the rank(B) linearly independent rows of B, the corresponding component in y can be set to zero or one freely. There are then N_y = 2^rank(B) possible vectors y which give solutions to Eq. (1). This means the entanglement entropy between the variable and parity qubits is $S({|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}}}) \sim {\log }_{2}{N}_{{{{{{{{\bf{y}}}}}}}}}={{{{{{{\rm{rank}}}}}}}}(B)$.

Entanglement phase transition

Each time we measure a qubit in a generic volume-law state, the system loses one bit of entanglement. This happens with every measurement, reducing the number of superposed configurations until just a single classical state remains.

The states ${|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}}}$ behave in a manifestly different manner. To prove this, we start by compiling enough linear equations into our circuits such that an output state is a superposition of N_y = 2^rank(B) = 2^L vectors y. Then, we calculate the entanglement entropy after measuring a subset M of the parity qubits in our system in the computational basis. We label the measurement outcome y_out,M and the partially measured state ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}},M}$. The state ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}},M}$ still contains an equal superposition of solutions to Eq. (1), but the elements of the parity vector y that correspond to the subset M are fixed to the measurement outcome y_out,M. For the same reasoning as in the previous section, there are ${N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}={2}^{{{{{{{{\rm{rank}}}}}}}}({B}_{M})}$ equally probable measurement outcomes y_out,M, where B_M are the rows of B that correspond to the M parity qubits. (In Fig. 1c, measuring the first three parity qubits would mean B_M is the first three rows of B.) Since measuring the output state determines y_out,M, there are ${N}_{{{{{{{{\bf{y}}}}}}}}}/{N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}={2}^{L-{{{{{{{\rm{rank}}}}}}}}({B}_{M})}$ vectors y remaining in the state ${|\psi \rangle }_{{{{{{\rm{out}}}}}},M}$. The entanglement entropy is then $S({|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}},M}) \sim {\log }_{2}({N}_{{{{{{{{\bf{y}}}}}}}}}/{N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}})=L-{{{{{{{\rm{rank}}}}}}}}({B}_{M})$. Therefore, the evolution of the entanglement entropy is given by rank(B_M) as a function of M. To obtain a size-independent control parameter, we define α ≡ ∣M∣/L, the ratio of measured parity qubits to variable qubits.

A unique feature of states ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}}}$ is that we can obtain an exact result for the evolution of the entanglement entropy under measurements through a mapping to a classical spin model¹. In this description, ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}}}$ is an entangled superposition of ground states x to a spin Hamiltonian with couplings defined by y (see “Methods”). We can then use the characteristics of the spin model both to predict the behaviour of the entanglement and, more importantly, to measure it on a quantum simulator, as we discuss in the next section.

To get concrete predictions for our experiments, we now specify the distribution we sample to populate the matrix B and define the ensemble of states ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}}}$. We pick three distinct variables uniformly at random for each equation in Eq. (1), and we ensure there are no repeated equations. With this choice, we get an exact correspondence between our output state and the ground states to an instance of the unfrustrated 3-spin model, both of which are given by the solutions to Eq. (1)^14,15. This model exhibits a phase transition at α_c ≈ 0.918. For α < α_c, our system corresponds to a paramagnet in the 3-spin model. We thus get a paramagnetic phase of entanglement. Here, rank(B_M) = ∣M∣ = Lα. This happens because there are few rows in B_M, which makes it highly probable that they are linearly independent. The entanglement entropy of the quantum system after measuring ∣M∣ parity qubits scales linearly in both L and α: $S({\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}},M}) \sim L\left(1-\alpha \right)$, i.e., the output state obeys a volume law. For α > α_c, the system enters a spin glass phase. The qubits vitrify, turning into an entangled superposition of spin glass ground states. We thus get a glassy phase of entanglement. Now, rank(B_M) < ∣M∣ because there are many rows in B_M and linear independence is lost. The entanglement entropy still scales linearly in L but decreases slower than linearly with increasing α.

We sketch how the measurement of parity qubits collapses the output state in Fig. 2. The two entanglement phases give rise to different behaviours. In the paramagnetic phase, measuring a parity qubit collapses its state to one of two equally probable values (Fig. 2b). This occurs when a parity qubit is independent of the other measured parity qubits, which is the case when B_M has full rank. Measuring the parity qubit halves the number of possible vectors y remaining in the superposition, so the system loses one bit of entanglement entropy. In the spin glass phase, there is a finite probability that a parity qubit has a definite value before measurement (Fig. 2c). This occurs when previous measurement outcomes determine the measurement outcome for the next parity qubit, which begins when rank(B_M) < ∣M∣. In this case, measuring the parity qubit does not change the number of possible vectors y, so the entanglement entropy remains the same.

Fig. 2: Measurement and collapse in the two entanglement phases.

a Prior to measurement, each parity qubit is either in a superposition as in (b), or already collapsed from previous measurements as in (c). In the paramagnetic phase, measurement collapses the parity qubit, decreasing the entanglement entropy by one bit. In the spin glass phase, a finite fraction of the parity qubits are already collapsed from previous measurements, so measuring them does not affect the entanglement. An abrupt change between the two behaviours occurs at ∣M∣ = Lα_c.

Measuring entanglement

The exact correspondence between spin glass physics and entanglement established above and in ref. [1] gives us an efficient way to detect the entanglement phases and criticality on a quantum simulator. In our setup, the entanglement entropy is characterized by the spin glass order parameter ^16,17,18

$$q({B}_{M})=\frac{1}{L}\mathop{\sum }\limits_{i=1}^{L}{\langle {(-1)}^{{x}_{i}}\rangle }^{2},$$

(2)

where 〈…〉 is an average over all the solutions x for Eq. (1) with matrix B_M and a given parity vector y_out,M, and x_i is the i-th variable in x. We derive the identity that links this order parameter to the entanglement entropy in the “Methods”. Therefore, while in most quantum systems quantifying entanglement is intractable, here we have direct access to the entanglement entropy through the spin glass order parameter, which is efficiently measurable (see “Methods”).

The order parameter describes the tendency of the solutions x to take the same value on each variable. It is zero in the paramagnetic phase, which means a given variable does not have correlations across solutions. The outcome of the measurement of a parity qubit is independent of previous parity measurements in this phase. At α = α_c, the variables abruptly become correlated across solutions, and the order parameter jumps to a finite value, eventually saturating to one. This implies solutions of the system are almost identical, differing on only a few variables. The outcome of the measurement of a parity qubit now depends on previous parity measurements.

In the language of physics, each variable can be thought of as one of L spins in a many-body system with ∣M∣ interactions. Then, each basis state in ${\left|\psi \right\rangle }_{{{{{{{{\rm{out}}}}}}}},M}$ of the variable qubits represents a spin configuration. In the paramagnetic phase where there are few interactions, these configurations have no correlation, leading to no order (q = 0). However, in the spin glass phase where ∣M∣ > Lα_c, the interactions induce correlations across the configurations, leading to spin glass order (q > 0).

Using arrays of up to 48 qubits, our experimental results (Fig. 3) clearly reveal the two entanglement phases and the transition between them, and are in agreement with theory. The transition at a critical value α_c becomes more abrupt with increasing system size, exactly as spin glass physics dictates. Finite-size scaling (see “Methods”) of the data using the scaling form $q(\alpha )=f((\alpha -{\alpha }_{c,\exp }){L}^{1/{\nu }_{\exp }})$ yields the experimental values for the critical measurement ratio ${\alpha }_{c,\exp }=0.95\pm 0.06$, which agrees with the theoretical value α_c ≈ 0.918¹⁵, and the critical exponent ${\nu }_{\exp }=2.5\pm 0.5$.

Fig. 3: Experimental results for the order parameter q as a function of measurement ratio α, using the ibm_washington, ibmq_brooklyn and ibm_hanoi quantum processors²⁴.

Each data point represents an average of the order parameter over 900 matrices B_M, except for L = 24, where we average over 50 matrices. Error bars indicate the standard error of the mean. We study system sizes L = 8 (light blue circles), L = 16 (medium blue triangles), and L = 24 (dark blue squares). Dashed curves provide a reference and connect points from a classical simulation of q, where we average over 10, 000 matrices per α. The inset shows the collapsed data around α_c ± 0.5 with the critical exponent ${\nu }_{\exp }$ and critical point ${\alpha }_{{{{{{{{\rm{c}}}}}}}},\exp }$ (light grey vertical line), where the curves in the main plot sharpen with increasing system size. For details on data collection, averaging, and finite-size effects, see “Methods”.

We note that while there are some similarities between the spin glass order we find and those in other work on entanglement criticality^19,20,21, there are a few differences. First, previous work focuses on states and circuits that respect certain symmetries, which play a role in creating a spin glass phase. In contrast, we do not need to impose symmetries. Second, there is an exact relation between the spin glass order parameter in our work and the entanglement entropy, which is not there in other works. Third, previous work focuses on spin glass order as a steady state property of the system, whereas our system goes into a spin glass state immediately after applying our circuit and measuring. Finally, as our results demonstrate, we can observe this spin glass order on existing quantum hardware.

Discussion

Since entanglement is a key resource for quantum computation, precise predictions and experimental verification of its possible behaviours in quantum devices under measurement are sought-after. Our findings demonstrate that partial measurements of quantum states can alone give rise to intricate phenomena related to entanglement. Measurements can force qubits to vitrify, and hence realize the celebrated¹³ spin glass phase of matter inside a quantum processor.

The spin glass quantum states implemented here are a subset of stabilizer states, an important class of states for quantum computation. Moreover, we already know that spin glass entanglement criticality is also present in more general classes of states than the ones studied here¹. It is interesting to ask whether similar physics applies to entanglement in monitored quantum systems at large, giving rise to different types of nonanalyticity for the entanglement entropy.

Methods

Spin Hamiltonian

The output of our circuits provides x and y from Eq. (1), which we can relate to the ground states and couplings of a p-spin model²² (with p = 3). The model consists of L spins, with R interactions encoded by the matrix B (see the example in Fig. 1a). The indices of the nonzero elements in each row a ∈ B correspond to the spins which are part of an interaction. The Hamiltonian is:

$$H(B,\, {{{{{{{\boldsymbol{\sigma }}}}}}}},\, {{{{{{{\bf{J}}}}}}}})=\frac{1}{2}\mathop{\sum}\limits_{a\in B}\left(1-{J}_{a}{\sigma }_{{a}_{1}}{\sigma }_{{a}_{2}}{\sigma }_{{a}_{3}}\right),$$

(3)

where a₁, a₂, a₃ refer to three distinct spins (the indices of the nonzero elements in a), J_a are the couplings of the interaction vector ${{{{{{{\bf{J}}}}}}}}\in {\left\{\pm 1\right\}}^{R}$, and the spins σ_i form the spin vector ${{{{{{{\boldsymbol{\sigma }}}}}}}}\in {\left\{\pm 1\right\}}^{L}$. The ground-state energy for this Hamiltonian is zero, which corresponds to ${J}_{a}{\sigma }_{{a}_{1}}{\sigma }_{{a}_{2}}{\sigma }_{{a}_{3}}=1$ for all a.

Using the mapping ${J}_{a}={(-1)}^{{y}_{a}}$ and ${\sigma }_{i}={(-1)}^{{x}_{i}}$, we see that Eq. (3) is zero whenever ${y}_{a}={x}_{{a}_{1}}+{x}_{{a}_{2}}+{x}_{{a}_{3}}\,{{{{{{\mathrm{mod}}}}}}}\,\,2$ for all a, which is Eq. (1). This lets us express the number of ground states ${{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}$ to Eq. (3) in terms of y and B. Because each of the N_y = 2^rank(B) vectors y has ${{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}$ ground states (out of a possible 2^L), we find ${{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}={2}^{L-{{\mbox{rank}}}(B)}$. The ground-state entropy is ${S}_{{{{{{{{\rm{GS}}}}}}}}}(B)\equiv \log {{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}=\left[L-\,{{\mbox{rank}}}\,(B)\right]\log \!2$ (we take the natural logarithm).

Quantum state and entanglement entropy

After applying the circuit given by B (Fig. 1c) to our input state ${\left|\psi \right\rangle }_{{{{{{{{\rm{in}}}}}}}}}$, we get the following state (see ref. 1 for more details):

$${|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}}}=\frac{1}{\sqrt{{N}_{{{{{{{{\bf{y}}}}}}}}}}}\mathop{\sum}\limits_{{{{{{{{\bf{y}}}}}}}}}|{{{{{{{\bf{y}}}}}}}} \rangle|\{{{{{{{{\bf{x}}}}}}}}\,:\,B{{{{{{{\bf{x}}}}}}}}={{{{{{{\bf{y}}}}}}}}\} \rangle \,.$$

(4)

This is a superposition of solutions $\left\{{{{{{{{\bf{x}}}}}}}}\right\}$ for each of the N_y possible y. We then measure the state of the first ∣M∣ parity qubits to be y_out,M. The resulting state is

$${|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}},M}=\frac{1}{\sqrt{{N}_{{{{{{{{\bf{y}}}}}}}}}/{N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}}}\mathop{\sum}\limits_{\{{{{{{{{\bf{y}}}}}}}}\,:\,{{{{{{{{\bf{y}}}}}}}}}_{|M|}={{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}\}}|{{{{{{{\bf{y}}}}}}}} \rangle|\{{{{{{{{\bf{x}}}}}}}}\,:\,B{{{{{{{\bf{x}}}}}}}}={{{{{{{\bf{y}}}}}}}}\}\rangle \,,$$

(5)

where the first ∣M∣ components of y are y_∣M∣ = y_out,M. The state is still an equal superposition of solutions for each y, but now there are only ${N}_{{{{{{{{\bf{y}}}}}}}}}/{N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}$ terms in the sum, with ${N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}={2}^{{{{{{{{\rm{rank}}}}}}}}({B}_{M})}$. The coefficient ${\lambda }_{{{{{{{{\bf{y}}}}}}}}}=\sqrt{{N}_{{{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}}/{N}_{{{{{{{{\bf{y}}}}}}}}}}$ determines the entanglement entropy for such a state:

$$S({|\psi \rangle }_{{{{{{{{\rm{out}}}}}}}},M})\equiv -\mathop{\sum}\limits_{\{{{{{{{{\bf{y}}}}}}}}\,:\,{{{{{{{{\bf{y}}}}}}}}}_{|M|}={{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M}\}}{\lambda }_{{{{{{{{\bf{y}}}}}}}}}^{2}\log ({\lambda }_{{{{{{{{\bf{y}}}}}}}}}^{2})=\left[\,{{\mbox{rank}}}(B)-{{\mbox{rank}}}\,({B}_{M})\right]\log \!2.$$

(6)

The entanglement entropy coincides with S_GS(B_M)—the ground-state entropy of the classical spin model—when we choose our initial B to satisfy rank(B) = L. In the limit of large system sizes, the expression for the averaged entropy density²³ is:

$$\mathop{\lim }\limits_{L\to \infty }\frac{\left\langle S({|\psi \rangle }_{{{{{{\rm{out}}}}}},M})\right\rangle }{L}=\left[(1-q(\alpha ))(1-\log (1-q(\alpha )))-\alpha (1-{q}^{3}(\alpha ))\right]\log \! 2,$$

(7)

where q(α) is the spin glass order parameter after performing the ensemble average using Eq. (2). Equation (7) establishes the exact correspondence between the entanglement entropy and the spin glass order parameter.

Quantum hardware

We used the ibm_washington, ibmq_brooklyn and ibm_hanoi quantum processors²⁴ for our experiments. We set the repetition delay to 0.00025s. We chose connected lines of qubits to take advantage of the one-dimensional structure of our circuits, while also having low measurement readout error and CNOT error rates at the time of job submission.

Circuit optimization

Because errors dominate the output in current quantum processors, we optimize our circuits to use as few gates as possible. As the CNOT is the native two-qubit gate on the IBM Q processors, we use the CNOT count N_CNOT as our metric (we ignore the L single-qubit Hadamard gates we always need). We build our circuits using the matrix B_M since it generates the solutions we use in Eq. (2) and requires fewer gates to implement than B. We then decompose SWAP gates in our circuits as

$$\,{{\mbox{SWAP}}}(i,\, j)={{\mbox{CNOT}}}(i,\, j)\times {{\mbox{CNOT}}}(j,\, i)\times {{\mbox{CNOT}}}\,(i,\, j),$$

(8)

where i and j are the qubits participating in the gate and CNOT(i, j) means qubit i controls the target qubit j. For the “1” gate in Fig. 1, using Eq. (8) reveals two consecutive CNOT gates with the same control and target, which we remove because they have no overall effect. As such, a one in the matrix requires two CNOTs while a zero requires three.

How many qubits and CNOT gates do we need to build circuits such as in Fig. 1c using B_M? There are ∣M∣ = Lα linear equations, so we require N = L + ∣M∣ = L(1 + α) qubits. To calculate N_CNOT, note that each row of B_M contains p ones and L − p zeros. There are then 2p + 3(L − p) CNOTs per row of B_M, where p = 3 for our model. Summing the CNOT count over all rows, we find ${N}_{{{{{{{{\rm{CNOT}}}}}}}}}({B}_{M})=|M|\left[2p+3(L-p)\right]=3L\left(L-1\right)\alpha$.

By transforming B_M using matrix row operations, we can reduce N and N_CNOT. We begin by putting B_M into row echelon form. Then for each row of the matrix (starting from the second), we find the index of the leading one, and add this row to the rows above it which have a zero at that index. These row additions generate more ones in the matrix, which we prefer because they require less gates than zeros to implement. We call the resulting matrix ${B}_{M^{\prime} }$. Note that solutions to Eq. (1) for a given parity vector remain unchanged under row operations.

For example, consider the following matrix:

$${B}_{M}=\left(\begin{array}{llllll}1&0&1&0&0&1\\ 0&1&0&1&0&1\\ 0&1&1&1&0&0\\ 0&1&1&0&1&0\\ 1&1&0&0&0&1\end{array}\right).$$

(9)

Applying the operations described gives

$${B}_{M^{\prime} }=\left(\begin{array}{llllll}1&1&1&1&1&0\\ &1&1&1&1&0\\ &&1&1&1&1\\ &&&1&1&0\\ &&&&1&0\end{array}\right),$$

(10)

where the omitted entries are zeros.

The form of ${B}_{M^{\prime} }$ helps us save qubits and gates. First, the matrix has size $|M^{\prime}|\times L$, with $|M^{\prime} |={{{{{{{\rm{rank}}}}}}}}({B}_{M})$. The resulting circuit requires N = L + rank(B_M) ≤ L + ∣M∣ qubits, fewer than the circuits built from B_M when rank(B_M) < ∣M∣. Second, notice that in Fig. 1c, there are gates for each entry of the matrix. By interspersing the parity and variable qubits instead of separating them, we can avoid including gates for the zeros to the left of the leading ones in ${B}_{M^{\prime} }$.

Each gate in the primitive circuit (Fig. 1b) exchanges the positions of the qubits it acts upon. The result is that a parity qubit exchanges positions with every variable qubit. However, only “1” gates contribute to the parity we want to measure. Therefore, once a parity qubit encounters all the “1” gates for its linear equation, we measure it in that position. We take advantage of this by altering the primitive circuit: we start the parity qubit at the top, reverse the gate sequence, and invert the control and target of the CNOTs in each “1” gate. Then, the locations of the leading ones in ${B}_{M{\prime} }$ provide the end positions for measuring the parity qubits. For example, the parity qubit for the first row in Eq. (10) exchanges positions with all variable qubits before we measure it. The parity qubit for the second row exchanges positions with variable qubits 6, 5, 4, 3, and 2 before measuring, and so on.

We provide an upper bound for ${N}_{{{{{{{{\rm{CNOT}}}}}}}}}({B}_{M{\prime} })$. We count the entries in ${B}_{M{\prime} }$ to the right of (and including) the main diagonal. We assume the leading ones are all part of the main diagonal. With this assumption, the number of entries to the right of (and including) the leading ones in a P × Q row echelon matrix is:

$$U(P,\, Q)=\mathop{\sum }\limits_{i=1}^{P}\left(Q-[i-1]\right)=P\left(Q-\frac{1}{2}\left(P-1\right)\right).$$

(11)

In our case, P = rank(B_M) and Q = L. There is at least a one per row (else the row would not be a part of ${B}_{M{\prime} }$). Since zeros contribute more to N_CNOT, we take the worst-case scenario where all other entries are zero. This implies there are P ones in the matrix, so there are Z = U(P, Q) − P zeros. Using these results, we get the upper bound:

$${N}_{{{\mbox{CNOT}}}}({B}_{M^{\prime} })\le 3Z+2P=\frac{3}{2}{{{{{{{\rm{rank}}}}}}}}({B}_{M})\left[2L-{{{{{{{\rm{rank}}}}}}}}({B}_{M})+\frac{1}{3}\right]\le \frac{1}{2}L\left(3L+1\right),$$

(12)

where we get the final inequality by maximizing the previous expression with respect to rank(B_M).

Finally, rather than putting a variable qubit in its initial superposition $\left |+\right\rangle=H\left|0\right\rangle$ as an input to the circuit, we apply the Hadamard gate H only when the corresponding variable first participates in a linear equation. (For example, in Eq. (10) variable 6 is first part of an equation in row 3.) Doing so reduces errors from trying to maintain superpositions for too long in current quantum processors. It also simplifies any SWAP gate involving a variable qubit in the state $\left|0\right\rangle$. If we have qubits i and j with the latter in the state $\left|0\right\rangle$, Eq. (8) reduces to

$$\,{{\mbox{SWAP}}}(i,\, \, j){|}_{{{\mbox{}}}j \, {{{{{{{\rm{in}}}}}}}} \, \left|0\right\rangle {{{}}}}={{\mbox{CNOT}}}(i,\ j)\times {{\mbox{CNOT}}}\,(j,\ i).$$

(13)

Our circuit optimization provides a significant savings compared to the circuits built from B_M, which require N_CNOT(B_M) = 3L(L − 1)α gates and $N=L\left(1+\alpha \right)$ qubits. In practice, our largest experiments (L = 24 and α ≳ 1) required an average of ${N}_{{{{{{{{\rm{CNOT}}}}}}}}}({B}_{M^{\prime} })\approx 600$ gates, which is much less than the N_CNOT(B_M) ≳ 1600 gates we would need if we used the B_M matrices instead.

Error mitigation and shot count

In our experiments, we only keep measurement results (shots) x and ${{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$ which satisfy ${B}_{M^{\prime} }{{{{{{{\bf{x}}}}}}}}={{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$. This provides significant error mitigation (see Supplementary Fig. 1) as we increase the number of qubits. For $L=\left\{8,\,16,\,24\right\}$, we took $\left\{10000,\,25000,\,750000\right\}$ shots per sample (see next section) to obtain our data.

Data collection

1. For the desired number of matrix samples:

• (a) Generate a random matrix B as described in the “Entanglement phase transition” section with L columns and $L{\alpha }_{\max }$ rows. Each B is a sample and provides data for $\alpha \in \left(0,\, {\alpha }_{\max }\right]$.

• (b) For each $\alpha \in \left(0,\, {\alpha }_{\max }\right]$:

  • i. Take the submatrix B_M, consisting of the first ∣M∣ = Lα rows of B.

  • ii. Put B_M into row echelon form and perform row operations as described in the “Circuit optimization” section. The result is ${B}_{M^{\prime} }$.

  • iii. Build the circuit corresponding to the matrix ${B}_{M^{\prime} }$ using the techniques described in the “Circuit optimization” section.

  • iv. Execute the circuit on the quantum processor a sufficient number of times. Here, sufficient means measuring several pairs $({{{{{{{\bf{x}}}}}}}},\, {{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} })$ that pass the test in the following step. We always had at least 18 pairs.

  • v. Test measurements $({{{{{{{\bf{x}}}}}}}},\, {{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} })$ by verifying if ${B}_{M^{\prime} }{{{{{{{\bf{x}}}}}}}}={{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$.

  • vi. Save the variable and parity vectors x and ${{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$ that pass the test.

Calculating the order parameter

1. For each $\alpha \in \left(0,\, {\alpha }_{\max }\right]$:

• (a) For each matrix B from the previous section:

  • i. Compute ${B}_{M^{\prime} }$ as in the previous section using B_M with ∣M∣ = Lα.

  • ii. For each saved parity vector ${{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$ associated to ${B}_{M^{\prime} }$:

    • A. Fix a reference solution z that maps solutions from the parity vector ${{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$ to the parity vector 0. We chose our reference to be the solution to ${B}_{M^{\prime} }{{{{{{{\bf{z}}}}}}}}={{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M^{\prime} }$ whose binary form represents the smallest integer. Note that finding a reference is efficient.

    • B. Map the saved solutions x associated with ${{{{{{{{\bf{y}}}}}}}}}_{{{{{{{{\rm{out}}}}}}}},M{\prime} }$ to ${{{{{{{\bf{x}}}}}}}}^{\prime}={{{{{{{\bf{x}}}}}}}}+{{{{{{{\bf{z}}}}}}}}$. Now, ${B}_{M^{\prime} }{{{{{{{\bf{x}}}}}}}}^{\prime}={{{{{{{\bf{0}}}}}}}}$. Remove any duplicates. Call this set $X=\left\{{{{{{{{\bf{x}}}}}}}}^{\prime} \right\}$.

  • iii. Compute $q({B}_{M^{\prime} })$ in Eq. (2) by uniformly sampling $\min \left(24,| X | \right)$ solutions from X, where we determined the number 24 yields a reasonable compromise between accuracy and quantum runtime.

• (b) Compute q(α) by averaging over $q({B}_{M^{\prime} })$ for all ${B}_{M^{\prime} }$.

Following the procedure for each L produces the curves in Fig. 3. To calculate the order parameter in step iii), we want as many solutions ${{{{{{{\bf{x}}}}}}}}^{\prime}$ as possible, but a finite number works, making the order parameter efficient to measure. For the classical simulation of q (the dashed lines in Fig. 3), we uniformly sample $\min \left(24,\, {{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}\right)$ solutions to the equation B_Mx = 0, where ${{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}={2}^{L-{{{{{{{\rm{rank}}}}}}}}({B}_{M})}$ is the total number of solutions. We did this by sampling random linear combinations of the basis vectors forming the null space of B_M. This is also efficient.

We note that the small size of the sample X leads to appreciable artefacts in Fig. 3, such as a deviation of the order parameter from the expected value of zero at small α. Concomitantly, we notice the onset of finite-size effects at values of L and α for which ${\min} \left({24},{\,} {{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}\right) \approx {{{{{{{{\mathcal{N}}}}}}}}}_{{{{{{{{\rm{GS}}}}}}}}}$. The dip of q at small α for L = 8 is such an effect. We nevertheless notice that, for all L and α, theory and experiment match well in Fig. 3, since these effects are present in both.

Finite-size scaling

The objective of finite-size scaling is to take the data in Fig. 3 and try to collapse it onto a common curve by finding suitable critical parameters. We follow the technique of ref. [25] and our previous work¹. In particular, we use the scaling form:

$$q=f\left(\left(\alpha -{\alpha }_{c,\exp }\right){L}^{1/{\nu }_{\exp }}\right),$$

(14)

and we minimize a cost function with the data and its associated error as input to find the critical parameters ${\alpha }_{c,\exp }$ and ${\nu }_{\exp }$.

We store our experimental data as a triple $\left(\alpha,q(\alpha ),e(\alpha )\right)$, where e(α) is the standard error of the mean for the data point. The standard error of the mean for n samples is:

$$e(\alpha )=\sqrt{\mathop{\sum }\limits_{i=1}^{n}\frac{{\left[{q}_{i}(\alpha )-\bar{q}(\alpha )\right]}^{2}}{n(n-1)}},$$

(15)

where q_i(α) is the order parameter for a given matrix B and α, and $\bar{q}(\alpha )$ is the mean of q_i(α) over i.

We then transform the triple according to the scaling form:

$$\left({t}_{i},\, {g}_{i},\, {e}_{i}\right)=\left(\left[\alpha -{\alpha }_{c,\, \exp }\right]{L}^{1/{\nu }_{\exp }},\, q(\alpha ),\, e(\alpha )\right).$$

(16)

We sort these triples by their t-values and then compute the cost function:

$$C({\alpha }_{c,\, \exp },\, {\nu }_{\exp })=\frac{1}{T-2}\mathop{\sum }\limits_{i=2}^{T-1}w\left({t}_{i},\, {g}_{i},\, {e}_{i}|{t}_{i-1},\, {g}_{i-1},\, {e}_{i-1},\, {t}_{i+1},\, {g}_{i+1},\, {e}_{i+1}\right),$$

(17)

with T being the number of data points. The quantity in the summation is:

$$w\left({t}_{i},\ {g}_{i},\ {e}_{i}|{t}_{i-1},\ {g}_{i-1},\ {e}_{i-1},\ {t}_{i+1},\ {g}_{i+1},\ {e}_{i+1}\right)={\left(\frac{{g}_{i}-g}{{{\Delta }}\left({g}_{i}-g\right)}\right)}^{2},$$

(18)

$$g=\frac{\left({t}_{i+1}-{t}_{i}\right){g}_{i-1}-\left({t}_{i-1}-{t}_{i}\right){g}_{i+1}}{\left({t}_{i+1}-{t}_{i-1}\right)},$$

(19)

$${\left[{{\Delta }}\left({g}_{i}-g\right)\right]}^{2}={e}_{i}^{2}+{\left(\frac{{t}_{i+1}-{t}_{i}}{{t}_{i+1}-{t}_{i-1}}\right)}^{2}{e}_{i-1}^{2}+{\left(\frac{{t}_{i-1}-{t}_{i}}{{t}_{i+1}-{t}_{i-1}}\right)}^{2}{e}_{i+1}^{2}.$$

(20)

The cost function $C({\alpha }_{c,\, \exp },\, {\nu }_{\exp })$ measures, for each index i, the squared deviation of the point $\left({t}_{i},\, {g}_{i}\right)$ from the linear interpolation ḡ between the points $\left({t}_{i-1},\, {g}_{i-1}\right)$ and $\left({t}_{i+1},\, {g}_{i+1}\right)$ on either side of the sorted sequence. We exclude the first and last points in the sequence since they have no neighbouring points to the left or right, respectively. The uncertainty (Eq. (20)) is a weighted sum of the squared error of the current point $\left({t}_{i},\, {g}_{i}\right)$ and the squared error from the linear interpolation (Eq. (19)). We skip over any three identical t-values in a row in Eq. (17) because of the division by zero in Equations (19) and (20) (though this only happens for isolated values of ${\alpha }_{c,\exp }$). When this happens, we reduce the denominator of the fraction in front of Eq. (17) by the number of skips.

We plot the cost function over a grid of values near the critical parameters from the literature (Supplementary Fig. 2). This allows us to visualize both the minimum and the uncertainty around it. The collapse for finite-size scaling works best when there are finite-size effects, so we restricted our data for the cost function to the region α_c ± 0.5 (the black connector linking the main plot with the inset in Fig. 3).

We chose our grid for the critical parameters to be ${\alpha }_{c,\exp }\in \left[0.85,\, 1.10\right]$, with a step size of 0.001, and ${\nu }_{\exp }\in \left[1.5,4.0\right]$, with a step size of 0.01. We chose a finer step size for ${\alpha }_{c,\exp }$ because we know the critical threshold. We used a larger step size for ${\nu }_{\exp }$ because there is less precision in the literature for ν.

To estimate our uncertainty, we plot a contour at the level $\left(1+r\right){C}_{{{\mbox{min}}}}$, where r is the size of the maximum deviation we allow in the minimum value. We chose r = 0.25, which means we remain uncertain about the minimum for values that are up to 25% larger. Changing r will grow or shrink the contour. We note in Supplementary Fig. 2 that the cost function’s minimum resides roughly in the centre of the contour. To quantify our uncertainty, we compute the width and height of the rectangle circumscribing the contour. Then, we take the uncertainty in ${\alpha }_{c,\exp }$ to be half the width and the uncertainty in ${\nu }_{\exp }$ to be half the height.

This gives us the following experimental values for the critical point and critical exponent:

$${\alpha }_{c,\exp }=0.95\pm 0.06,\,\,\,{\nu }_{\exp }=2.5\pm 0.5.$$

(21)